谱相关密度

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

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

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

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

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

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

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

${\displaystyle X=x\\ x\\ x\\ .\\ .\\ .\\ ],}$

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

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

${\displaystyle X^{\prime }=\otimes <!--\; \otimes \; -->X,}$

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

${\displaystyle W=)FFT(X)\\ FFT(X)\\ .\\ .\\ .\\ ].}$

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

${\displaystyle W^{\prime }=e^{j\mathrm{\omega <!--\; \omega \; -->}L}\otimes <!--\; \otimes \; -->W\\ e^{j\mathrm{\omega <!--\; \omega \; -->}2L}\otimes <!--\; \otimes \; -->W\\ e^{j\mathrm{\omega <!--\; \omega \; -->}3L}\otimes <!--\; \otimes \; -->W\\ .\\ .\\ .\\ ],}$

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

${\displaystyle \mathrm{\omega <!--\; \omega \; -->}=.}$

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

${\displaystyle S=W^{\prime }{W}^{\prime }\ast <!--\; \ast \; -->,}$

其中${\displaystyle {}^{\ast <!--\; \ast \; -->}}$表示复共轭。换言之，若记${\displaystyle W_i=W^{\prime }}$是${\displaystyle (1,N^{\prime })}$ 的矩阵，${\displaystyle S}$可改写为

${\displaystyle S=W_i^H{W}_i,}$

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

${\displaystyle S^{\prime }=FFT(S).}$

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

${\displaystyle A=\frac{j-<!--\; -\; -->k}{N^{\prime }}+(i-<!--\; -\; -->\frac{P}{2})\mathrm{\Delta <!--\; \Delta \; -->}\mathrm{\alpha <!--\; \alpha \; -->}.}$

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

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

${\displaystyle s^{\prime }=\frac{1}{P}\underset{i=0}{\overset{P-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}{S}^{\prime }.}$

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

${\displaystyle s^{\prime }=\frac{1}{P}\underset{i=0}{\overset{P-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}S,}$

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