基于兩階段加窗插值的多音信號頻率估計算法

柏果，程郁凡，唐萬斌

(電子科技大學通信抗干擾技術國家級重點實驗室 成都 611731)

多音信號的頻率估計是許多領域的基本問題，包括雷達[1]和無線通信[2]。如在正交頻分復用(orthogonal frequency division multiplexing,OFDM)通信系統的載波頻偏估計[2]和多音干擾消除[3]中會遇到多音信號的頻率估計。現有的用于多音信號頻率估計算法主要分為兩類：時域估計算法和頻域估計算法。

時域估計算法通常會需要矩陣求逆、特征值分解、奇異值分解[4]，因此其計算復雜度為O(LN3)，其中L為單音分量的數量，N為采樣點數。

頻域估計算法主要基于離散傅里葉變換(discrete fourier transform,DFT)和一些插值形式，因此他們的計算復雜度一般為O(LNlog2N)，這遠小于時域估計算法，因此頻域估計算法得到了更多的關注和應用。針對單音信號，已有很多插值器被提出[5]，但在處理多音信號時，頻譜泄露造成的分量間互干擾會降低這些插值器的估計精度。為了降低分量間互干擾，許多基于加窗的插值器被提出。文獻[6]提出基于賴夫?文森特類I 窗的插值器。文獻[7]提出漢寧窗，文獻[8]提出余弦窗。文獻[9]提出了一種適用于經典窗函數的補零插值器，然而由于其中的三次方程，該插值器的估計性能較差。文獻[10]提出了一種高階多項式插值算法來補償特定于窗函數的估計器偏差，支持任意窗函數。基于文獻[10]，文獻[11]提出了一個新的插值器，該插值器僅根據窗函數調整偏差補償因子，因此其復雜度更低。雖然這些非矩形窗可以減少分量間的互干擾，但也會造成一定的功率損失。文獻[12]提出了牛頓化正交匹配追蹤(newtonalized orthogonal matching pursuit,NOMP)算法，該算法通過牛頓迭代法和反饋提高頻率估計的精度，然而其計算復雜度較高。文獻[13]提出了一種迭代估計算法，該算法采用了A&M 插值器[14]。文獻[15]提出了一種兩階段估計算法，該算法采用了拋物線插值器[16]來克服分量間互干擾，然而當分量間頻率間隔較小時，其估計性能仍然較差。

為了降低分量間的相互干擾，提高頻率估計性能，本文提出了一種基于兩階段加窗插值(twostage windowed interpolation,TSWI)的多音信號頻率估計算法，該算法包含一個新的支持任意窗函數的插值器。在第一個階段，通過選用可以降低頻譜泄露的窗函數獲得初始估計值；在第二個階段，由于分量間的相互干擾可以通過初始估計值消除，通過選用矩形窗函數彌補第一階段非矩形窗帶來的信噪比(signal-to-noise ratio,SNR)損失并得到精估計值。數值仿真結果表明，本文算法具有比現有算法更優的估計性能。

1 系統模型

含噪聲的多音信號可以表示為：

式中，n=0,1···,N?1，N為樣點數；L為分量數；Al、fl和 φl為第l個分量sl[n]的復幅度、頻率和初相，fs為采樣率；w[n]為加性高斯白噪聲(additive white gaussian noise,AWGN)，w[n]～CN(0,2σ2)。

為了獲得接收信號頻譜的更多細節，可以利用帶補零的DFT 變換，即：

式中，k=0,1···,N′?1；N′=N+P0，P0為補零長度。

將式(1)代入式(2)中，可得：

根據式(4)可得：

通過分析式(7)，可得：

2 基于TSWI 的多音頻率估計算法

2.1 新加窗插值器

基于式(10)，本文提出一種加窗插值器，其推導過程如下。

為了簡化分析，假設L=1。加窗的接收信號可以表示為：

式中，h[n]為窗函數。類似地，[n]的頻域信號為：

根據式(10)，本文提出一種假設，即存在一個系數CN′，使得式(14)為真。

根據式(13)可知：

根據泰勒展開可知，f(x+δl)可展開為：

根據式(15)、式(16)和式(18)可知，式(14)中的分子可以展開為：

類似地，式(14)中的分母可以展開為：

因為a0、a1、b0和b1都是實數，所以有：

將式(20)和式(23)代入式(26)可得：

根據式(27)可知，本小節所提出的假設是真，且系數CN′可以表示為：

特別地，當窗函數為矩形窗時，有CN′=1。

基于式(27)，本文提出了一種加窗插值器，即：

2.2 TSWI 的第一階段

不失一般性地，假設|A1|≥|A2|≥···≥|AL|。在這個階段，l被依次設置為1,2,···,L，當前l?1個分量被估計后，可以根據式(31)將這些分量從接收信號中消除以降低分量間互干擾：

根據文獻[16]可知，kl可估計為：

為了降低分量間互干擾和高階項O()對頻率估計精度的影響，可以通過將第l個分量的頻率移至零附近以減小高階項，并對接收信號采用能降低頻譜泄露的非矩形窗，即：

根據上一節的分析可知，fl可被重新估計為：

此外，Al可估計為：

在這一階段的最后，可得到所有參數的估計值。

2.3 TSWI 的第二階段

雖然第一階段的非矩形窗可以降低頻譜泄露，但是也會造成一定的SNR 損失，因此第一階段的頻率估計值的精度可以被進一步提升。

在這一階段中，l仍然被依次設置為1,2,···,L。因為在第一階段所有分量的參數估計值都已得到，所以分量間的互干擾可以被消除為：

第二階段采用矩形窗，并可將CN′=1、h[n]=1、=和代入式(35)～式(38)更新，然后將和代入式(39)更新。

在第二階段的最后，所有參數估計值都可更新。基于TSWI 的多音信號頻率估計算法的詳細步驟如下一節所示。注意，為了降低計算復雜度，補零長度P0首先設置為?N，此時DFT 可以通過快速傅里葉變換(fast Fourier transform,FFT)實現，然后為了提高估計精度，P0設置為M，M為一個遠大于N的整數。

3 性能分析

根據基于TSWI 的多音信號頻率估計算法可知，頻率估計性能主要由決定。

根據式(40)消除其他分量，得到[n]；

當SNR 較高時，認為?fl/fs≈0，S簡化為：

同理，式(41)可以簡化為：

將式(49)代入式(50)可得：

根據式(46)可知：

式中，SNRl=|Al|2/(2σ2)。

由于第二階段采用矩形窗，式(53)可以簡化為：

4 數值仿真

本節仿真了不同多音信號頻率估計算法，并對其性能進行了對比分析，仿真參數如下：L=2，|A1|=|A2|=1，f1～U[0,fs)，f2=f1+α，α為頻率間隔，∠A1～U[?π,π)，∠A2～U[?π,π)，M=3N，N=256，沿用文獻[13]算法的迭代次數為2。

圖1 給出了TSWI 第一階段采用不同窗函數時的兩個分量的平均均方MSE(root MSE，RMSE)性能，其中α=2.4fs/N，一個DFT 柵格為fs/N，SNR為所有分量的信噪比，其與SNRl的關系可以表示為：

從圖1 中可以看出，采用非矩形窗時的頻率估計性能優于采用矩形窗時的頻率估計性能，這是因為這些非矩形窗可以降低頻率泄露，減少分量間的互干擾。此外，采用漢明窗時TSWI 具有最優的頻率估計性能，當SNR ≥?5 dB時，其RMSE 仿真曲線與理論RMSE 曲線基本重合，且與CRLB 也基本重合，在之后的仿真中，TSWI 第一階段都采用漢明窗。

圖1 TSWI 第一階段采用不同窗函數時的兩分量平均RMSE

圖2 展示了α=2.4fs/N時不同多音信號頻率估計算法的平均RMSE，由于頻率間隔較小，分量間互干擾較大，文獻[11]算法、文獻[13]算法和拋物線插值算法都出現估計誤差平層。與這些算法不同，TSWI 未出現估計誤差平層，當S NR ≥?5 dB時，其RMSE 基本達到了CRLB，這是因為TSWI 在第一階段通過漢明窗抑制了頻譜泄露，并在第二階段通過矩形窗彌補了第一階段漢明窗造成的SNR 損失。

圖2 α=2.4 fs/N時不同多音頻率估計算法的平均RMSE

在圖3 中，α被設置為31.9fs/N，由于頻率間隔較大，所有算法都未出現頻率估計誤差平層，然而，TSWI 的估計性能仍然優于其他算法。

圖3 α=31.9 fs/N時不同多音頻率估計算法的平均RMSE

圖4 對比了不同頻率間隔下各多音信號頻率估計算法的平均RMSE，其中SNR=30dB。可以發現，TSWI 具有最優的頻率估計性能，當α>2fs/N時，TSWI 的RMSE 基本達到了CRLB，其他算法都需要更大的頻率間隔才能不再受分量間互干擾的影響。

圖4 SNR=30dB 時不同多音頻率估計算法的平均RMSE

圖5 將L增大到4，各分量的頻率分別為5.1fs/N、7.2fs/N、16.3fs/N和35.8fs/N，各分量的幅度都為1，且初相都服從區間為 ?π～π的均勻分布。由于L增大到4，所有算法都出現了頻率估計平層，然而TSWI 在SNR>30dB 時才出現誤差平層，其他算法在SNR>20dB 時就出現了誤差平層，此外，TSWI 的誤差平層遠小于其他算法。

圖5 L=4 時不同多音頻率估計算法的平均RMSE

5 結束語

本文提出了一種基于TSWI 的多音信號頻率估計算法，本算法包含一個新的插值器，其插值器支持任意窗函數，通過在第一階段選擇可以抑制頻譜泄露的非矩形窗并在第二階段選擇矩形窗，在降低分量間互干擾的同時不損失SNR。仿真結果表明，本文算法具有比現有算法更優的頻率估計性能。