譜相關密度

譜相關密度 (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 度的版本 。