頻率估計的多段正弦信號快速頻譜融合算法

肖瑋，涂亞慶，沈艷林，蘇丹，張磊

1.后勤工程學院后勤信息工程系，重慶 401311

2.涿州綜合倉庫，河北保定 611730

◎信號處理◎

頻率估計的多段正弦信號快速頻譜融合算法

肖瑋1，涂亞慶1，沈艷林1，蘇丹1，張磊2

1.后勤工程學院后勤信息工程系，重慶 401311

2.涿州綜合倉庫，河北保定 611730

多段正弦信號頻譜融合法（簡稱“原融合算法”）是提高低信噪比條件下正弦信號頻率估計精度的一條有效途徑，具有重要研究意義和應用價值。為滿足雷達、聲納、電子對抗等實時性要求較高的頻率估計應用需求，提出多段正弦信號快速頻譜融合算法。該方法通過設計離散時間傅里葉變換（Discrete Time Fourier Transform，DTFT）快速算法、降維處理加權融合頻譜矩陣和1/3主瓣相關性分析處理等措施來降低算法計算量，提高實時性。重點對上述三項措施的原理進行了闡述與分析。計算量對比和仿真實驗表明，多段正弦信號快速頻譜融合算法在精度損失極小的前提下，能夠大幅降低計算量；在信噪比極低的情況下(SNR≤-13 dB)，其性能略優于原融合算法。

快速頻譜融合；多段正弦信號；離散時間傅里葉變換（DTFT）快速算法；降維；1/3主瓣相關

1 引言

多段正弦信號頻譜融合法（簡稱“原融合算法”）是一類新興的正弦信號頻率估計法。該方法通過對多段信號頻譜進行加權融合實現多段信號的信息積累，能夠近似成倍延長被測信號時寬，從而解決單段信號包含信息量不足、抗噪性差的問題，從源頭上為提高低信噪比條件下正弦信號頻率估計精度創造條件[1-4]。因此其頻率估計精度較高，除應用于頻率估計外，還可用于LFM（Linear Frequency Modulated）信號參數估計[1]、VCO（Voltage Controlled Oscillator）非線性度校正[5]等方面，具有重要的理論研究意義和工程實用價值。

多段正弦信號頻譜融合法的計算量主要集中在多段正弦信號頻譜計算、最優加權融合頻譜生成和頻域相關性分析三部分[4]。為滿足雷達、聲納、電子對抗等實時性要求較高的應用需求[2，6]，本文通過設計離散時間傅里葉變換（Discrete Time Fourier Transform，DTFT）快速算法、降維處理加權融合頻譜矩陣和1/3主瓣相關性分析處理等措施來分別降低上述三個環節的計算量，提出頻率估計的多段正弦信號快速頻譜融合算法（簡稱“快速融合算法”）。

2 算法原理

2.1 DTFT快速算法設計

2.1.1 DTFT快速算法計算量分析

在多段正弦信號頻譜融合法中，利用DTFT計算M (M≥2)段單段信號頻譜Xm(fP)所需的計算量較大[4]。文獻[7]提出了一種滑動DTFT快速算法，該方法的實質是通過不斷增加計算序列的長度來實現指定頻率處傅里葉系數的快速計算，并未對DTFT算法本身進行改進，因此不能適用于多段正弦信號頻譜融合法中序列長度不能增加的情況。

在多數數字信號處理的書籍和文獻中，DTFT是針對離散信號的[8-9]，一般僅考慮n=0，1，…，N-1的有限長序列，其定義如式（1）所示：

根據需要確定被測頻率f的取值范圍fscope= [fmin，fmax]、頻率分辨率Δf、采樣頻率fs和輸出點數K= (fmax-fmin)/Δf+1，代入式（1）得：

其中，k=0，1，…，K-1。

利用布魯斯坦（Bulestein）所提出的式（3）[8-9]為DTFT設計快速算法：

則式（2）可以表示為：

由式（4）可知，計算x(n)的K點DTFT可以通過計算序列g(n)和h(n)的線性卷積與表達式e-jπk2Δf/fs的乘積來實現。由于g(n)為長度為N的序列，h(n)為一無窮長的序列；為計算g(n)和h(n)的線性卷積，只需取h(n)在-N+1≤n≤K-1內的值即可，因此可以把h(n)看成一個長度為L=K+N-1的有限長序列。由于時域計算線性卷積效率低，當循環卷積的循環長度大于或等于L+N-1時，計算循環卷積和線性卷積的結果相同，且計算循環卷積可以利用FFT算法[8-9]，因此可將g(n)和h(n)的線性卷積轉換為循環卷積。又因為計算K點DTFT只需要輸出0，1，…，K-1共K個點的卷積值，即使后面N-1個點發生混疊也不會影響前面的值，因此可將循環卷積的周期縮短為L。

綜上所述，DTFT快速算法的基本思想如圖1所示，實現步驟如下：

步驟1確定g(n)和h(n)循環卷積的周期L，L≥K+N-1，且L為2的整數次冪。

步驟2將h(n)按式（7）轉化為一個長度為L點的新序列hL(n)：

步驟3根據已知參數按式（5）計算序列g(n)。

步驟4將g(n)補零成為長度為L的序列g1(n)，利用FFT算法計算g1(n)的L點離散傅里葉變換（Discrete Fourier Transform，DFT），記為G(l)。

步驟5利用FFT算法計算hL(n)的L點的DFT，記為H(l)。

步驟6將G(l)與H(l)相乘，得到g(n)和h(n)的循環卷積，記為Y(l)，并對其進行L點快速傅里葉逆變換（Inverse Fast Fourier Transform，IFFT）得到序列y′(n)。取y′(n)在0≤n≤K-1范圍內的值，即為g(n)和h(n)的線性卷積，記為y″(n)。

步驟7將y″(n)與表達式e-jπk2Δf/fs相乘，得到K點DTFT的值。

圖1 DTFT快速算法的基本思想

2.1.2 DTFT快速算法計算量分析

（1）計算g(n)：根據式（5），需要N次復數乘法。

（2）3次計算L點FFT（步驟4、5和6）：共需要1.5L lbL次復數乘法和3L lbL次復數加法。

（3）計算g(n)和h(n)的循環卷積Y(l)：需要L次復數乘法。

（4）計算y″(n)與e-jπk2Δf/fs的乘積：需要K次復數乘法。

故在DTFT快速算法中，計算K點DTFT的計算量約為S′DTFT×次復數乘法和S′DTFT+次復數加法：

直接K點DTFT需要復數乘法次數SDTFT×=KN，需要復數加法次數SDTFT+=(N-1)K。當K，N取值變大時，DTFT快速算法的計算量將遠小于常規DTFT算法。

2.2 加權融合頻譜矩陣降維處理

由文獻[4]可知，加權融合頻譜矩陣X′(fP)是由P×1階矩陣Xm(fP)經M×P×Q(M≥2，P≥1，Q≥1)階加權因子e-jD加權融合生成，X′(fP)是一個P×Q階的矩陣，可用圖2所示的圖形表示。圖2是由P根橫線和Q根豎線編制的一個網格圖，X′(fP)中第(p，q)個元素X′q[fP(p)]用圖中第p條橫線和第q條豎線的交點表示，最優加權融合頻譜Xq0'(fP)由第q0條豎線上所有交點表示。

圖2 加權融合頻譜矩陣降維原理示意圖

（1）當P≠Q時，由DTFT頻率估計原理[10]可知，Xq0'(fP)可能位于X′(fP)任一列上。只能通過峰值搜索abs[Xq'(fP)]，根據其峰值元素所在列q0來確定Xq0'(fP)。因此必須計算出X′(fP)中P×Q個元素，計算量較大。

（2）當P=Q時，①當被測頻率f恰好等于序列fP中第p1(p1∈[1，P])個元素，即f=fP(p1)。其中，fP表示將fscope線性等分構成的一個長度為P的序列。在無噪聲條件下，abs[Xq'(fP)]的峰值元素必定位于X′(fP)的 (p1，q1)處，此時p1=q1，即在X′(fP)的對角線元素上。②當f等于fP中第p1個元素和第p2(p2∈[1，P])個元素之間的某個值時，假設f更接近fP(p1)，即f≈fP(p1)，在無噪聲情況下，abs[Xq'(fP)]的峰值元素也必定位于矩陣X′(fP)的(p1，q1)處，也在X′(fP)的對角線元素上。綜上所述，只需設置P=Q，abs[Xq'(fP)]的峰值元素必定會位于X′(fP)的對角線元素上。為確定最優加權融合頻譜Xq0'(fP)只需計算由X′(fP)對角線上元素構成的P×1階矩陣Xpq'(fP)即可。通過峰值搜索abs[Xpq'(fP)]即可確定Xq0'(fP)在X′(fP)的列數q0。按此方法可將原融合算法中計算P×Q階X′(fP)降維為計算P×1階的Xpq'(fP)，從而大幅度減少計算量。

Xpq'(fP)={X1'[fP(1)]，X2'[fP(2)]，…，XQ'[fP(Q)]}（10）

同時，對X′(fP)進行降維處理還可提高算法的抗噪性。因為當P=Q時，雖然由理論分析Xq0'(fP)對應的q0應出現在X′(fP)的對角線元素上，但受噪聲干擾，若通過直接峰值搜索abs[Xq'(fP)]來確定的q0很可能不在X′(fP)的對角線上。但采用上述降維處理后，即使在噪聲干擾情況下，也能確保q0來源于X′(fP)對角線元素，因此從一定程度上提高了算法的抗噪性。

2.3 1/3主瓣相關性分析

在原融合算法中，通過計算同一信號的兩個不同類型頻譜——多段正弦信號的最優加權融合頻譜Xq0'(fP)和其累加頻譜X(fP)的頻域相關譜r(fP)，來抑制虛假譜峰和噪聲干擾，提高信號參數估計精度。由圖3可知，Xq0'(fP)僅在其1/3主瓣對應的頻段內與其X(fP)對應部分具有最大的相似度，因此沒有必要分析Xq0'(fP)和X(fP)在整個頻段內的相關性，僅需計算圖3中灰色區域的1/3主瓣相關譜r13(fP)即可，從一定程度上消減原融合算法的計算量，其數學表達式如式（11）和式（12）所示：

圖31 /3主瓣相關性分析示意圖

表1 多段正弦頻譜融合法計算量

表2 多段正弦快速頻譜融合法計算量

3 算法計算量分析

將算法所需的復數乘法和復數加法次數按如下關系轉化為實數乘法和實數加法：一次復數乘法需要四次實數乘法和兩次實數加法，一次復數加法需兩次實數加法。設置M=4，P=150，Q=150，Nm=50，N′=200，Lm=256(m∈[1，M])，Pl1l2=P/5，采用DSP2812[11]為核心處理器，設置其主頻為40 MHz（其峰值可達150 MHz）。以計算一次實數乘法運算需要四個機器周期，計算一次實數加法運算需要一個機器周期為例對算法計算量進行分析，分析結果如表3所示。原融合法的總耗時均約為22.868 ms，能夠滿足普通應用所需。本文提出的快速融合算法的總耗時約為3.575 ms，較原融合算法的運算速度提高了約84.37%。當提高DSP的主頻時，上述算法的耗時將進一步降低。

表3 算法計算量和耗時說明表

4 仿真實驗

（1）信噪比變化條件下的仿真實驗

為測試快速融合算法在不同信噪比條件下的性能，針對信號b1和信號b2，分別進行了11組實驗，每組包括1 000次Monte_Carlo實驗，各組實驗的信噪比設置如表5第一列所示，其余參數見表4。

表4 實驗參數設置

表5 不同信噪比條件下的頻率估計均方根誤差kHz

由表5所示的仿真結果可知，快速融合算法在估計信號b2頻率時的RMSE低于估計信號b1的頻率。在信噪比為-11～5 dB的范圍內，快速融合算法頻率估計的RMSE略高于原融合算法；在信噪比為-15～-13 dB的實驗范圍內，快速融合算法頻率估計的RMSE略低于原融合算法；這是因為快速融合算法中對加權融合頻譜矩陣的降維處理能夠進一步提高算法抗噪性，因此在信噪比較低時其性能略優于原融合算法。

（2）單段信號長度變化條件下的仿真實驗

為測試快速融合算法在單段信號長度變化條件下的性能，針對信號b1和信號b2，分別進行了10組實驗，每組包括1 000次Monte_Carlo實驗，每組b1和b2中第m段信號的長度Nb1m和Nb2m設置如表6所示，其余參數參見表4。

表6 多段正弦信號中單段信號長度設置表

由表7所示的仿真結果可知，隨著單段信號長度在表6所示范圍內變化?？焖偃诤纤惴ǖ念l率估計RMSE隨著單段信號長度的增加而遞減，略大于原融合算法頻率估計的RMSE，整體性能穩定。

表7 不同單段信號長度條件下的頻率估計均方根誤差kHz

（3）被測頻率變化條件下的仿真實驗

為測試快速融合算法在被測頻率f變化條件下的性能，進行了10組實驗，每組實驗中f設置如表8第一列所示，每組包括1 000次Monte_Carlo實驗，其余參數見表4。

表8 不同被測頻率f條件下的頻率估計均方根誤差kHz

由表8所示的仿真結果可知，隨著被測頻率f在7.5～12 MHz內變化，在估計信號b1頻率時，快速融合算法和原融合算法的頻率估計平均RMSE分別為38.27 kHz，37.41 kHz；在估計信號b2頻率時，快速融合算法和原融合算法的頻率估計平均RMSE分別為37.35 kHz，36.22 kHz。上述仿真結果表明快速融合算法的性能較為接近原融合算法，較穩定；在滿足奈特斯特采樣定理的前提下，能夠適用于不同的被測頻率對象，具有良好的普適性。

（4）采樣頻率變化條件下的仿真實驗

為測試快速融合算法在采樣頻率fs變化條件下的性能，進行了10組實驗，每組實驗中fs設置如表9第一列所示，每組包括1 000次Monte_Carlo實驗，其余參數設置見表4。

表9 不同采樣頻率fs條件下的頻率估計均方根誤差kHz

由表9所示的仿真結果可知，隨著采樣頻率fs在38～42.5 MHz內變化，在估計信號b1頻率時，快速融合算法和原融合算法的頻率估計平均RMSE分別為38.44 kHz，37.54 kHz；在估計信號b2頻率時，快速融合算法和原融合算法的頻率估計平均RMSE分別為36.99 kHz，36.04 kHz。上述仿真結果表明快速融合算法的性能較為接近原融合算法，在滿足奈特斯特采樣定理的前提下，其性能基本不依賴于采樣頻率。

5 總結

為進一步提高多段正弦信號頻譜融合算法的實時性，滿足雷達、聲納、電子對抗等應用需求，本文通過設計DTFT快速算法、降維處理加權融合頻譜矩陣和1/3主瓣相關性分析處理等措施來降低算法計算量，提出頻率估計的多段正弦信號快速頻譜融合算法。計算量對比和仿真實驗表明，多段正弦信號快速頻譜融合算法在精度損失極小的前提下，能夠大幅降低計算量；在信噪比極低的情況下(SNR≤-13 dB)，其性能略優于原融合算法。

[1]Xiao Wei，Tu Yaqing，Liu Lingbing.Parameters estimation of LFM signal based on fusion of signals with the same length and known frequency-difference[C]//IEEE Proceedings of 8th World Congress on Control and Automation（WCICA’2010），Jinan，China，2010：6776-6781.

[2]Wu Guibo，Zhang Lijun，Lu Hui，et al.A real-time andhigh-precision algorithm for frequency estimation byfusion multi-segment signals[J].Procedia Engineering，2011，15：2313-2317.

[3]劉良兵，涂亞慶.基于多段分頻等長信號融合的頻率估計方法[J].計算機工程與應用，2008，44（11）：170-175.

[4]肖瑋，涂亞慶，劉良兵，等.頻率估計的多段同頻正弦信號頻譜相關算法[J].電子與信息學報，2012，34（3）：564-570.

[5]Liu Liangbing，Tu Yaqing，Xu Baosong.A division ratiovariable delay method for VCO FM nonlinearity correction[C]//IEEE Proceedings of 7th World Congress on Control and Automation（WCICA’2008），Chongqing，China，2008：2692-2695.

[6]閻振華，黃建國，韓晶.改進的低信噪比短序列快速頻率估計算法[J].系統工程與電子技術，2009，31（8）：1990-1992.

[7]張海濤，涂亞慶.計及負頻率影響的科里奧利質量流量計信號處理方法[J].儀器儀表學報，2007，28（3）：539-544.

[8]Prioakis J G，Manolakis D G.Digital signal processing：principle，algorithms，and application[M].New Jersey：Prentice Hall，2006.

[9]胡廣書.數字信號處理——理論、算法與實現[M].北京：清華大學出版社，2003.

[10]Ayhan B，Trussell H J，Chow M Y，et al.On the use of a lower sampling rate for broken rotor bar detection with DTFT and AR-based spectrum methods[J].IEEE Transactions on Industrial Electronics，2008，55（3）：1421-1434.

[11]Li Fengpan，Lu Wenhua.Turbine governor system frequency measurement based on DSP[C]//2011 2nd International Conference on Control，Instrumentation and Automation（ICCIA），2011：459-462.

XIAO Wei1,TU Yaqing1,SHEN Yanlin1,SU Dan1,ZHANG Lei2

1.Department of Logistical Information Engineering,Logistical Engineering University,Chongqing 401311,China
2.Zhuozhou Comprehensive Storehouse,Baoding,Hebei 611730,China

The spectra fusion method for multi-sections sinusoids（called“the fore-fusion method”for short）is an effective way of estimation frequency for the sinusoids in low Signal-to-Noise Ratio（SNR）,which has an important theoretical significance and practical value.In order to meet the high real-time demand in some fields such as radar,sonar and electronic countermeasures,a fast spectra fusion algorithm for multi-section sinusoids is put forward.This proposed algorithm can reduce the calculation and improve the real-time characteristic by the following techniques：design a fast DTFT algorithm,reduce dimensions of the weighted fusion spectrum matrix,and analyse the correlation of the 1/3 main-lodes of the optimization weighted-accumulation spectrum and the accumulation spectrum.The principles of the above techniques are expatiated.Calculation analyses and simulations demonstrate that the proposed algorithm can reduce most calculation of the fore-fusion method with lower little precision,and it works better in very low SNR(SNR≤-13 dB).

fast spectra fusion;multi-sections sinusoids;fast Discrete Time Fourier Transform（DTFT）algorithm;dimensions reduction;correlation of 1/3 main-lodes

A

TN957

10.3778/j.issn.1002-8331.1207-0435

XIAO Wei,TU Yaqing,SHEN Yanlin,et al.Fast frequency estimation algorithm based on spectra fusion of multi-section sinusoids.Computer Engineering and Applications,2014,50（6）：186-191.

國家自然科學基金（No.61271449，No.61302175）；重慶市自然科學基金（No.cstc2012jjA0877，No.cstc2011BA2015）；后勤工程學院青年科研基金項目資助課題。

肖瑋（1982—），女，博士，研究領域為信號處理，嵌入式系統；涂亞慶（1963—），男，博士，教授，博士生導師，研究領域為信號處理；沈艷林（1987—），男，碩士研究生，研究領域為信號處理；蘇丹（1983—），女，博士研究生，研究領域為信號處理；張磊（1980—），男，博士，研究領域為軍事物流與信號處理。E-mail：WZWRY@163.com

2012-07-29

2012-11-23

1002-8331（2014）06-0186-06

CNKI網絡優先出版：2012-12-13，http://www.cnki.net/kcms/detail/11.2127.TP.20121213.1629.003.html