馬羅文，劉寧，胡心語(yǔ)，倪麗花

（1.中國(guó)人民解放軍95865部隊(duì)，北京 102200；2.中國(guó)人民解放軍32802部隊(duì)，北京 100082；3.電子科技大學(xué)，四川 成都 611731）

0 引言

在雷達(dá)［1-3］、通信［4-5］、聲吶［6］等領(lǐng)域?qū)π盘?hào)頻率進(jìn)行快速、準(zhǔn)確的測(cè)量具有重要意義。目前，常用的頻率估計(jì)算法按照工作原理可以分為：計(jì)數(shù)法、相位法、快速傅里葉變換算法、子空間分解類算法等。其中，單頻信號(hào)［7］的估計(jì)問(wèn)題一般采用計(jì)數(shù)法和相位法；而多頻信號(hào)估計(jì)問(wèn)題一般采用的方法為子空間分解類算法和快速傅里葉變換算法。子空間分解類算法［8-9］中典型代表為多重信號(hào)分類（Multiple Signal Classification，MUSIC）算法［10］和旋轉(zhuǎn)不變技術(shù)估計(jì)信號(hào)參數(shù)（Estimating Signal Parameter Via Rotational Invariance Techniques，ESPRIT）算法［11］，此類子空間分解算法具有較高的估計(jì)精度，但由于需要搜索幅度譜峰值位置或特征值分解導(dǎo)致計(jì)算量相對(duì)較大。快速傅里葉變換算法［12］只能對(duì)有限點(diǎn)序列進(jìn)行處理，會(huì)造成能量泄露，同時(shí)存在柵欄效應(yīng)，當(dāng)信號(hào)頻譜峰值譜線不能與主瓣的中心重合時(shí)，信號(hào)真實(shí)頻率落于主瓣內(nèi)的兩根譜線之間，便會(huì)產(chǎn)生較大的頻率測(cè)量誤差［13-16］。

由于快速傅里葉變換的幅度譜頻率搜索間隔與信號(hào)采樣數(shù)或快速傅里葉變換的長(zhǎng)度成反比，為了降低快速傅里葉變換法幅度譜頻率搜索間隔的量化誤差對(duì)頻率測(cè)量精度的影響，常采用填零的方法增加信號(hào)采樣數(shù)或快速傅里葉變換的長(zhǎng)度。但是，快速傅里葉變換的計(jì)算量與信號(hào)采樣數(shù)或快速傅里葉變換長(zhǎng)度成正比，這種提高測(cè)頻精度的方法也會(huì)使頻率測(cè)量的計(jì)算量成倍地增加，導(dǎo)致其難以應(yīng)用于實(shí)時(shí)性要求較高的場(chǎng)景。

高精度和實(shí)時(shí)性要求的矛盾一直都是頻率估計(jì)在高速全數(shù)字解調(diào)器中應(yīng)用的瓶頸。因此，研究一種既保證高精度，又滿足低計(jì)算量的頻率估計(jì)算法，具有重要的應(yīng)用前景［17］。本文提出同時(shí)利用信號(hào)采樣的快速傅里葉變換的幅度譜峰值位置和頻率搜索間隔量化誤差校正值確定信號(hào)頻率的方法，能實(shí)現(xiàn)高精度快速測(cè)頻的目的。

1 信號(hào)模型

暫不考慮噪聲時(shí)，單頻復(fù)諧波信號(hào)表示為

式中：A為復(fù)正信號(hào)的幅度；f0為頻率，Hz；θ0為初相，rad。在不做特殊說(shuō)明的情況下，本文基于該信號(hào)展開(kāi)算法理論推導(dǎo)與實(shí)驗(yàn)分析。

經(jīng)過(guò)N點(diǎn)采樣后的離散信號(hào)表示為

對(duì)x(n)進(jìn)行N點(diǎn)快速傅里葉變換，得到

幅度最大值處的離散頻率索引值為m，根據(jù)X[k]幅度最大值處的位置m可以得到信號(hào)頻率的估計(jì)值為由于柵欄效應(yīng)的存在，當(dāng)信號(hào)頻譜峰值譜線不能與主瓣的中心重合時(shí)，信號(hào)真實(shí)頻率落于主瓣內(nèi)的兩根譜線之間，便會(huì)產(chǎn)生較大的頻率測(cè)量誤差。

因此，本文利用信號(hào)相位信息，消除了傳統(tǒng)方法僅使用幅度信息受到的搜索間隔的限制，提出了FFT幅相聯(lián)合的快速高精度頻率估計(jì)方法。

2 FFT幅相聯(lián)合的高精度頻率估計(jì)方法

首先對(duì)前N-1個(gè)信號(hào)采樣x(n)，n=0，1，2，…，N-2，采用快速傅里葉變換，得到頻譜序列y1(k)，為

然后經(jīng)過(guò)一位滑窗平移得到后N-1個(gè)信號(hào)采樣x(n)，n=1，2，…，N-1，并對(duì)其采用快速傅里葉變換，得到頻譜序列y2(k)，為

由于頻譜序列y1(k)和y2(k)是對(duì)具有相同頻率的采樣信號(hào)進(jìn)行FFT運(yùn)算得到的，那么y1(k)和y2(k)幅度譜的峰值位置相同。根據(jù)頻譜序列y1(k)和y2(k)，得到幅度譜z(k)為

根據(jù)式（6），信號(hào)頻率估計(jì)值所對(duì)應(yīng)的序列號(hào)，即幅度譜的譜峰位置m為

根據(jù)式（3），式（4）和式（5）在幅度譜的譜峰位置m處的比值為

FFT幅相聯(lián)合的高精度頻率估計(jì)f0為

式中：angle表示復(fù)數(shù)的相位角。

由式（9）可知，利用了頻譜序列y1(m)和y2(m)的相位差對(duì)頻率進(jìn)行估計(jì)，頻率的估計(jì)結(jié)果不再受限于頻率搜索間隔的量化誤差，從而實(shí)現(xiàn)高精度頻率估計(jì)。

3 FFT幅相聯(lián)合的快速高精度頻率估計(jì)方法

根據(jù)式（4）和式（5）可知，為了求解y1(k)與y2(k)，需要進(jìn)行兩次N-1點(diǎn)FFT變換，對(duì)此，提出FFT幅相聯(lián)合的快速高精度頻率估計(jì)方法，利用迭代，僅需要一次FFT變換求解y1(k)，然后通過(guò)對(duì)y1(k)進(jìn)行簡(jiǎn)單的加法和乘法迭代運(yùn)算就可以得到y(tǒng)2(k)，避免了FFT變換重復(fù)計(jì)算的問(wèn)題可實(shí)現(xiàn)快速高精度頻率估計(jì)。

由式（4）和式（5）得

式中：v為信號(hào)采樣的尾首樣本差。

根據(jù)頻譜序列y1(k)和信號(hào)采樣的尾首樣本差v，即可快速得到幅度譜z(k)為

由式（7）知信號(hào)頻率估計(jì)值所對(duì)應(yīng)的序列號(hào)m。

進(jìn)一步利用頻譜序列y1(k)和y1(k)+v在幅度譜峰位置m處的相位差，消除FFT算法受限于量化間隔的影響。根據(jù)式（8）和式（10），在幅度譜的譜峰位置m處，可以得到

頻率搜索間隔的量化誤差校正值u為

得到FFT幅相聯(lián)合的快速高精度頻率估計(jì)f0為

上述算法流程圖如圖1所示。

圖1 FFT幅相聯(lián)合快速高精度算法流程圖Fig.1 Flow chart of fast high precision algorithm based on FFT amplitude and phase combination

綜上，F(xiàn)FT幅相聯(lián)合的快速高精度頻率估計(jì)算法，一方面通過(guò)使用滑窗迭代算法，只需要進(jìn)行一次FFT變換就可以得到頻譜序列y1(k)和y2(k)，減少了冗余計(jì)算；另一方面，利用幅度譜峰y1(m)和y1(m)+v的相位差對(duì)頻率搜索間隔量化誤差進(jìn)行校正，達(dá)到了高精度快速測(cè)頻的目的。

4 計(jì)算復(fù)雜度分析

為了求解y1(k)與y2(k)，需要進(jìn)行2次N-1點(diǎn)FFT變換，計(jì)算復(fù)雜度為

由式（10）可知，通過(guò)一次FFT變換求解y1(k)，然后對(duì)y1(k)進(jìn)行簡(jiǎn)單的加法和乘法迭代運(yùn)算就可以得到y(tǒng)2(k)，從而避免了FFT變換重復(fù)計(jì)算的問(wèn)題，因此，與傳統(tǒng)FFT算法相比，本文所提的FFT幅相聯(lián)合的快速高精度頻率估計(jì)算法的計(jì)算復(fù)雜度節(jié)省了O((N-1)log(N-1)-2)，從而實(shí)現(xiàn)快速高精度頻率估計(jì)。

5 仿真結(jié)果與分析

為了驗(yàn)證FFT幅相聯(lián)合的快速高精度頻率估計(jì)算法的性能，進(jìn)行多次蒙特卡洛仿真實(shí)驗(yàn)，并與經(jīng)典的FFT算法進(jìn)行比較。測(cè)量信號(hào)頻率時(shí)，F(xiàn)FT算法只利用幅度譜的譜峰位置測(cè)量信號(hào)頻率，F(xiàn)FT幅相聯(lián)合的快速高精度頻率估計(jì)算法不僅利用了幅度譜的譜峰位置，還利用了頻率搜索間隔的量化誤差校正值。

仿真過(guò)程中，信號(hào)采樣數(shù)分四種情況，即N=32，64，128，256，信號(hào)采樣中包括歸一化頻率（相對(duì)于采樣頻率）等于0.2140，信噪比為18 dB的信號(hào)，以及歸一化頻率等于0.1070、干信比為11 dB的干擾。

當(dāng)N=32，64，128，256時(shí)，頻譜序列y1(k)和y2(k)的幅度譜及其相位如圖2至圖5所示。可見(jiàn)，頻譜序列y1(k)和y2(k)的幅度譜y1，y2的譜峰位置相同，驗(yàn)證了通過(guò)式（6）估計(jì)譜峰位置是合理的；且隨著采樣點(diǎn)數(shù)增加，這一個(gè)特征越明顯。

圖2 N=32時(shí)，y1(k)和y2(k)的幅度譜及相位Fig.2 Amplitude spectrums and phases of y1(k)and y2(k) for N=32

圖3 N=64時(shí)，y1(k)和y2(k)的幅度譜及相位Fig.3 Amplitude spectrums and phases of y1(k)r and y2(k) for N=64

圖5 N=256時(shí)，y1(k)和y2(k)的幅度譜及相位Fig.5 Amplitude spectrums and phases of y1(k)and y2(k) for N=256

當(dāng)N= 32，64，128，256時(shí)，F(xiàn)FT算法與FFT幅相聯(lián)合的快速高精度頻率估計(jì)算法的信號(hào)歸一化測(cè)頻均方根誤差的累計(jì)概率分布函數(shù)（CDF）如圖6所示?？梢?jiàn)，F(xiàn)FT幅相聯(lián)合的快速高精度頻率估計(jì)算法在信號(hào)采樣數(shù)N= 32時(shí)的CDF接近FFT算法在信號(hào)采樣數(shù)N= 256時(shí)的CDF，達(dá)到了高精度測(cè)頻的目的。

圖4 N=128時(shí)，y1(k)和y2(k)的幅度譜及相位Fig.4 Amplitude spectrums and phases of y1(k)and y2(k) for N=128

圖6 不同采樣數(shù)，F(xiàn)FT算法與FFT幅相聯(lián)合快速算法的歸一化測(cè)頻誤差的CDF及局部圖Fig.6 CDFs of normalized frequency measurement error of FFT algorithm and fast FFT amplitude and phase combination algorithm with different sampling numbers

由仿真結(jié)果可見(jiàn)，不同采樣數(shù)下，F(xiàn)FT算法的信號(hào)歸一化測(cè)頻均方根誤差分別為0.0091，0.0044，0.0023和0.0011；FFT幅相聯(lián)合的快速高精度頻率估計(jì)算法的信號(hào)歸一化測(cè)頻均方根誤差分別為0.0014，0.0007，0.0003和0.0002。可以得出，與FFT算法相比，利用同樣的采樣信號(hào)，當(dāng)N=32，64，128，256時(shí)，F(xiàn)FT幅相聯(lián)合的快速高精度頻率估計(jì)算法的歸一化均方根誤差分別降低了84.15%，84.57%，85.24%和84.73%，F(xiàn)FT幅相聯(lián)合的快速高精度頻率估計(jì)算法在信號(hào)采樣數(shù)N=32時(shí)的測(cè)頻精度接近FFT算法在信號(hào)采樣數(shù)N=256時(shí)的測(cè)頻精度，達(dá)到了高精度測(cè)頻的目的。

6 結(jié)論

本文提出了一種FFT幅相聯(lián)合的快速高精度頻率估計(jì)算法，一方面，只需要進(jìn)行一次FFT變換就可以得到兩個(gè)頻譜序列，從而減少了冗余計(jì)算；另一方面，利用兩個(gè)頻譜序列的相位差對(duì)頻率搜索間隔量化誤差進(jìn)行校正，達(dá)到了高精度快速測(cè)頻的目的。實(shí)驗(yàn)結(jié)果和分析表明，由于同時(shí)利用了幅度譜的譜峰位置與頻率搜索間隔的量化誤差校正值，該算法在信號(hào)采樣數(shù)N=32時(shí)的測(cè)頻精度接近FFT算法在信號(hào)采樣數(shù)N=256時(shí)的測(cè)頻精度，適用于需要快速、高精度測(cè)頻的應(yīng)用場(chǎng)合。未來(lái)的研究還可以利用更多的頻譜序列，進(jìn)一步提高頻率估計(jì)的精度。