谱相关密度

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

定义

时间序列的循环自相关函数${\textstyle x(t)}$计算如下：

$$R_{x}^{\alpha }(\tau )=\int _{-\infty }^{\infty }x\left(t-{\frac {\tau }{2}}\right)x^{*}\left(t+{\frac {\tau }{2}}\right)e^{-i2\pi \alpha t}\,dt$$

其中 (*) 表示复数的共轭。根据Wiener-Khinchin 定理[有疑问，需讨论]，谱相关密度为：

$$S_{x}^{\alpha }(f)=\int _{-\infty }^{\infty }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'$的数据帧，并将其组合成矩阵形式，记为$X$。

$X={\begin{bmatrix}x[0:N']\\x[L:L+N']\\x[2L:2L+N']\\x[3L:3L+N']\\.\\.\\.\end{bmatrix}},$

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

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

$X'={\begin{bmatrix}a(N')\\a(N')\\a(N')\\.\\.\\.\\\end{bmatrix}}\otimes X,$

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

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

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

$W'={\begin{bmatrix}W[0,:]\\e^{j\omega L}\otimes W[1,:]\\e^{j\omega 2L}\otimes W[2,:]\\e^{j\omega 3L}\otimes W[3,:]\\.\\.\\.\\\end{bmatrix}},$

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

$\omega ={\bigg [}-\pi ,...,\pi -{\frac {2\pi }{N'}}{\bigg ]}.$

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

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

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

$S[i,:,:]=W_{i}^{H}W_{i},$

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

$S'[:,i,j]=FFT(S[:,i,j]).$

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

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

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

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

跳过第二次 FFT 直接估计 SCD

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

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

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

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

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