基于時域平滑FFT累加算法的循環譜實現方法

宋鵬，張麟兮，張曼，唐勇

（1.西北工業大學電子信息學院，陜西西安710129；2.湖北省京山縣第一高級中學湖北荊門431800）

無人機數據鏈所采用的通信信號都是有調制的信號，調制信號的解調多數情況下需要恢復載波信息，載波恢復性能對信號解調有著很大的影響。當信號的信噪比較低時，一般用于載波提取的鎖相環性能嚴重下降或失鎖，通常以環路信噪比6 dB作為環路的失鎖門限[1]。為了在低信噪比條件下，保證無人機順利的完成預定任務，有效載波頻率提取成為無人機數據鏈的關鍵技術之一。

自上世紀80年代，Gardner.W.A.等人提出了循環譜[2-3]這一概念后，循環平穩信號處理技術已在通信、雷達、聲納等信號處理方面得到了廣泛的應用，如：擴頻通信信號的檢測與識別，雷達信號檢測與參數估計以及陣列信號處理等領域。循環譜密度函數將通常的功率譜定義域從頻率軸推廣到了頻譜頻率——循環頻率雙頻率平面，且對于平穩噪聲而言，在非零循環頻率處循環譜密度函數為零，即在非零循環頻率處不呈現譜相關。因此，本文使用循環譜技術，可以有效的擺脫平穩噪聲對參數估計的影響，實現低信噪比條件下，對信號的載頻估計。

1 循環譜相關理論

1.1 循環自相關函數

通常把統計特性呈周期或多周期平穩變化的信號稱為循環平穩或周期平穩信號。設x（t）為循環平穩信號，T為循環周期，循環自相關函數（τ）可以表示為[3]：

其中，τ為時間間隔。我們把Rαx（τ）≠0的頻率α稱為信號x（t）的循環頻率。從式（1）可以看出，循環自相關函數實際是在不同循環頻率上，對輸入信號的相關函數乘上不同的循環系數e-j2παt，使得信號的自相關函數在不同循環頻率上實現相關累積，因此輸入信號的循環譜自相關函數在某些循環頻率點處出現譜峰。

1.2 譜相關密度函數

平穩信號的自相關函數和功率譜密度函數是一對傅里葉變換對，同樣，循環平穩信號的循環自相關函數和譜相關密度函數也是一對傅里葉變換對[4]。信號x（t）的循環自相關函數Rαx（τ）的傅里葉變換稱為循環譜密度函數或譜相關密度函數，其中，XT（f）=即時域信號x（u）的傅里葉變換。

當α=0時，式（1）和式（2）退化為通常的自相關函數和功率譜密度函數，由此可見，循環譜自相關函數Rαx（τ）和循環譜密度函數Sαx（f）是傳統自相關函數和功率譜密度函數Sx（f）的推廣。由于大多數無人機采用BPSK調制信號，本文將以BPSK信號為例，實現循環譜對載頻的估計。BPSK信號的循環譜表達式為：

φ0為初始相位，Tb為碼元寬度，fc為載波頻率。

BPSK信號的循環譜在α=±2fc+N/Tb、α=N/Tb、（N為整數）處出現譜峰，其中α=0及α=±2fc的譜峰值相等且最大。

1.3 譜相關密度函數的估計

通過時域或頻域平滑的方法，可以得到對譜相關密度函數的估計。由于對信號的分析是在有限時間間隔Δt內進行，因此，譜相關密度函數是一個估計值。下面以時域平滑的方法對譜相關函數進行估計[5]。式（2）可重新寫為

其中：

式中

Δt為信號分析的時間間隔，TW是短時傅里葉加窗長度，且為滑動短時傅里葉變換。如圖1所示。在信號分析的有限時間間隔Δt內，頻率分辨率Δt由短時傅里葉加窗長度TW決定，即Δf=1/TW。L表示在每次短時傅里葉變換中的重疊部分。為了避免頻譜混疊和泄露，L的取值范圍[6]應滿足：L≤TW/4。

此外，該算法中還涉及到一個參數M=Δt/TW，即Grenander不確定條件[6]，M的取值應遠遠大于1。該不確定條件意味著信號分析的有限時間間隔Δt應遠遠大于短時傅里葉加窗長度。循環譜的循環頻率分辨率由信號分析的時間間隔Δt所決定，即Δα=1/Δt。

圖1 使用短時傅里葉變換對譜密度函數的估計Fig.1 Using short-time Fourier transform to estimate the spectral density function

2 時域平滑FFT累加算法

由于循環譜的計算是在頻譜頻率——循環頻率雙頻率平面上進行的，因此循環譜的計算量非常大。為了計算的有效性以及硬件實現的響應速度，有效的循環譜實現算法變的不可或缺。循環譜的計算有兩種方法，一種是循環自相關函數的傅里葉變換的方法，一種是信號譜相關法，采用譜相關的數值計算法又可以分為時域平滑FFT累加算法和頻域平滑FFT累加算法[5]。本節將詳細介紹循環譜的時域平滑FFT累加算法。

我們將式（4）和式（5）改寫為

式中

XN′（n，k）為時域信號x（n）在經過數據加窗函數處理后的離散傅里葉變換，w（n）為數據加窗函數，N′為加窗后截取的信號長度，N為數據的總長度。該算法的實現步驟流程如圖2所示。

圖2 時域平滑FFT累加算法步驟流程圖Fig.2 The flow chart of time-smoothing FFT accumulation method

由圖2可以看出，時域平滑FFT累加算法主要由三部分組成。首選，采用加窗函數對輸入的時域信號進行加窗，并對每一個加窗后的截短序列，做N′點的快速傅里葉變換；然后將XN′（n，k）分為上支路和下支路兩個支路分量，分別乘以循環系數e-jπαn/N′和ejπαn/N′，并將下支路信號取共軛，與上支路信號做乘積計算。最后，對所有的加窗截短序列求均值，即可得到信號在數據總長度的循環譜結果。

3 計算機仿真結果

根據前面循環譜的相關理論結果，并采用時域平滑FFT累加算法，對BPSK信號循環譜進行MATLAB仿真。試驗數據的信號參數如下：信號碼元個數為1 024；信息碼元速率為20 Mbps；載波fc=2Rb，即4×107Hz；采樣頻率fs=4 fc；頻率分辨率Δf=fs/64；Grenander不確定條件M=128，滿足遠遠大于1的條件；循環頻率分辨率為Δα=Δf/128；每次短時傅里葉變換中的重疊部分L=16。為了體現循環譜進行調制信號的參數估計可以擺脫平穩噪聲的影響，我們設計了兩組對比試驗，試驗一是理想無噪聲污染信號，試驗二為信噪比是-10 dB的加噪信號。噪聲選用高斯白噪聲。仿真結果分別為圖3和圖4所示。

圖3 無噪聲時，BPSK信號循環譜仿真結果圖Fig.3 The simulation results of BPSK cyclic spectrum when no noise

圖4 信噪比-10 dB，BPSK信號循環譜仿真結果圖Fig.4 The simulation results of BPSK cyclic spectrum when SNR is-10 dB

由前面的分析可知，BPSK信號的循環譜在α=±2fc+N/Tb、α=N/Tb（N為整數）處出現譜峰，其中α=0及α=±2fc的譜峰值相等且最大。且當循環頻率α=0時，譜相關函數即為功率譜密度函數。

由圖3（a）可以看出，在無噪聲影響時，有4個非常明顯的尖峰出現。對f=0做循環頻率的切片圖，可以準確的顯示，在α=±2fc處，即8×107Hz，存在兩個明顯的峰值。對α=0做頻率的切片圖，和傳統的功率譜分析一樣，在載波fc處，存在2個明顯的峰值。

當信噪比惡化到-10 dB時，如圖4（c）所示的循環頻率切片圖上，在正半軸循環頻率α=8×107和負半軸循環頻率α=-8×107處依然有兩個明顯的尖峰存在，且對信號載波與干擾雜波進行信雜比（載波與最高干擾雜波譜線的功率比）分析，經計算可以得到其信雜比為10lg（0.3865/0.2097）=2.656 dB，而在頻率切片圖上，如圖4（d）所示，用傳統的功率譜分析可得信雜比為10lg（1/0.7113）=1.479 dB。可以看出在信噪比為-10 dB時，使用循環頻率對信號載波的檢測比功率譜檢測，信雜比提高了1.177 dB。

從頻譜帶寬范圍來看，采用循環譜估計方法得到的載波頻率非常精確，是一根兩倍的載波譜線。而傳統功率譜估計方法會出現一定的帶寬范圍，這對由多普勒引起的載波頻譜偏差，將出現一定的估計誤差。

4 結論

本文從循環譜的基本理論開始論述，通過使用時域平滑FFT累加算法，實現了對BPSK信號的計算機仿真驗證，證明了在低信噪比條件下，循環譜依然可以有效的檢測信號的載波頻率。因此，將循環譜檢測技術應用到無人機數據鏈低信噪比信號的檢測是一種非常有效的新方法。

[1] Roland E.Best.鎖相環設計、仿真與應用[M].5版.李永明等譯.北京：清華大學出版社，2007.

[2] Gardner W A.Measurement of spectral correlation[J].IEEE Trans.ASSP，1986，ASSP-34（10）:1111-1123.

[3] Gardner W A.The spectral correlation theory of cyclostationary time-series[J].Signal Processing，1986，11（1）:13-36.

[4] Gardner W A，Brown W A，Chen C K.Spectral correlation of modulated signal:part1-analog modulation[J].IEEE Transactions on Communication，1987，35（6）:584-594.

[5] Phillip E P.Detecting and classifying low probability of intercept radar second edition[M].Artech House，2009：513-530.

[6] Lima，A.F.，Jr.，Analysis of low probability of intercept radar signals using cyclostationary processing[D].Naval Postgraduate School Master’s thesis，2002.