谱相关密度

Template:Orphan谱相关密度 (spectral correlation density , SCD)，有时也称为循环谱密度(cyclic spectral density)或频谱相关函数(spectral correlation function)，是描述时间序列的所有频移版本对的交叉频谱密度的函数。谱相关密度仅适用于周期平稳过程，或称为循环平稳过程，普通平稳过程不具备谱相关性。 谱相关被广泛用于信号检测和信号分类。 谱相关密度与每个双线性时频分布密切相关，但不被认为是 Cohen 类分布。

定义

时间序列的循环自相关函数$x(t)$计算如下：$$R_x^{\alpha }(\tau )=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t-\frac{\tau }{2})x^{*}(t+\frac{\tau }{2})e^{-i2\pi \alpha t}dt$$其中 (*) 表示复数的共轭。根据Wiener-Khinchin 定理[有疑问，需讨论]，谱相关密度为：$$S_x^{\alpha }(f)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{R}_x^{\alpha }(\tau )e^{-i2\pi f\tau }d\tau$$

估计方法

对数字信号而言，SCD 可按任意频率和时间分辨率进行估计。由于直接计算SCD具有较高的计算复杂性，为满足信号实时分析的需求，有几类较为有效的信号谱相关估计方法被提出。

目前常用的算法是 FFT 累加法 (FFT Accumulation Method, FAM) 和带状谱相关法 (Strip-Spectral Correlation Algorithm)， 近日，又有一种新的快速频谱相关 (fast-spectral-correlation, FSC) 算法被提出 。

FFT累加法（FAM）

在本节中，我们将介绍实际在计算机上估计SCD的方法。如使用MATLAB或Python中的NumPy库，以下步骤的实现将相当简单。

FFT累加法 (FAM) 是一种计算 SCD 的数字方法。它的输入是一组 IQ 样本矩阵，输出是复值图像（或者说是一复值矩阵），即目标 SCD。FAM输入的信号、或说是 IQ 样本矩阵 $x$，应为复值张量的形式，或者是尺寸为$(N,)$的多维数组的形式 ，其中数组中的每个元素都是一个 IQ 样本点。

FAM的第一步，是将所输入信号$x$分为多个相互重叠且长度为$N^{\prime }$的数据帧，并将其组合成矩阵形式，记为$X$。

$X=\left[\begin{array}{c}x[0:N^{\prime }]\\ x[L:L+N^{\prime }]\\ x[2L:2L+N^{\prime }]\\ x[3L:3L+N^{\prime }]\\ .\\ .\\ .\end{array}\right],$

其中，$L$两数据帧间起始位置相距的长度。为实现重叠，应有$L<N^{\prime }$ 。$X$是形状为$(P,N^{\prime })$ 的张量， $P$取决于$x$能够容纳多少帧 。

随后，将一形状为$(N^{\prime },)$的窗函数$a(N^{\prime })$ ，应用于$X$的每一行 （如汉明窗等），得到$X^{\prime }$。

$X^{\prime }=\left[\begin{array}{c}a(N^{\prime })\\ a(N^{\prime })\\ a(N^{\prime })\\ .\\ .\\ .\end{array}\right]\otimes X,$

其中$\otimes$是逐元素乘法，也就是将矩阵中的每个元素分别与对应位置的窗函数相乘。接下来，要对中的每一行进行 FFT ，得到$W$。

$W=\left[\begin{array}{c}FFT(X[0,:])\\ FFT(X[1,:])\\ FFT(X[2,:])\\ .\\ .\\ .\end{array}\right].$

$W$就是通常称为瀑布图或频谱图的矩阵。 FAM 的下一步是校正FFT后数据帧的相位延迟。

$W^{\prime }=\left[\begin{array}{c}W[0,:]\\ e^{j\omega L}\otimes W[1,:]\\ e^{j\omega 2L}\otimes W[2,:]\\ e^{j\omega 3L}\otimes W[3,:]\\ .\\ .\\ .\end{array}\right],$

其中$\omega$对应于 FFT 结果中的每个数字频率，是形状为$(N^{\prime },)$张量。

$\omega =[-\pi ,...,\pi -\frac{2\pi }{N^{\prime }}].$

随后，通过求经 FFT 后结果的自相关，得到形状为$(P,N^{\prime },N^{\prime })$张量$S$ 。

$S[i,j,k]=W^{\prime }[i,j]W^{\prime }[i,k{]}^{*},$

其中${*}^{}$表示复共轭。换言之，若记$W_i=W^{\prime }[i,:]$是$(1,N^{\prime })$ 的矩阵，$S$可改写为

$S[i,:,:]=W_i^H{W}_i,$

其中 H 表示矩阵的Hermitian （共轭转置）矩阵。接下来的一步，是将 $S$ 沿着第一维分别进行 FFT。

$S^{\prime }[:,i,j]=FFT(S[:,i,j]).$

$S^{\prime }$是一个包含完整 SCD 信息的三维张量，但我们的目标是构建形状为$(P{N}^{\prime },N^{\prime })$的二维张量，即矩阵或着图像的形式，张量的两个维度分别对应特定频率$f$和循环频率$\alpha$。$S^{\prime }$中所有$\alpha$的值可以通过张量$A$的到，而所有的频率值$f$则记录在张量$F$中。这里的 $f\in [-.5,+.5)$和$\alpha \in [-1,+1)$是归一化频率。

$F[i,j,k]=\frac{j+k}{2N^{\prime }}-.5.$

$A[i,j,k]=\frac{j-k}{N^{\prime }}+(i-\frac{P}{2})\mathrm{\Delta }\alpha .$

上式中，$\mathrm{\Delta }\alpha =1/N$ 。至此，SCD 可以退化为一个二位的图像或矩阵$s$，$S^{\prime }$中的$(f,\alpha )$对都可以赋为0，有效值可以通过$A$和$F$获取 。

跳过第二次 FFT 直接估计 SCD

完整计算一次 SCD 具有相当大的复杂度，复杂度的主要来源是第二轮 FFT。幸运的是，从$S^{\prime }$估计$s^{\prime }$ SCD 的计算公式为

$s^{\prime }=\frac{1}{P}\sum _{i=0}^{P-1}{S}^{\prime }[i,:,:].$

为了更小的计算复杂度，我们可以通过下式，直接从$S$计算$s^{\prime }$，因为在 FFT 前或后计算FFT中所有数值的均值是等效的。

$s^{\prime }=\frac{1}{P}\sum _{i=0}^{P-1}S[i,:,:],$

需要注意的是，$s^{\prime }$将看起来像真正 SCD 的 $s$ 旋转45 度的版本 。