查看“︁谱相关密度”︁的源代码

{{Orphan|time=2023-03-16T00:12:27+00:00}}'''谱相关密度''' (spectral correlation density , SCD)，有时也称为'''循环谱密度'''(cyclic spectral density)或'''频谱相关函数'''(spectral correlation function)，是描述时间序列的所有频移版本对的[[谱密度|交叉频谱密度]]的函数。谱相关密度仅适用于周期平稳过程，或称为[[循环平稳过程]]，普通平稳过程不具备谱相关性。 <ref>{{Cite journal |last=Gardner |first=W.A. |date=1986-10-01 |title=Measurement of spectral correlation |journal=IEEE Transactions on Acoustics, Speech, and Signal Processing |volume=34 |issue=5 |page=1111–1123 |doi=10.1109/TASSP.1986.1164951 |issn=0096-3518}}</ref>谱相关被广泛用于[[信号检测理论|信号检测]]和[[信号情报|信号分类]]。 <ref>{{Cite journal |last=Yoo |first=Do-Sik |last2=Lim |first2=Jongtae |last3=Kang |first3=Min-Hong |date=2014-12-01 |title=ATSC digital television signal detection with spectral correlation density |journal=Journal of Communications and Networks |volume=16 |issue=6 |page=600–612 |doi=10.1109/JCN.2014.000106 |issn=1229-2370 |s2cid=757095}}</ref> <ref>{{Cite journal |last=Hong |first=S. |last2=Like |first2=E. |last3=Wu |first3=Zhiqiang |last4=Tekin |first4=C. |date=2010-01-01 |title=Multi-User Signal Classification via Spectral Correlation |journal=2010 7th IEEE Consumer Communications and Networking Conference (CCNC) |page=1–5 |doi=10.1109/CCNC.2010.5421830 |isbn=978-1-4244-5175-3 |s2cid=17126519}}</ref>谱相关密度与每个[[科恩系列分佈|双线性时频分布]]密切相关，但不被认为是 Cohen 类分布。

== 定义 ==
时间序列的[[模稜函數|循环自相关]]函数<math display="inline">x(t)</math>计算如下：<math display="block">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</math>其中 (*) 表示复数的共轭。根据[[维纳－辛钦定理|Wiener-Khinchin 定理]][有疑问，需讨论]，谱相关密度为：<math display="block">S_x^\alpha(f) = \int_{-\infty}^\infty R_x^\alpha(\tau) e^{-i2\pi f\tau} \, d\tau</math>

== 估计方法 ==
<!-- 檔案不存在 [[File:Spectral_Correlation_Density_of_BPSK_and_QPSK.png|thumb|常见通信信号的 SCD 估计]] ，可從英文維基百科取得 -->
对数字信号而言，SCD 可按任意频率和时间分辨率进行估计。由于直接计算SCD具有较高的计算复杂性，为满足信号实时分析的需求，有几类较为有效的信号谱相关估计方法被提出。

目前常用的算法是 FFT 累加法 (FFT Accumulation Method, FAM) 和带状谱相关法 (Strip-Spectral Correlation Algorithm)， <ref>{{Cite journal |last=Borghesani |first=P. |last2=Antoni |first2=J. |date=October 2018 |title=A faster algorithm for the calculation of the fast spectral correlation |journal=Mechanical Systems and Signal Processing |volume=111 |page=113–118 |bibcode=2018MSSP..111..113B |doi=10.1016/j.ymssp.2018.03.059 |issn=0888-3270 |s2cid=125098069}}</ref>近日，又有一种新的快速频谱相关 (fast-spectral-correlation, FSC) 算法<ref>{{Cite journal |last=Roberts |first=R.S. |last2=Brown |first2=W.A. |last3=Loomis |first3=H.H. |date=1991-04-01 |title=Computationally efficient algorithms for cyclic spectral analysis |journal=IEEE Signal Processing Magazine |volume=8 |issue=2 |page=38–49 |bibcode=1991ISPM....8...38R |doi=10.1109/79.81008 |issn=1053-5888 |s2cid=1763992}}</ref>被提出 。

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

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

FAM的第一步，是将所输入信号<math>x</math>分为多个相互重叠且长度为<math>N'</math>的数据帧，并将其组合成矩阵形式，记为<math>X</math>。

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

</math>

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

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

<math>

X' = \begin{bmatrix}

a(N') \\

a(N') \\

a(N') \\

. \\

. \\

. \\

\end{bmatrix} \otimes X
,
</math>

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

</math>。

<math>W = \begin{bmatrix}

FFT(X[0,:]) \\

FFT(X[1,:]) \\

FFT(X[2,:]) \\

. \\

. \\

. \\

\end{bmatrix}
.

</math>

<math>W

</math>就是通常称为瀑布图或[[时频谱|频谱图]]的矩阵。 FAM 的下一步是校正FFT后数据帧的相位延迟。

<math>  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}
    ,</math>

其中<math>\omega

</math>对应于 FFT 结果中的每个数字频率，是形状为<math>(N',)

</math>张量。

<math>\omega = \bigg[-\pi,...,\pi-\frac{2\pi}{N'}\bigg]
.</math>

随后，通过求经 FFT 后结果的自相关，得到形状为<math>(P,N',N')

</math>张量<math>S

</math> 。

<math>S[i,j,k] = W'[i,j]W'[i,k]^*,</math>

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

<math>S[i,:,:] = W_i^HW_i
,</math>

其中 H 表示矩阵的[[埃尔米特矩阵|Hermitian]] （共轭转置）矩阵。接下来的一步，是将 <math>S</math> 沿着第一维分别进行 FFT。

<math>S'[:,i,j] = FFT(S[:,i,j]).</math>

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

<math>F[i,j,k] = \frac{j+k}{2N'} -.5
.
</math>

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

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

=== 跳过第二次 FFT 直接估计 SCD ===
完整计算一次 SCD 具有相当大的复杂度，复杂度的主要来源是第二轮 FFT。幸运的是，从<math>S'</math>估计<math>s'</math> SCD 的计算公式为

<math>s' = \frac{1}{P}\sum_{i=0}^{P-1}S'[i,:,:].</math>

为了更小的计算复杂度，我们可以通过下式，直接从<math>S</math>计算<math>s'

</math>，因为在 FFT 前或后计算FFT中所有数值的均值是等效的。

<math>s' = \frac{1}{P}\sum_{i=0}^{P-1}S[i,:,:],</math>

需要注意的是，<math>s'

</math>将看起来像真正 SCD 的 <math>s

</math> 旋转45 度的版本 。

== 参考文献 ==
<references group="" responsive="1"></references>

*

[[Category:時頻分析]]
[[Category:信号处理]]