基于時域平滑F(xiàn)FT累加算法的循環(huán)譜實現(xiàn)方法

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

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

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

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

1 循環(huán)譜相關(guān)理論

1.1 循環(huán)自相關(guān)函數(shù)

通常把統(tǒng)計特性呈周期或多周期平穩(wěn)變化的信號稱為循環(huán)平穩(wěn)或周期平穩(wěn)信號。設(shè)x（t）為循環(huán)平穩(wěn)信號，T為循環(huán)周期，循環(huán)自相關(guān)函數(shù)（τ）可以表示為[3]：

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

1.2 譜相關(guān)密度函數(shù)

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

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

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

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

1.3 譜相關(guān)密度函數(shù)的估計

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

其中：

式中

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

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

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

2 時域平滑F(xiàn)FT累加算法

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

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

式中

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

圖2 時域平滑F(xiàn)FT累加算法步驟流程圖Fig.2 The flow chart of time-smoothing FFT accumulation method

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

3 計算機仿真結(jié)果

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

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

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

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

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

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

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

4 結(jié)論

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

[1] Roland E.Best.鎖相環(huán)設(shè)計、仿真與應(yīng)用[M].5版.李永明等譯.北京：清華大學(xué)出版社，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.