含噪多頻信號頻率估計算法

陳 鵬， 劉春華， 蘇 欣， 涂亞慶， 趙少美

(1. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000；2. 中國空氣動力研究與發展中心 設備設計及測試技術研究所，四川 綿陽 621000；3. 陸軍勤務學院 軍事物流系，重慶 401311)

多頻信號可理解為多個單頻信號的疊加，其頻率估計廣泛應用于線性系統識別、低頻機械光譜學、電力系統、核磁共振波譜分析以及無損檢測等諸多領域[1-3]。

目前，針對多頻信號的頻率估計算法大致可分為兩大類：時域法和頻域法。時域法如線性預測法(linear prediction, LP)[4]、計算量更加復雜的多重信號分類法(multiple signal classification, MUSIC)[5]、旋轉不變估計法(estimating signal parameter via rotational invariance techniques, ESPRIT)[6]。時域法分辨率高，但所需采樣序列長，導致計算量大，且抗噪性差，不利于實際應用。頻域法主要利用離散傅里葉變換法(discrete Fourier transform，DFT)對信號進行頻譜分析，可利用DSP(digital signal processing)和FPGA(field programmable gate array)等硬件直接實現、計算速度快，抗噪性好，是研究多頻信號頻率估計的主要方向[7]。

單頻信號是多頻信號的一種特例，也是研究重點。文獻[8]通過頻譜兩點插值和迭代計算，實現了頻率高精度估計，后續稱為AM法。文獻[9]利用頻譜三點插值得到頻率估計值，后續稱為Candan法。針對單頻信號，AM法和Candan法分別是迭代類算法和非迭代類算法中綜合性能最好的算法，但處理多頻信號時，二者均受其他非待估計頻率分量頻譜泄漏的影響，頻率估計精度差，不能直接使用。

為抑制頻譜泄漏影響，文獻[10-11]在AM法的基礎上，采用頻譜相減策略實現了頻譜泄漏校正，提高了頻率估計精度，但在信號頻率間隔較近和中高信噪比條件下的頻率估計精度有待提高。

文獻[12]首先用改進的AM法求取多頻信號各頻率分量的頻率粗估計值，然后構造信號中非待估計的所有頻率分量，并通過相減策略以抑制頻譜泄漏的影響，與Ye等研究中的YA法具有相當的估計性能。但該算法減去非待估計頻率分量這一步驟只進行了一次，導致頻率估計精度受頻率粗估計值的影響。特別是在信號頻率間隔較近、頻譜粗估計存在較大偏差時，由于信號中的非待估計頻率分量減去不徹底，使得算法仍受頻譜泄漏的影響，降低了頻率估計精度。

Djukanovi[13]通過濾除信號中非待估計頻率分量的方式來提高頻率估計精度，針對單頻實信號，提出了DFE-1法。首先利用Candan法對信號進行預處理，得到頻率粗估計值，并生成參考信號；然后將信號與參考信號相乘以實現頻率搬移，通過濾除直流分量的形式降低頻譜泄漏的影響；最后采用AM法對濾除了非待估計頻率分量的信號進行迭代計算，提高了信號頻率估計精度。針對多頻信號，在DFE-1法的基礎上進行了推廣，提出了DFE-2法，實現了多頻信號頻率的高精度估計[14]。

DFE法對提高多頻信號的頻率估計精度提供了思路，但存在設計缺陷，影響了頻率估計精度。本文在分析頻譜泄漏影響和DFE法缺陷的基礎上，提出了新的頻率估計算法，并通過仿真實驗驗證了所提算法的有效性。

1 問題描述

根據信號參數是否隨時間變化，采樣信號可分為平穩信號和非平穩信號，本文以平穩信號為模型進行分析。不失一般性，具有M個頻率分量的平穩多頻信號如式(1)所示。

(1)

特別說明：頻率分量個數M可以是已知的，也可以是未知的。在實際應用中，M一般是未知的，可通過廣義阿卡克信息準則和最小長度描述法等方式計算得到。對頻率分量個數M進行估計屬于信號檢測范疇，不屬于信號參數估計領域，是一個獨立的問題。本文直接利用已有算法得到信號頻率分量個數，然后將M作為先驗知識進行處理，后續不再對其進行討論。

特別地，當M=2，且ω1=-ω2，a1=a2，θ1=θ2時，采樣信號稱為單頻實信號。

(2)

單頻實信號是一種特殊的多頻信號，含有正頻率分量和負頻率分量。設計多頻信號頻率估計算法時，對單頻實信號的頻率估計進行討論，是有必要且具有代表性的。

在對多頻信號進行頻譜分析時，可將其頻率表示為

(3)

式中：km=[ωmN/2π]為第m分量在頻譜中能量最大值點的索引； [·]為取最接近于·的整數；-0.5≤δm≤0.5為第m分量的頻譜偏移量。因此，要得到精確的頻率估計值，需要得到準確的頻譜索引和精確的頻譜偏移量。

在對信號第m分量進行頻率估計時，頻率估計精度受其他頻率分量頻譜泄漏和噪聲的疊加影響。由于DFT法具有很強的抗噪性，噪聲的影響相對較小。為直觀理解頻譜泄漏的影響，在無噪條件下，對多頻信號和單頻信號進行頻譜分析，如圖1所示。

圖1 單頻信號和多頻信號頻譜

可以看出，受其他頻率分量頻譜泄漏的疊加影響，多頻信號待估計頻率分量的頻譜值大于同頻單頻信號的頻譜值，從而影響由頻譜分析得到的頻率估計值，與真實的頻率值存在估計偏差。

2 DFE算法

針對單頻實信號，Djukanovic提出了DFE-1算法，流程如表1所示。

表1 DFE-1算法流程

在此基礎上，Djukanovic提出了針對多頻信號的頻率估計算法，后續稱為DFE-2法。與DFT-1法相比較，兩個算法的流程一致。DFE-2法沿用了DFE-1法濾除信號中非待估計頻率分量的思路，但在進行頻率粗估計和精估計時，均結合了Candan法和三點周期圖最大化法，沒再使用AM法。

DFE法抑制了非待估計頻率分量頻譜泄漏的影響，但存在設計缺陷：

(1) DFE-1法分別利用非迭代類Candan法和迭代類AM法對單頻實信號進行頻率粗估計和精估計，但將粗估計和精估計兩個步驟完全分開，既增加了計算量也降低了頻率估計精度。特別是在信號頻率較低，即信號正頻率和負頻率相隔很近、頻譜泄漏嚴重時，容易導致負頻率分量濾除不徹底，使得頻率估計精度受頻率粗估計的影響更加嚴重。

(2) 相位是關于頻率的函數，但DFE-1法在進行頻率搬移時，只考慮了信號頻率，沒有考慮相位，降低了頻率估計精度。

(3) 針對多頻信號，DFE-2法分別進行了一次頻率粗估計和頻率精估計，可以理解為進行了迭代計算。但在處理濾除了非待估計頻率分量的單頻信號時，選用了非迭代類Candan法和三點周期圖最大化法，導致頻率估計精度受頻率粗估計的影響，且三點周期圖最大化法計算更加復雜，增加了計算量，降低了算法實時性。

(4) DFE-2法也沒有考慮相位的影響。

3 所提算法

為提高多頻信號的頻率估計精度，利用DFE法的思路，提出了新的參數估計算法，詳細步驟如下：

步驟1對采樣信號進行FFT計算，并求信號頻譜中最大的M個極大值。

(4)

(5)

式中：f(·)為求函數·中最大M個極大值的索引；ki為第i分量的索引，i根據各頻率分量頻譜最大值遞減的順序依次排列，即k1和kM分別為信號頻譜M個極大值中最大值和最小值的索引。

(6)

(7)

式中，|·|和∠·分別為取復數·的模和角度。

步驟3對信號進行頻率搬移，濾除信號中的其他非待估計頻率分量。

分析時，按照信號頻譜能量最大到最小的順序進行分析，即從k1依次分析到kM。

首先構造其他非待估計頻率分量的參考信號

(8)

式中，m根據非待估計分量頻譜最大值遞減的順序依次排列。

其次將參考信號與多頻信號相乘，將非待估計中的第m分量(m=1)搬移到0頻附近，得到搬移信號。

ym(n)=x(n)rm(n)

(9)

然后將搬移信號中0頻附近的信號能量視為直流分量，利用式(10)濾除直流分量，并將信號搬移回原頻率處，得到抑制了第m分量的降頻信號。

(10)

最后將降頻信號xM-1(n)代入式(9)和式(10)，濾除非待估計頻率分量中的第m分量(m=2)，得到降頻信號xM-2(n)。按照m的取值順序，重復式(9)、式(10)，依次濾除信號中非待估計頻率分量，最終得到只含有第i分量的降頻信號x1(n)。

步驟4采用AM法對降頻信號x1(n)進行分析，在索引ki兩邊插值，間隔為0.5。利用式(11)計算插值點的頻譜值，并根據兩個插值點的頻譜由式(12)計算頻譜偏移量。

(11)

(12)

步驟5按照i的取值順序，循環計算步驟2～步驟4，得到每個頻率分量的參數粗估計值。

步驟6迭代計算步驟2～步驟5，進一步提高各頻率分量的參數估計精度，得到各頻率分量的幅值和初相位估計值，并利用式(3)得到各頻率分量的頻率估計值。

綜上分析，算法的具體流程如表2所示。

表2 算法流程

此外，針對單頻實信號，根據表2的算法流程，即可得到濾除了負頻率分量頻譜泄漏影響的頻率估計值。

與DFE法的不同處在于：

(1) 不單獨區分信號頻率粗估計和精估計，均在迭代中進行計算，既可以更加有效地濾除非待估計頻率分量的影響，提高頻率估計精度，也可以降低計算量，提升算法實時性。

(2) 濾除非待估計頻率分量時，考慮了相位的影響，有利于更加徹底地濾除非待估計頻率分量，從而提高頻率估計精度。

所提算法和DFE法、以及其他現有優秀算法的性能將在仿真驗證中進行對比分析，檢驗所提算法的優越性。

4 仿真驗證

為檢驗所提算法的有效性，利用MATLAB軟件在不同條件下，對多頻信號和單頻實信號進行頻率估計，且主要對多頻信號進行分析。為降低計算時的隨機誤差，每組仿真進行2 000次蒙特卡羅實驗。為方便分析，將估計結果轉換為均方誤差(mean square errors, MSEs)，并用對數表示。

(13)

式中，L式蒙特卡羅實驗次數。

實驗時，以含有3個頻率分量的信號為例進行頻率估計，并與YA法、AK法、DFE-1法、DFE-2法、Ye法[15]以及克拉美羅下限(Cramer-Rao lower bound，CRLB)[16]進行對比分析。設采樣信號x(n)=1.5ej(3.1ωsn+θ1)+1.4ej((4.7+Δk)ωsn+θ2)+1.2ej((8.3+Δk)ωsn+θ3)+z(n)，長度為128，索引間隔Δk以1為步長從1增加到32，初相位θ1，θ2和θ3獨立隨機取值。

(1) 不同迭代次數

所提算法屬于迭代類算法，首先分析算法在不同迭代次數下的頻率估計性能。仿真時，設SNR為30 dB，迭代次數為1，2，3和4，結果如圖2所示。

圖2 不同迭代次數的頻率估計結果

經由1次或2次迭代時，算法整體估計效果較差，但隨著頻率間隔增大而逐漸變好，且2次迭代的頻率估計精度高于1次迭代的頻率估計精度。當頻率間隔較大(Δk≥5)時，3次和4次迭代具有相當的估計性能，信號各頻率分量頻譜間相互泄漏對所提算法的頻率估計精度影響非常小，頻率估計結果靠近CRLB。當頻率間隔相距非常近時，4次迭代效果更好，考慮到算法的估計精度，特別是在信號頻率相隔較近、頻譜泄漏嚴重時頻率估計精度，后續實驗均采用4次迭代計算。

(2) 無噪聲

由前文分析可知，利用頻域法對多頻信號進行頻率估計時，受頻譜泄漏和噪聲的影響。因此，在無噪聲環境下，可以檢驗各算法對頻譜泄漏的抑制能力，仿真結果如圖3所示。

圖3 無噪聲條件下的頻率估計結果

在不同頻譜泄漏程度下，AK法的頻率估計精度優于YA法。當頻率相隔較近、頻譜泄漏嚴重時，AK法優于DFE-2法；當頻率間的間隔較大、頻譜泄漏影響減小時，DFE-2法的效果優于AK法。相比其他算法，本文算法采用頻率搬移策略濾除了信號中其他非待估計頻率分量，抑制了多頻信號中其他頻率頻譜間的泄漏影響，具有更強的頻譜泄漏抑制能力，明顯強于其他幾種算法，達到了算法的設計目的。

(3) 不同信噪比a

為說明算法在不同信噪比下的估計性能，設x(n)=1.5ej(3.1ωsn+θ1)+1.35ej(7.2ωsn+θ2)+1.2ej(10.3ωsn+θ3)+z(n)，在SNR以1 dB為步長從0增加到30 dB的條件下進行了仿真實驗，結果如圖4所示。

圖4 不同信噪比條件下的多頻信號頻率估計結果

當信號頻率間隔小且各頻率分量的能量大小相近時，算法受噪聲影響大，在SNR<5 dB時的頻率估計精度很差。在5 dB15 dB、即頻率估計精度受頻譜泄漏影響更嚴重時，各算法均出現了頻率估計偏差，逐漸偏離CRLB，但相比于YA法、AK法和DFE-2法，本文算法的頻譜泄漏抑制能力更強，頻率估計精度最高，最接近CRLB，優于其他幾種算法。

(4) 不同信噪比b

圖5 不同信噪比條件下的單頻信號頻率估計結果

當信號頻率較低，即信號正頻率和負頻率分量相隔較近、頻譜泄漏嚴重時，本文算法優勢明顯。當SNR<5 dB時，Ye法、DFE-1法和本文算法具有相當的估計精度，均靠近CRLB。當SNR>5 dB時，Ye法和DFE-1法的估計精度逐漸降低、逐漸偏離CRLB，且DFE-1法優于Ye法。當-5 dB

5 結 論

本文算法通過信號預處理、構造參考信號、頻率搬移、濾除非待估計頻率分量等方式抑制了多頻信號中非待估計頻率分量頻譜泄漏的影響，并經迭代計算得到了各頻率分量精確的頻率、幅值和初相位估計值。

仿真實驗結果表明，所提算法有效地抑制了多頻信號中頻譜泄漏的影響，具有更高的頻率估計精度，頻率估計值的均方誤差比DFE法和其他現有優秀算法的頻率估計結果更靠近CRLB。