频率估计的多段正弦信号快速频谱融合算法

肖玮，涂亚庆，沈艳林，苏丹，张磊

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所示范围内变化。快速融合算法的频率估计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