基于傅里葉系數(shù)插值的短波信道迭代頻偏估計*

游行遠(yuǎn) 楊平 徐彬彬

（1.哈爾濱工程大學(xué) 信息與通信工程學(xué)院， 黑龍江 哈爾濱150001；2.武漢船舶通信研究所， 湖北 武漢430079）

在短波突發(fā)通信系統(tǒng)中，發(fā)射機與接收機兩端電臺振蕩器不穩(wěn)定［1］，導(dǎo)致發(fā)送端與接收端載波之間存在偏差.短波傳輸環(huán)境中信道衰落導(dǎo)致的多普勒頻移，在經(jīng)過相干解調(diào)下變頻之后接收到的基帶信號存在頻率偏差.現(xiàn)有短波調(diào)解器的波形格式參照了北約STANAG4285、4529、5066 標(biāo)準(zhǔn)，一般通過在波形格式起始位置或數(shù)據(jù)塊結(jié)構(gòu)之間插入已知的8PSK 符號序列（訓(xùn)練序列）來進(jìn)行信號同步捕獲、信道參數(shù)估計、自適應(yīng)均衡.同時波形負(fù)荷的數(shù)據(jù)部分會根據(jù)信息速率的需求采用不同的調(diào)制方式，如QPSK、8PSK、16QAM、64QAM 等.

頻率偏移會直接導(dǎo)致接收信號中有用的信號成分功率衰減，特別對PSK、QAM 調(diào)制來說，頻偏會導(dǎo)致數(shù)據(jù)符號相位旋轉(zhuǎn)與擴散，影響數(shù)據(jù)信息正確解調(diào)，增大系統(tǒng)誤碼率［2］.因此，需要在接收端對頻偏進(jìn)行有效的估計，以修正其對接收信號的影響，提高接收機的性能.

目前國內(nèi)外已經(jīng)提出了不少頻率估計方法，時域方法中基于自相關(guān)的方法一直廣受關(guān)注，其利用信號的自相關(guān)相位來提取頻率.其中，文獻(xiàn)［3］中通過將不同延遲的自相關(guān)相位進(jìn)行加權(quán)平均來提取頻率，但算法的頻率估計范圍受限于參與計算的最大自相關(guān)延遲.文獻(xiàn)［4］中利用自相關(guān)相位的差分來實現(xiàn)頻率估計，避免了頻率受限問題，但低信噪比時性能下降，計算復(fù)雜度增加.在文獻(xiàn)［4］算法的基礎(chǔ)之上，文獻(xiàn)［5］中通過對自相關(guān)函數(shù)進(jìn)行歸一化來改進(jìn)算法，減少了算法求自相關(guān)相位的次數(shù).頻域方法大多是基于離散傅里葉變換（DFT）的方法，此類方法物理意義明確，計算量小，因而得到了廣泛的應(yīng)用［6］.現(xiàn)有基于DFT 的頻偏估計方法主要分為粗估計與精估計兩個步驟，粗估計利用DFT 對信號頻譜峰值進(jìn)行定位，精估計是目前此類方法的研究熱點.如文獻(xiàn)［7-9］中利用DFT 結(jié)果峰值譜線與相鄰譜線進(jìn)行一次插值來得到估計結(jié)果.此類插值方法簡單易行，但性能相對較差.文獻(xiàn)［10-11］中提出了利用DFT 結(jié)果進(jìn)行二分搜索的方法，但需要在粗估計時通過補零的方式（即增加DFT 點數(shù)）來獲取更精細(xì)的頻率分辨率［12］，導(dǎo)致計算復(fù)雜度大幅提升.文獻(xiàn)［13］中利用DFT 峰值譜線±1/2 刻度的傅里葉系數(shù)進(jìn)行迭代，算法理論均方誤差（MSE）為Cramer-Rao下界（CRLB）的1.0147 倍.

為適應(yīng)不同速率調(diào)解器的波形格式，保證不同長度訓(xùn)練序列的頻偏估計精度，文中提出了一種基于傅里葉系數(shù)插值的迭代頻偏估計算法，首先進(jìn)行非補零DFT 的粗估計（即DFT 點數(shù)與訓(xùn)練序列長度保持一致），同時計算信號的信噪比，然后根據(jù)信噪比信息、訓(xùn)練序列長度、對插值算法的誤差分析來確定頻偏搜索區(qū)間，最后在搜索區(qū)間內(nèi)進(jìn)行迭代，并利用拋物線插值獲取更精確的結(jié)果.

1 系統(tǒng)模型

假設(shè)長度為L 的8PSK 符號訓(xùn)練序列為ck（k=1，2，…，L），理想同步條件下，經(jīng)過匹配濾波后，接收到的采樣信號可以表示為

式中：A 為信號幅度；Ts為符號采樣間隔；f0為頻率偏移；φ 為相位偏移，一般由信道時延與采樣偏差產(chǎn)生；n（k）為高斯白噪聲（均值為0，方差為σ2）.接收信號與本地訓(xùn)練序列逐點共軛相乘，可得

式中，（·）*表示求共軛，w（k）=n（k）c*k.不失一般性，設(shè)ckc*k=1，同時將頻偏進(jìn)行歸一化，得到F=f0Ts（以下不做特殊描述時均指歸一化頻偏），式（2）可簡化輸出信號為

式中，w（k）與n（k）有相同的統(tǒng)計特性，服從高斯分布.因此，在短波突發(fā)通信中的頻偏估計可以等效于高斯白噪聲條件下復(fù)正弦信號的頻率估計［1］.

在粗估計步驟中，對式（3）進(jìn)行N 點DFT 變換（N≥L），可得

式中，W（n）為噪聲的傅里葉系數(shù)（服從均值為0、方差為Nσ2的高斯分布）［14］，DFT 的歸一化分辨率為1/N.在實際系統(tǒng)中，真實頻偏一般表示為

式中：m 為整數(shù)部分，即DFT 峰值譜線的位置；δ∈［-1/2，1/2］，為小數(shù)部分，即真實頻偏相對于DFT 峰值譜線的偏差.文中的目標(biāo)是低信噪比條件下對δ 的精確估計.

同時利用接收信號與本地訓(xùn)練序列之間的相關(guān)性，可得

由式（6）可求得接收信號的信噪比:

2 算法原理

2.1 搜索區(qū)間的確定

根據(jù)頻偏粗估計DFT 結(jié)果，其峰值譜線與相鄰譜線可表示為

式中，p=-1，0，1，分別對應(yīng)譜線的位置.為了方便描述，定義Yp=Y（m+p），Wp=W（m+p）.

當(dāng)N?p 時，

將式（9）代入式（8），可得

式中，B=Aexp（jφ） [1-exp（j2 δ） ],通過DFT 峰值譜線與相鄰譜線進(jìn)行Jacobsen 插值［8］，得

式中，WN=W-1-W1，WD=2W0-W-1-W1.高斯白噪聲的傅里葉系數(shù)Wp為，而BN 為O（N），其中O（·）表示函數(shù)階符號［13］.因此，當(dāng)N 取較大值時，式（11）可近似于

由于噪聲均值E［Wp］=0，故為δ 的無偏估計.由于噪聲傅里葉系數(shù)的方差var［Wp］=Nσ2，噪聲的虛部方差為var［Im（Wp）］=Nσ2/2，可得插值結(jié)果的近似方差為

式中：Q（·）為高斯Q 函數(shù)；△＞0，為搜索區(qū)間兩端到δ 的距離.

置信度可以反映置信區(qū)間的可靠程度，置信度越大，用以估計實際頻偏值出錯的概率越小.然后在實際的應(yīng)用中，除了置信區(qū)間的可靠以外，還需要保證足夠的精確.文中取置信度為1-a，a=10-5，得△=，則δ 的置信區(qū)間為，將其作為精頻偏估計的搜索區(qū)間，得到初始狀態(tài)為

2.2 精頻偏估計

精頻偏估計主要利用迭代方法在區(qū)間［m0-△0，m0+△0］內(nèi)進(jìn)行搜索.

首先利用搜索區(qū)間的中心點及兩端，求得對應(yīng)的Y（n）（n=m0-△0，m0，m0+△0），根據(jù)這3 點的大小關(guān)系進(jìn)行如下判決:

因此，下一次迭代的初始狀態(tài)為m1=m0+δ0△0，△1=△0/4.

（2）當(dāng)兩端S（m0-△0）或S（m0+△0）為最大值時，真實頻偏屬于區(qū)間［m0-△0，m0-△0/2］或［m0+△0/2，m0+△0］，即，為了方便算法迭代，文中取δ0=-3/4或δ0=3/4，則下一次迭代的初始狀態(tài)為m1=m0±3△0/4，△1=△0/4.

利用上述判決方法進(jìn)行I 次迭代，可得精頻偏估計結(jié)果：

算法的具體步驟如下：

1）粗頻偏估計Y，并確定頻譜峰值m，即同時根據(jù)式（7）計算信噪比.

2）確定搜索區(qū)間.利用式（11）對頻譜峰值譜線與相鄰譜線進(jìn)行插值，得到，并利用插值算法誤差分析的結(jié)果（式（13））來得到插值結(jié)果的近似方差，最后利用△=4.42來確定精頻偏估計的搜索區(qū)間

3）精頻偏估計.利用I 次（i=0，1，...，I-1）迭代插值，根據(jù)式（4）計算并根據(jù)其大小關(guān)系進(jìn)行判決.

本次迭代結(jié)果為mi+1=mi+δi△i，△i+1=△i/4.I 次迭代結(jié)束后，最終結(jié)果為

3 算法性能仿真分析

為進(jìn)一步驗證所提算法的性能，文中對其進(jìn)行仿真分析，仿真中訓(xùn)練序列采用8PSK 調(diào)制，加性高斯白噪聲（AWGN）信道條件，粗頻偏估計采用非補零的DFT，即N=L，通過Nmc=106次蒙特卡洛實驗得到MSE，計算公式如下：

因短波通信中的頻偏估計可以等效為高斯白噪聲條件下復(fù)正弦信號的頻率估計，故可以利用其CRLB 對算法的性能進(jìn)行驗證，CRLB 可表示為［16］

對于不同長度的訓(xùn)練序列，文中算法的MSE 及CRLB 理論線如圖1所示，其中設(shè)定的δ 均勻分布于區(qū)間［-1/2，1/2］.從圖中可知，不同長度的訓(xùn)練序列存在各自的SNR 閾值：L為128、256、512 的訓(xùn)練序列的SNR 閾值分別為-4、-6、-9dB 左右，即粗估計時譜線峰值準(zhǔn)確定位的概率下降導(dǎo)致了低于SNR 閾值時性能急劇惡化，而高于SNR 閾值時性能趨于平穩(wěn)［17］.當(dāng)I=0，即未使用迭代算法僅使用一次插值獲取的結(jié)果時，算法的MSE 無法逼近CRLB，而使用I=1，2 的迭代插值后，算法對不同長度的訓(xùn)練序列均可逼近CRLB.

當(dāng)L=512、δ 取不同值時文中算法的MSE 如圖2所示.當(dāng)δ 為0.05、0.20、0.40 時，SNR 閾值分別為-11、-10、-9dB；當(dāng)δ 較小時，算法能獲得更小的閾值.

圖1 文中算法對不同長度的訓(xùn)練序列的MSEFig.1 MSE of the proposed algorithm for training sequence with different length

圖2 L=512、δ 取不同值時文中算法的MSEFig.2 MSE of the proposed algorithm with different values of δ and L=512

當(dāng)L=512 時文中算法與現(xiàn)有算法的MSE 比較見圖3，其中設(shè)定的δ 均勻分布于區(qū)間［-1/2，1/2］.時域算法集中利用不同延遲的自相關(guān)相位arg［Ryy（m）］來提取頻率，其中arg［·］表示取幅角，Ryy（m）表示為

式中，M 為最大延遲.文獻(xiàn)［3］算法的估計范圍存在局限性，其估計范圍為0≤≤1/M，文中取M=20.文獻(xiàn)［4］算法解決了估計范圍的問題，其M 需小于或等于L/2，為了提高算法性能，文中取M 為128 和256作比較.文獻(xiàn)［5］算法先歸一化Ryy（m）得到，再利用與進(jìn)行差分相乘，并對其結(jié)果進(jìn)行加窗后，將自相關(guān)相位作為估計結(jié)果.仿真結(jié)果表明：在時域類方法中，文獻(xiàn)［3］算法在信噪比ρ∈［-20，10］時皆可對頻率進(jìn)行跟蹤，但MSE 性能較差；文獻(xiàn)［4］算法的MSE 性能在ρ∈［-20，-10］時較差，但隨著SNR 的增加性能逐步改善，當(dāng)M 為128和256 時，分別在SNR 為10、-7 dB 處逼近CRLB.當(dāng)M=256 時文獻(xiàn)［5］算法在SNR 為-5 dB 處逼近CRLB，其性能趨近于文獻(xiàn)［4］算法，但因僅需要一次求相位運算而降低了計算復(fù)雜度.頻域類方法中，文獻(xiàn)［3］算法無法逼近CRLB;文獻(xiàn)［10-11］算法在粗估計使用N=2L 的DFT 時，SNR 閾值為-11 dB，表明適當(dāng)?shù)卦黾覦FT 長度可以提高DFT 頻譜的精細(xì)程度，提高頻譜峰值準(zhǔn)確定位的概率，但會提升計算復(fù)雜度；文獻(xiàn)［13］算法在I=2 時可以逼近CRLB.

圖3 L=512 時不同算法的MSE 比較Fig.3 Comparison of MSE among different algorithms with L=512

當(dāng)L=512、SNR 為0dB 時文中算法與現(xiàn)有算法的MSE 比較見圖4，其中MSE 與CRLB 之比R 為

SNR 為0dB 時CRLB 的理論值為CRLB0.時域類方法在δ 取不同值時性能較穩(wěn)定，而頻域類方法在δ 取不同值時性能會起伏.其中文獻(xiàn)［10-11］算法在使用非補零DFT 時無法逼近CRLB，而在使用N=2L 的DFT 時，通過10 次迭代搜索，針對不同的δ 均可以獲得較好的性能.文獻(xiàn)［13］算法在I=1、δ 越趨近于0時性能越好，而在I=2 時可逼近CRLB;文中算法在I=1，2 時皆可逼近CRLB.

SNR 為0 dB 時文中算法與文獻(xiàn)［13］算法在不同插值次數(shù)（I′）下的MSE 與CRLB 之比R 如表2所示.由于文獻(xiàn)［13］算法未使用粗估計的DFT 結(jié)果進(jìn)行插值，僅利用其對頻偏進(jìn)行定位，故以插值次數(shù)為對比標(biāo)準(zhǔn)，相應(yīng)地文中算法的迭代次數(shù)I=I′-1，文獻(xiàn)［13］算法的I=I′.從表中可知，文中算法在δ=0.0，0.1，0.2 時通過兩次插值即可獲得較穩(wěn)定的性能，得到的R 可達(dá)1.005 左右，而δ 較大時需要3 次插值才能獲得較穩(wěn)定的性能.文獻(xiàn)［13］算法在δ 較小時需要3 次插值才能獲得趨于穩(wěn)定的性能，而δ較大時需要4 次插值才能獲得趨于穩(wěn)定的性能，其R 在1.01～1.02 之間.

圖4 L=512、SNR=0 dB、δ 取不同值時幾種算法的MSE比較Fig.4 Comparison of MSE among several algorithms with L=512，SNR=0dB and different values of δ

計算復(fù)雜度也是評價算法性能的重要指標(biāo)［18］，時域類方法主要集中于自相關(guān)性Ryy（m）的計算，其計算復(fù)雜度與最大延遲有關(guān)，需要M（2L-M-1）/2 次復(fù)數(shù)相乘運算，當(dāng)M=L/2 時，其計算復(fù)雜度為3L2/8，為O（L2）.頻域類方法包括粗頻偏估計與精頻偏估計兩部分，粗頻偏估計主要為N=L 的DFT，需要進(jìn)行（L/2）log2L 次復(fù)數(shù)運算；而精頻偏估計主要利用傅里葉系數(shù)的相互關(guān)系進(jìn)行迭代，主要集中于Y（n）的計算，計算量為CY=L（8CM+6CP），為O（L），CM為復(fù)數(shù)乘法運算，CP為復(fù)數(shù)加法運算，插值部分的計算復(fù)雜度多為O（1），可以忽略.文中算法通常迭代兩次即可獲得較好的性能，單次迭代計算量為3CY.文獻(xiàn)［13］算法需要迭代4 次才能保證算法的性能趨于穩(wěn)定，單次迭代計算量為2CY.文獻(xiàn)［10-11］算法在粗估計時的補零DFT 需要Llog2（2L）次復(fù)數(shù)運算，同時其精頻偏估計需要迭代10 次，單次迭代計算量為CY.

表2 兩種算法的R 值比較Table2 Comparison of R between two algorithms

4 結(jié)論

通過對現(xiàn)有基于DFT 類方法的研究，文中分析了Jacobsen 插值算法的誤差，建立了估計值方差與信噪比及序列長度的一種函數(shù)關(guān)系，據(jù)此來提取估計值的置信區(qū)間并作為精頻偏搜索區(qū)間，之后利用迭代在區(qū)間內(nèi)進(jìn)行插值的精頻偏估計.仿真結(jié)果表明：在確立的置信區(qū)間內(nèi)進(jìn)行迭代精頻偏估計，能有效地提高算法的估計精度，減少迭代次數(shù)；在低信噪比條件下，通過兩次迭代文中算法的MSE 即可逼近CRLB 理論線， 故適用于短波突發(fā)通信的信道環(huán)境.文中算法也適用于高斯白噪聲環(huán)境下對復(fù)正弦信號的頻率估計.

［1］Palmer J，Rice M.Low-complexity frequency estimation usingmultipledisjointpilotblocksinburst-modecommunications［J］.IEEE Transactions on Communications，2011，59（11）：3135-3145.

［2］Pollett T，Moeneclaey M.The effect of carrier frequency offset on the performance of band limited single carrier and OFDM signals ［C］//Proceedings of IEEE Global Telecommunications Conference.London：the Key to Global Prosperity，1996：719-723.

［3］Fitz M P.Further results in the fast estimation of a single frequency ［J］.IEEE Transactions on Communications，1994，42（2/3/4）：862-864.

［4］Mengali U,Morelli M.Data-aided frequency estimation for burst digital transmission ［J］.IEEE Transactions on Communications,1997,45（1）：23-25.

［5］Awoseyila A B,Kasparis C,Evans B G.Improved single frequency estimation with wide acquisition range［J］.Electronics Letters，2008，44（3）：245-247.

［6］曹燕.含噪實信號頻率估計算法研究［D］.廣州：華南理工大學(xué)電子與信息學(xué)院，2012.

［7］Quinn B G.Estimating frequency by interpolation using Fourier coefficients ［J］.IEEE Transactions on Signal Processing，1994，42（5）：1264-1268.

［8］Jacobsen E，Kootsookos P.Fast，accurate frequency estimators ［J］.IEEE Signal Processing Magazine，2007，24（3）：123-125.

［9］CandanC.Amethodforfineresolutionfrequencyestimation from three DFT samples ［J］.IEEE Signal Processing Letters，2011，18（6）：351-354.

［10］Zakharov Y V，Tozer T C.Frequency estimator with dichotomous search of periodogram peak［J］.Electronics Letters，1999，35（19）：1608-1609.

［11］Zakharov Y V，Baronkin V M，Tozer T C.DFT-based frequency estimators with narrow acquisition range ［J］IEE Proceedings：Communications，2001，148（1）：1-7.

［12］Fang L，Duan D，Yang L.A new DFT-based frequency estimator for single-tone complex sinusoidal signals ［C］//Proceedings of Military Communications Conference.Orlando：IEEE，2012：1-6.

［13］Aboutanios E，Mulgrew B.Iterative frequency estimation by interpolation on Fourier coefficients ［J］.IEEE Transactions on Signal Processing，2005，53（4）：1237-1242.

［14］Provencher S.Estimation of complex single-tone parameters in the DFT domain ［J］.IEEE Transactions on Signal Processing，2010，58（7）：3879-3883.

［15］StoerJ，BulirschR.Introductiontonumericalanalysis［M］.New York：Springer，2002：43-49.

［16］韋崗，李道遠(yuǎn)，陳芳炯.基于逆積分方程的正弦波頻率估計［J］.華南理工大學(xué)學(xué)報:自然科學(xué)版，2007，35（10）：147-151.Wei Gang，Li Dao-yuan，Chen Fang-jiong.Sinusoidal frequency estimation based on inversion integral equation［J］.Journal of South China University of Technology：Natural Science Edition，2007，35（10）：147-151.

［17］Aboutanios E.A modified dichotomous search frequency estimator ［J］.IEEE Signal Processing Letters，2004，11（2）：186-188.

［18］Quinn B G.Recent advances in rapid frequency estimation ［J］.Digital Signal Processing，2009，19（6）：942-948.