基于立方Farrow 結(jié)構(gòu)的插值細(xì)搜頻偏估計(jì)算法

陸俊宇， 李 鈺

（華東理工大學(xué)信息科學(xué)與工程學(xué)院, 上海 200237）

高斯最小頻移鍵控（GMSK）是調(diào)制指數(shù)為1/2 和頻率脈沖為高斯脈沖的連續(xù)相位調(diào)制（Continuous Phase Modulation，CPM）信號(hào)，具有包絡(luò)恒定、相位連續(xù)的特性，適用于功率和頻帶資源有限的通信場景[1-3]。GMSK 信號(hào)相干解調(diào)的前提是要求接收機(jī)能夠產(chǎn)生和接收信號(hào)同頻、同相的載波，但在許多具有高動(dòng)態(tài)相對(duì)運(yùn)動(dòng)的通信環(huán)境中，如衛(wèi)星、飛機(jī)、列車等，由相對(duì)運(yùn)動(dòng)引起的多普勒頻移將會(huì)在接收機(jī)下變頻所得的基帶信號(hào)中引入殘留載波頻偏，致使基帶信號(hào)的相位軌跡發(fā)生偏離，最終造成解調(diào)誤碼[4-5]。因此，精確的載波同步是解調(diào)前的重要一環(huán)。載波同步的目的是從接收信號(hào)中提取載波參數(shù)，之后根據(jù)提取的參數(shù)值對(duì)接收信號(hào)進(jìn)行修正。載波參數(shù)的捕獲算法主要分為閉環(huán)和開環(huán)兩種。閉環(huán)算法一般基于鎖相環(huán)（PLL）結(jié)構(gòu)，具有較高的跟蹤精度，缺點(diǎn)是捕獲速率慢、收斂時(shí)間長，不適合突發(fā)信號(hào)載波參數(shù)的捕獲。開環(huán)載波參數(shù)捕獲算法不需要反饋誤差信息，可直接根據(jù)接收信號(hào)的特征來提取載波參數(shù)。開環(huán)算法又分為數(shù)據(jù)輔助（Data-Aided，DA）[6-7]和非數(shù)據(jù)輔助（Non-Data-Aided，NDA）[8-9]。DA 開環(huán)算法采用了添加前導(dǎo)序列的方法來估計(jì)載波參數(shù)，具有捕獲精度高、抗噪性能好的優(yōu)點(diǎn)，但缺點(diǎn)是添加的前導(dǎo)序列需要占用額外的帶寬；NDA 開環(huán)算法不需要添加先驗(yàn)信息，能保持較高的頻譜利用率，但是在低信噪比下的載波參數(shù)捕獲精度較低。

本文在短時(shí)突發(fā)通信[10]的應(yīng)用背景下，采用了DA 開環(huán)載波參數(shù)提取算法，引用了文獻(xiàn)[11]提出的最優(yōu)前導(dǎo)序列結(jié)構(gòu)和克拉美羅界（CRB），在H&P（Hosseini & Perrins）算法[12]的基礎(chǔ)上，通過對(duì)快速傅里葉變換（FFT）所得頻譜進(jìn)行一維搜索來估計(jì)載波頻偏值。受限于柵欄效應(yīng)，僅通過FFT 后直接掃頻的方法只能估計(jì)整數(shù)倍頻偏，存在較大的估計(jì)誤差[13-14]。為了進(jìn)一步提高頻偏估計(jì)的精度，本文提出了一種基于Farrow 結(jié)構(gòu)的立方插值細(xì)搜索算法。在FFT 有限補(bǔ)零的條件下，通過粗搜索確定細(xì)搜索范圍，使用立方Farrow 插值結(jié)構(gòu)計(jì)算細(xì)搜范圍內(nèi)的小數(shù)倍頻點(diǎn)的插值譜線，以逼近真實(shí)的頻偏譜線。通過仿真驗(yàn)證，本文所提出的基于Farrow 結(jié)構(gòu)的立方插值細(xì)搜算法能在H&P 算法的基礎(chǔ)上進(jìn)一步提高載波頻偏的估計(jì)精度，具有良好的估計(jì)性能和工程應(yīng)用價(jià)值。

1 頻偏估計(jì)的實(shí)現(xiàn)框架

頻偏估計(jì)的實(shí)現(xiàn)框架如圖1 所示。框架主要由H&P 算法和插值細(xì)搜（IDS）算法構(gòu)成，其基本的設(shè)計(jì)思路是在H&P 算法的基礎(chǔ)上引入本文提出的插值細(xì)搜算法以獲得更高精度的載波頻偏估計(jì)值。在H&P 算法部分，將含有前導(dǎo)序列部分的基帶數(shù)字信號(hào)輸入預(yù)處理模塊（Pre），得到2 路預(yù)處理的輸出信號(hào)并分別送入2 個(gè)FFT 模塊，計(jì)算這2 路預(yù)處理信號(hào)的頻譜序列，依次取模（Abs）、相加得到2 路預(yù)處理信號(hào)的幅度譜序列。之后將該幅度譜序列進(jìn)行分流，一路送入IDS 算法部分，另一路送入H&P 算法部分的Max 模塊進(jìn)行掃頻找到粗搜索的最高譜線。在IDS 算法部分，插值模塊根據(jù)輸入幅度譜序列確定細(xì)搜區(qū)域，然后在細(xì)搜區(qū)域內(nèi)進(jìn)行固定倍數(shù)的插值，再經(jīng)由Max 模塊掃頻得到最高插值譜線，最后將粗搜索最高譜線和插值最高譜線進(jìn)行比較，選最高譜線所對(duì)應(yīng)的頻點(diǎn)作為載波頻偏的估計(jì)值。

圖1 頻偏估計(jì)的實(shí)現(xiàn)框架Fig.1 Implementation framework of frequency offset estimation

2 H&P 算法的基本原理

H&P 算法是面向短時(shí)突發(fā)通信場景提出的載波頻偏、載波相偏、定時(shí)偏差聯(lián)合同步算法，由于信號(hào)的持續(xù)時(shí)間較短，在接收機(jī)已經(jīng)知道基帶信號(hào)中前導(dǎo)序列位置的條件下，可以假設(shè)在前導(dǎo)序列持續(xù)時(shí)間內(nèi)接收基帶信號(hào)所攜帶的載波頻偏為固定常數(shù)[15]。在加性高斯白噪聲（Additive White Gaussian Noise，AWGN）信道條件下，GMSK 接收信號(hào)經(jīng)過下變頻后的復(fù)基帶接收信號(hào)[16]可以表示為

式中Es為碼元能量；T為碼元周期；w(t) 為復(fù)基帶高斯白噪；誤差參數(shù)有載波頻偏fd、載波相偏 θ 、定時(shí)偏差 ε ； φ 為攜帶碼元信息的信號(hào)相位，其表達(dá)式為

其中，調(diào)制指數(shù)h=1/2 ；前導(dǎo)序列 α={αk} ，k=0,1,···,L0?1，L0為序列的碼元個(gè)數(shù)。最優(yōu)前導(dǎo)序列的結(jié)構(gòu)如圖2 所示，序列主要分為3 段，前L0/4個(gè)碼元和后L0/4的碼元取值為?1，中L0間/2個(gè)碼元取值為1。q(t)為相位響應(yīng)，是頻率脈g沖(t)的積分形式，即

圖2 最優(yōu)前導(dǎo)序列結(jié)構(gòu)Fig.2 Optimal leading sequence structure

式中，L為關(guān)聯(lián)長度；Bb為3 分貝（Decibel，dB）帶寬。

隨即，接收機(jī)可通過對(duì)前導(dǎo)序列持續(xù)時(shí)間內(nèi)的基帶信號(hào)r(t;α) 進(jìn)行傅里葉變換，掃描整個(gè)頻譜找到最高譜線，并以最高譜線所對(duì)應(yīng)的頻點(diǎn)作為載波頻偏的估計(jì)值。最后，將頻偏估計(jì)值代入相應(yīng)的閉合表達(dá)式，依次求出定時(shí)偏差、載波相偏。H&P 算法的載波頻偏、定時(shí)偏差、載波相偏的估計(jì)式依次為

其中，

式中T0表示前導(dǎo)序列的持續(xù)時(shí)間。

因?yàn)槎〞r(shí)偏差和載波相偏是通過將載波頻偏估計(jì)值分別代入它們的閉合表達(dá)式求取的，其估計(jì)精度同樣受到載波頻偏估計(jì)精度的影響，所以本文的關(guān)注重點(diǎn)在于如何在數(shù)字接收機(jī)系統(tǒng)中進(jìn)一步提高載波頻偏的估計(jì)精度。在數(shù)字接收機(jī)系統(tǒng)中，GMSK 接收信號(hào)需要經(jīng)下變頻、模數(shù)轉(zhuǎn)換后變?yōu)閿?shù)字基帶信號(hào)序列r(n;α) 。式（9）和（10）為連續(xù)傅里葉變換的形式，與之對(duì)應(yīng)的是數(shù)字信號(hào)的離散傅里葉變換（Discrete Fourier Transform，DFT），在變換之前可先對(duì)數(shù)字序列拆分后進(jìn)行補(bǔ)零處理，DFT 的輸入信號(hào)為

其中，N為每個(gè)碼元的樣本點(diǎn)數(shù)，P為補(bǔ)零倍數(shù)，表示補(bǔ)零后的序列長度為之前的P倍，則序列r1(n;α)和r2(n;α) 的DFT 形式寫作

實(shí)際上可以使用更高效的FFT 來實(shí)現(xiàn)DFT 的計(jì)算。取 λ1(k) 與 λ2(k) 的模并相加得離散譜線

式中，譜線相鄰頻點(diǎn)間隔 Δf=fs/(PNL0)，fs為采樣頻率。記最高譜線的頻率索引為km，則歸一化頻偏估計(jì)值為vd=kmΔf/fs。

3 基于立方Farrow 結(jié)構(gòu)的插值細(xì)搜算法

3.1 細(xì)搜算法的基本原理

圖3 為細(xì)搜算法的流程圖，算法的工作流程為：

圖3 細(xì)搜算法的流程Fig.3 Flow chart of the detailed search algorithm

(1) 系統(tǒng)對(duì)補(bǔ)零后的序列r1(n;α) 與r2(n;α) 進(jìn)行FFT，然后分別取模，相加得到粗搜索譜線X(k) ；

(2) 對(duì)譜線進(jìn)行第一次的掃頻，即粗搜索，找到最高譜線及其兩側(cè)的次高譜線，記最高譜線的頻率索引為km，次高譜線的索引為km?1 、km+1 ；

(3) 以次高譜線之間的頻率范圍作為細(xì)搜索范圍，分別在索引區(qū)間 [km?1,km] 和 [km,km+1] 進(jìn)行一定倍數(shù)的插值，由此在2 個(gè)離散整數(shù)頻點(diǎn)之間得到小數(shù)倍頻點(diǎn)的插值譜線；

(4) 在插值譜線中搜索最高譜線，即細(xì)搜索，并比較最高插值譜線和譜線X(km) 之間的大小，在兩者中取幅值最大譜線所對(duì)應(yīng)的頻點(diǎn)作為頻偏估計(jì)值。

3.2 立方插值濾波器的Farrow 結(jié)構(gòu)

Farrow 結(jié)構(gòu)是多項(xiàng)式插值算法的數(shù)字化實(shí)現(xiàn)，有助于插值算法的工程應(yīng)用[17]。插值濾波器可實(shí)現(xiàn)小數(shù)倍插值，其基本原理可概括為先將數(shù)字信號(hào)重構(gòu)為模擬信號(hào)，然后重采樣。文獻(xiàn)[18]給出了重采樣后的插值公式，根據(jù)插值細(xì)搜算法的設(shè)計(jì)要求將公式重新表達(dá)為

式中，Ti表示插值間隔；j=0,···,I?1，I為插值倍數(shù)；kb表示插值基點(diǎn)； μb,j表示以kb為起始點(diǎn)的第j個(gè)小數(shù)間隔，即 μb,j=jΔμ ，小數(shù)間隔 Δμ∈(0,1) ；Ts表示采樣周期；i為濾波器索引，取值為I1到I2之間的整數(shù)，其取值范圍決定了計(jì)算插值所需的輸入樣本數(shù)為I2?I1+1 ；h(t) 為插值濾波器的有限階脈沖響應(yīng)，可采用關(guān)于 μb,j的多項(xiàng)式逼近[19]，即

其中，M為濾波器抽頭數(shù)，bl(i) 是獨(dú)立于 μb,j的加權(quán)系數(shù)。將式（17）等號(hào)右側(cè)表示為插值多項(xiàng)式的拉格朗日形式，有

插值計(jì)算公式可重寫為

對(duì)于4 點(diǎn)立方插值濾波器[18]，設(shè)置I1= ?M/2 ，I2=M/2?1，M取=4代入式（18）得三次多項(xiàng)式方程組

將式（20）代入（19）得立方插值計(jì)算表達(dá)式

由此繪制出立方插值的Farrow 結(jié)構(gòu)，如圖4 所示。關(guān)于 μb,j的小數(shù)倍插值的計(jì)算復(fù)雜度，主要體現(xiàn)在延時(shí)次數(shù)、縮放因子數(shù)、加法次數(shù)、乘法次數(shù)這4 個(gè)方面。立方插值的Farrow 結(jié)構(gòu)共有4 列，每列計(jì)算是并行執(zhí)行的，前3 列計(jì)算需要延時(shí)3 個(gè)Ts，第4 列延時(shí)2 個(gè)Ts后等待前面幾列計(jì)算完成。小數(shù)間隔μb,j首先需要等待3 個(gè)Ts后才能和第1 列的計(jì)算結(jié)果相乘，之后 μb,j還需要再做2 次乘法和3 次加法，共6 次計(jì)算，后級(jí)計(jì)算需要等待前級(jí)計(jì)算輸出結(jié)果后才能執(zhí)行，等待的時(shí)間為Ts，故計(jì)算完1 次 μb,j的小數(shù)倍插值需要的延時(shí)數(shù)為8。立方插值Farrow 結(jié)構(gòu)的每列都相當(dāng)于1 個(gè)FIR 濾波器，并且它們擁有公共的抽頭系數(shù)1/6、1/3、1/2，因此縮放因子數(shù)為3，它們可存入查找表中待需要時(shí)再提取，不用重新計(jì)算。此外，通過直接觀察Farrow 結(jié)構(gòu)圖，可知完成1 次插值需要的加法次數(shù)為11，乘法次數(shù)為3?？傊?，計(jì)算1 次關(guān)于 μb,j的小數(shù)倍插值需要8 次延時(shí)、3 類縮放因子、11 次加法、3 次乘法，計(jì)算復(fù)雜度低，易于工程實(shí)現(xiàn)。

圖4 立方插值的Farrow 結(jié)構(gòu)Fig.4 Farrow structure for cubic interpolation

3.3 小數(shù)倍頻插值的基本原理

小數(shù)倍頻插值如圖5 所示，根據(jù)立方插值的Farrow 結(jié)構(gòu)，插值計(jì)算的所需樣點(diǎn)集為{X(kb?1) ，X(kb)，X(，kb+}1)，X(kb所+2)以可將索引j從零開始依次加一遞增，取個(gè)I?1μb,j，采用立方Farrow 結(jié)構(gòu)依次計(jì)算2 個(gè)離散頻點(diǎn)kb和kb+1之間的小數(shù)倍頻點(diǎn)的插值譜線。結(jié)合插值細(xì)搜示意圖（圖6），進(jìn)一步給出基于立方Farrow 結(jié)構(gòu)的IDS 算法實(shí)現(xiàn)小數(shù)倍頻插值并估計(jì)載波頻偏的具體步驟：

圖5 小數(shù)倍頻插值示意圖Fig.5 Schematic diagram of decimal frequency interpolation

圖6 插值細(xì)搜示意圖Fig.6 Schematic diagram of interpolation detailed search

(1) 對(duì)FFT 所得譜線X(k) 進(jìn)行第1 次掃頻，通過初步的粗搜索，找到最高譜線X(km) 以及兩側(cè)的次高譜線X(km?1)和X(km+1) 。

(2) 以次高譜線的頻率索引km?1 作為第1 個(gè)插值基頻點(diǎn)，則k=km?1 ，可以取b=1，將第1 個(gè)插值基點(diǎn)記為k，其小數(shù)間隔表示為 μ1,j。對(duì)于立方Farrow 插值，取離散值X(k1?1)、X(k1)、X(k1+1) 、X(k1+2)，在索引區(qū)間k1[?1，k1]依次計(jì)I?算1個(gè)小數(shù)倍頻率索引，k1+，μ1,1，k1+·μ1,2·k1+·μ1,3，k1+μ1,I?1處的y(k1插+μ1,1)值y(k1+μ，1,2)y(k，1+μ1,3)，...，。

(3) 以最高譜線的頻率索引km作為第2 個(gè)插值基頻點(diǎn)，則k=km，可以取b=2 ，并將第2 個(gè)插值基點(diǎn)記為k2，其小數(shù)間隔表示為u2,j。取離散值X(k2?1) ，X(k2)，X(，k2+，1)在X(k2索+2)引區(qū)間[k2k，2+1]之間依次計(jì)算個(gè)I?1小數(shù)倍頻率索k2+引μ2,1、k2+μ2,2，k2，+μ2,3··k2+·μ2,I?，1處的y插(k2+值μ2,1)，y(k2+μ2,2)，y(k，2+μ2·,3)·y(k·2+μ，2,I?1)。

(4) 留下區(qū)域[k1+μ1,1，···，k1+μ1,I?1，k2+μ2,1，···，k2+μ2,I?1]作為細(xì)搜區(qū)間，然后進(jìn)行第2 次掃頻，通過細(xì)搜索確定最高插值譜線。細(xì)搜的結(jié)果存在2 種情況：若插值y(kb+μb,j) 大于X(km) ，則留下kb+μb,j作為頻偏估計(jì)的索引，頻偏的估計(jì)值為(kb+μb,j)Δf；若X大(km于)任何一個(gè)插值，則留下km作為頻偏估計(jì)的索引，頻偏估計(jì)值為kmΔf。

4 算法的仿真驗(yàn)證

使用Matlab 仿真工具驗(yàn)證本文提出的基于立方Farrow 結(jié)構(gòu)的IDS 算法，仿真的參數(shù)設(shè)置如表1所示。仿真給出了采用了分段拋物線Farrow 結(jié)構(gòu)[18]、立方Farrow 結(jié)構(gòu)的IDS 算法的頻偏估計(jì)均方差曲線；同時(shí)，還和文獻(xiàn)[12]中使用的高斯插值算法進(jìn)行了對(duì)比，分析它們之間的頻偏估計(jì)性能；然后，分析了基于立方Farrow 結(jié)構(gòu)的IDS 算法的插值倍數(shù)I和頻偏估計(jì)精度之間的關(guān)系；最后，通過觀察頻偏修正后的星座圖的相位分布進(jìn)一步驗(yàn)證基于立方Farrow結(jié)構(gòu)的IDS 算法的效果。

表1 仿真參數(shù)設(shè)置Table 1 Simulation parameter setting

不同插值方法和不同補(bǔ)零倍數(shù)對(duì)載波頻偏估計(jì)精度的影響如圖7 所示。在無差值條件下的補(bǔ)零倍數(shù)（P）依次設(shè)置為1、2、4。通過觀察發(fā)現(xiàn)，增大補(bǔ)零倍數(shù)有助于提高頻偏的估計(jì)精度，但是僅通過補(bǔ)零的方式只能給估計(jì)精度帶來較少的信噪比增益。在P=4 的條件下引入IDS 算法，I=32 ：采用了分段拋物線 IDS 算法后，頻偏估計(jì)精度獲得了明顯的信噪比增益，均方差曲線隨信噪比的增大而下降，但是下降幅度卻逐漸減小，9 dB 后趨于平穩(wěn)；不同于分段拋物線插值算法，立方IDS 算法表現(xiàn)出了逼近CRB 的頻偏估計(jì)性能，相比高斯插值算法，頻偏估計(jì)性能僅損失約0.7 dB 的增益。

圖7 補(bǔ)零和插值對(duì)頻偏估計(jì)精度的影響Fig.7 Effect of zero padding and interpolation on the accuracy of frequency offset estimation

圖8 給出了基于立方Farrow 結(jié)構(gòu)的IDS 算法在不同插值倍數(shù)下的頻偏估計(jì)均方差，P=4，I由小到大依次取2、4、8、16、32。如圖8 所示，頻偏估計(jì)精度隨插值倍數(shù)的增大而增大，頻偏估計(jì)精度在插值倍數(shù)2 到4 之間有顯著提升，此后增大插值倍數(shù)對(duì)于頻偏估計(jì)精度的提升作用越來越小，16 倍插值之后的提升開始接近極限。表2 示出了16 倍和32 倍插值條件下本文所提算法的頻偏估計(jì)均方差的仿真值，兩者的頻偏估計(jì)性能非常相近，與CRB 相比，精度損失約1.2 dB。因此，實(shí)際工程中的插值倍數(shù)可取16，有助于進(jìn)一步減少插值倍數(shù)所帶來的計(jì)算量。

表2 頻偏估計(jì)均方差Table 2 Mean square error (MSE) of frequency offset estimates

圖8 不同插值倍數(shù)對(duì)頻偏估計(jì)精度的影響Fig.8 Effect of different interpolation multiples on the accuracy of frequency offset estimation

在16 倍插值條件下，采用基于立方Farrow 結(jié)構(gòu)的IDS 算法對(duì)基帶信號(hào)進(jìn)行頻偏補(bǔ)償。圖9 給出了GMSK 基帶信號(hào)頻偏補(bǔ)償前后的相位星座分布，仿真的信噪比為30 dB，頻偏分布 (fd) 為 0.1%fs， 1%fs。在圖中，紅圈為理論相位分布點(diǎn)，藍(lán)圈為實(shí)際的信號(hào)相位分布。通過觀察發(fā)現(xiàn)，頻偏補(bǔ)償前的信號(hào)相位雜散分布，同時(shí)沿單位圓發(fā)生旋轉(zhuǎn)，偏離理論分布點(diǎn)，而頻偏補(bǔ)償后的信號(hào)相位重新收斂于理論位置。由此可以得出，基于立方Farrow 結(jié)構(gòu)的IDS算法估計(jì)的載波頻偏值能夠有效地糾正GMSK 基帶信號(hào)的相位軌跡。

圖9 基帶信號(hào)頻偏補(bǔ)償前后的相位星座圖Fig.9 Phase constellation before and after frequency offset compensation of baseband signal

5 結(jié) 論

在GMSK 信號(hào)的載波同步問題上，為了進(jìn)一步提高H&P 算法的載波頻偏估計(jì)精度，本文提出了一種基于立方Farrow 結(jié)構(gòu)的插值細(xì)搜算法，在FFT 粗搜索的基礎(chǔ)上進(jìn)一步地細(xì)搜索，使用立方Farrow插值結(jié)構(gòu)在由FFT 粗搜索所確定的細(xì)搜區(qū)域內(nèi)進(jìn)行小數(shù)倍頻的插值，確定最高插值譜線，并將其和粗搜索的最高譜線進(jìn)行比較，取幅值最高譜線的頻點(diǎn)作為載波頻偏的估計(jì)值。通過MATLAB 仿真，在前導(dǎo)序列長度為64 的條件下，首先分析了基于立方Farrow結(jié)構(gòu)的插值細(xì)搜算法的頻偏估計(jì)性能，然后分析了頻偏估計(jì)性能和插值倍數(shù)的關(guān)系，最后再觀察頻偏修正后的相位星座收斂情況。仿真結(jié)果表明，本文所提出的基于立方Farrow 結(jié)構(gòu)的插值細(xì)搜算法能夠在H&P 算法的基礎(chǔ)上有效提升頻偏估計(jì)精度，在插值倍數(shù)為16 的時(shí)候具有逼近克拉美羅界的頻偏估計(jì)性能，信噪比估計(jì)損失僅為1.2 dB，具有良好的性能。