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

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

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

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

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

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

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

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

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

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

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

當(dāng)α=0時(shí)，式（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ù)無人機(jī)采用BPSK調(diào)制信號(hào)，本文將以BPSK信號(hào)為例，實(shí)現(xiàn)循環(huán)譜對(duì)載頻的估計(jì)。BPSK信號(hào)的循環(huán)譜表達(dá)式為：

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

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

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

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

其中：

式中

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

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

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

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

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

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

式中

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

[1] Roland E.Best.鎖相環(huán)設(shè)計(jì)、仿真與應(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.