溫兵會(huì) 毛衛(wèi)寧

(1 東南大學(xué)網(wǎng)絡(luò)空間安全學(xué)院 南京 210096)

(2 東南大學(xué)信息科學(xué)與工程學(xué)院 南京 210096)

0 引言

正弦波信號(hào)的頻率估計(jì)是數(shù)字信號(hào)處理領(lǐng)域的一個(gè)值得研究的經(jīng)典課題，在聲吶、雷達(dá)等領(lǐng)域應(yīng)用廣泛。采用離散傅里葉變換(Discrete Fourier transform,DFT)譜估計(jì)方法估計(jì)正弦波信號(hào)的頻率是一種常用方法，由于快速傅里葉變換(Fast Fourier transform,FFT)算法的高效性，該方法在工程上得到廣泛應(yīng)用[1?3]。利用FFT 譜估計(jì)法進(jìn)行正弦波信號(hào)頻率估計(jì)的方法很多，大體上可以分為兩類(lèi)。一類(lèi)是需要判別頻率修正方向的算法，如Rife算法[4?5]、M-Rife算法[6]和I-Rife算法[7]。Rife算法是正弦信號(hào)頻域頻率估計(jì)的經(jīng)典算法，通過(guò)頻譜插值對(duì)實(shí)際頻率相對(duì)于譜線最大值頻率的偏移量進(jìn)行估計(jì)，計(jì)算量小，但是存在兩個(gè)問(wèn)題：一是當(dāng)實(shí)際頻率在譜線最大值頻率附近時(shí)，頻率估計(jì)誤差相對(duì)較大；二是，其頻率估計(jì)精度易受噪聲的影響，低信噪比(Signal-to-noise ratio,SNR)時(shí)估計(jì)性能下降。M-Rife算法解決了Rife算法的第一個(gè)問(wèn)題，但增加了運(yùn)算量。I-Rife算法解決了Rife算法的兩個(gè)問(wèn)題，但進(jìn)一步增加了計(jì)算量。Rife算法、M-Rife算法和I-Rife算法存在的共同問(wèn)題是需要判別頻率修正方向，再計(jì)算頻率偏差。另一類(lèi)，基于FFT的頻率估計(jì)方法，不需要判別頻率修正方向，直接計(jì)算頻率偏差，如Candan算法[8]、Fang算法[9]和改進(jìn)的Fang算法[10?12]。Candan算法計(jì)算簡(jiǎn)單，但當(dāng)信噪比較低時(shí)，容易出現(xiàn)插值方向錯(cuò)誤，導(dǎo)致誤差較大。Fang算法通過(guò)對(duì)信號(hào)在時(shí)域補(bǔ)等信號(hào)長(zhǎng)度的零，采用FFT頻譜中最大譜線相鄰的兩根譜線幅度估計(jì)頻率偏差，提高了頻率估計(jì)性能，但增加了計(jì)算量。改進(jìn)的Fang算法在Fang算法基礎(chǔ)上提高了頻率估計(jì)性能，減少了計(jì)算量。Candan算法和Fang算法及其改進(jìn)算法存在的共同問(wèn)題是涉及到非線性函數(shù)的計(jì)算，增加了算法的復(fù)雜度。為此，本文提出了一種新的頻率估計(jì)算法，在N點(diǎn)FFT運(yùn)算基礎(chǔ)上，利用兩點(diǎn)細(xì)化的頻譜值估計(jì)頻率偏移量，分析了頻偏估計(jì)的偏差和算法的復(fù)雜度。相比于I-Rife和改進(jìn)的Fang算法，本文不需要判別頻率修正方向，降低了算法計(jì)算量和復(fù)雜度。

1 問(wèn)題描述

含有噪聲的復(fù)正弦信號(hào)為

其中，f0=(k0+δ)·?f為信號(hào)頻率，?f=fs/N為頻率分辨率，δ為數(shù)字頻率偏差，取值范圍在?0.5～0.5之間，fs為采樣頻率，k0為數(shù)字頻率；w(n)為高斯白噪聲，均值為0，方差為A和θ0分別為信號(hào)幅度和初相；N為信號(hào)長(zhǎng)度。

對(duì)式(1)所示的信號(hào)進(jìn)行N點(diǎn)DFT得到

式(2)中，W(k)表示高斯白噪聲的DFT，在不考慮噪聲的情況下信號(hào)傅里葉變換的幅值為

當(dāng)信號(hào)的真實(shí)頻率不等于頻率分辨率?f的整數(shù)倍時(shí)，利用FFT 估計(jì)頻率存在頻率偏差，數(shù)字頻率偏差δ取值范圍在?0.5～0.5之間，信號(hào)的真實(shí)頻率位于譜線最大值與次大值對(duì)應(yīng)的頻率之間。

假設(shè)信號(hào)經(jīng)過(guò)N點(diǎn)FFT 運(yùn)算后，頻譜幅度最大值為|Xk0|，對(duì)應(yīng)的數(shù)字頻率為k0，次大值為|Xk0+α|。Rife算法利用頻譜幅度的最大值和次大值估計(jì)頻率偏差，頻率估計(jì)公式為[4]

其中，頻率修正方向α=±1。當(dāng)|Xk0+1| >|Xk0?1|時(shí)，α=+1；當(dāng)|Xk0+1|<|Xk0?1|時(shí)，α=?1。

Rife算法計(jì)算量小，但在頻偏較小時(shí)，頻率估計(jì)誤差增大，此外，低信噪比時(shí)估計(jì)性能下降；I-Rife算法[7]是對(duì)Rife算法的改進(jìn)，利用頻譜細(xì)化技術(shù)計(jì)算譜值|X(k0±0.5)|作為頻率修正方向的判據(jù)，利用頻移技術(shù)將被估計(jì)頻率移至兩相鄰數(shù)字頻率的中心處，以確保Rife算法能夠準(zhǔn)確估計(jì)頻率。I-Rife算法在一定程度上克服了Rife算法的不足，但增加了計(jì)算量；改進(jìn)Fang算法(I-Fang)[10]利用DFT計(jì)算6個(gè)頻點(diǎn)的譜值來(lái)估計(jì)數(shù)字頻率偏差，避免了Fang算法中對(duì)信號(hào)補(bǔ)等長(zhǎng)度的零進(jìn)行2N點(diǎn)FFT運(yùn)算，減少了計(jì)算量，但利用I-Fang算法進(jìn)行頻率估計(jì)同樣涉及到非線性函數(shù)的計(jì)算，算法復(fù)雜度較高。I-Fang算法性能與I-Rife算法接近，但計(jì)算量進(jìn)一步增大，兩者在很低信噪比時(shí)，依然存在性能下降和不穩(wěn)定問(wèn)題。

2 本文算法和偏差分析

2.1 本文算法

為了提高低信噪比時(shí)的頻率估計(jì)性能，同時(shí)減小計(jì)算量，本文提出一種快速有效的頻率估計(jì)方法，運(yùn)用頻譜細(xì)化方法計(jì)算頻譜值|X(k0±0.5)|，利用這兩點(diǎn)譜值估計(jì)數(shù)字頻偏，不需要判別頻率修正方向，提高了頻率估計(jì)的穩(wěn)定性。

根據(jù)式(3)，k0點(diǎn)左右相鄰兩點(diǎn)頻譜幅度的比值為

式(5)中由于(π(k?δ)+π(k+δ))/2 =kπ，則|sin(π(k?δ))|=|sin(π(k+δ))|，假設(shè)π(k?δ)?N和π(k+δ)?N，則由公式(5)可得數(shù)字頻率偏差估計(jì)值：

當(dāng)k=1時(shí)，即利用FFT 頻譜中的|X(k0?1)|和|X(k0+1)|值來(lái)估計(jì)數(shù)字頻偏δ，對(duì)應(yīng)的表達(dá)式如下：

式(7)利用頻譜值|X(k0?1)|和|X(k0+1)|計(jì)算數(shù)字頻偏易受噪聲的影響，為了提高低信噪比下的頻率估計(jì)精度，利用頻譜細(xì)化獲得頻譜值|X(k0?0.5)|和|X(k0+0.5)|，代入式(6)得到

信號(hào)頻率的估計(jì)值為

本文提出的頻率估計(jì)算法的流程如下：

步驟1 對(duì)采樣信號(hào)進(jìn)行N點(diǎn)FFT 運(yùn)算，尋找幅度譜最大值對(duì)應(yīng)的數(shù)字頻率k0。

步驟2 利用頻譜細(xì)化方法計(jì)算|X(k0±0.5)|兩點(diǎn)譜值。

步驟3 利用式(8)和式(9)估計(jì)信號(hào)頻率。

2.2 偏差分析

本文算法得到的數(shù)字頻率偏差表示式簡(jiǎn)單，計(jì)算量小，進(jìn)一步分析表明，還可以得到無(wú)噪聲情況下頻偏的偏差閉合表示式，可用于頻率偏差校正，以進(jìn)一步提高測(cè)頻精度。

數(shù)字頻偏δ估計(jì)的偏差定義為

根據(jù)式(9)和式(10)，頻率估計(jì)的偏差為

因?yàn)閿?shù)字頻率偏差的偏差bδ與頻率偏差bf之間的關(guān)系是線性的，因此只需分析數(shù)字頻偏的偏差bδ。

sin(π(0.5+δ)/N)和sin(π(0.5?δ)/N)是關(guān)于δ的函數(shù)，對(duì)其進(jìn)行泰勒級(jí)數(shù)展開(kāi)，在無(wú)噪聲情況下，根據(jù)式(8)和式(10)可得

有噪聲時(shí)，定義與信噪比有關(guān)的數(shù)字頻偏的偏差為

可以證明，至少存在一個(gè)SNR值γ滿(mǎn)足bδ(γ)=0。證明過(guò)程如下：

根據(jù)式(8)和式(13)可得

其中，

根據(jù)公式(2)可得

其中，X+0.5=sin(π(0.5?δ))/sin(π(0.5?δ)/N)，W(k0+0.5)表示復(fù)高斯噪聲。

Wr(k0+0.5)～N(0,1)和Wi(k0+0.5)～N(0,1)，因此，可以得到包含高斯隨機(jī)變量的|X(k0+0.5)|的形式為

其中，K=cos(θ0?π(N?1)/N(0.5?δ))Wr(k0+0.5)+sin(θ0?π(N?1)/N(0.5?δ))Wi(k0+0.5)。

當(dāng)SNR滿(mǎn)足γ →0時(shí)，根據(jù)式(16)可得

可以證明W(k0+0.5)和W(k0?0.5)是不相關(guān)的。Wr(k0+0.5)、Wi(k0+0.5)、Wr(k0?0.5)和Wi(k0?0.5)不相關(guān)，且具有相同的分布，則根據(jù)式(17)和式(18)可得

另一方面，當(dāng)SNR取值滿(mǎn)足γ →+∞時(shí)，它將近似于無(wú)噪聲的情況，根據(jù)公式(12)，偏差bδ(+∞)和δ有相同的符號(hào)。

其中，sgn(·)表示符號(hào)函數(shù)。如上所述，高SNR和低SNR情況下，偏差的正負(fù)號(hào)是相反的。因此，存在至少一個(gè)SNR值γ0滿(mǎn)足bδ(γ0)=0。

3 計(jì)算量與性能分析

3.1 仿真實(shí)驗(yàn)及性能分析

信噪比定義為SNR =采樣頻率fs=48 kHz，信號(hào)長(zhǎng)度N=4800，頻率分辨率?f=fs/N=10?f=fs/N=10 Hz，信號(hào)頻率為f0=(256+0.10)?f/2。對(duì)于復(fù)正弦信號(hào)，在頻率、幅度、相位參數(shù)均未知的情況下頻率估計(jì)的克拉美羅界(Cramer-Rao lower bound,CRLB)[9]為

為了驗(yàn)證本文算法的頻率估計(jì)性能，把本文算法與Rife算法[5]、I-Rife算法[7]、Candan算法[8]和I-Fang算法[10]進(jìn)行頻率估計(jì)性能比較。1000次Monte Carlo仿真，信噪比變化范圍為?22～20 dB，數(shù)字頻偏δ=0.1，本文算法及以上幾種算法的頻率估計(jì)均方根誤差和頻率估計(jì)偏差隨信噪比的變化如圖1和圖2所示。

圖1和圖2表明，當(dāng)信噪比大于10 dB時(shí)，5種算法頻率估計(jì)性能接近；信噪比小于10 dB時(shí)，Rife算法的頻率估計(jì)性能最差，Candan算法性能次之，I-Rife算法、I-Fang算法和本文算法估計(jì)性能接近。

圖1 頻率估計(jì)均方根誤差隨信噪比變化Fig.1 Frequency estimation root mean square error varies with SNR

圖2 頻率估計(jì)偏差隨信噪比變化Fig.2 Frequency estimation bias varies with SNR

取信噪比SNR =?22 dB，數(shù)字頻率偏差變化范圍為?0.5～0.5，比較5種算法頻率估計(jì)的均方根誤差隨數(shù)字頻偏的變化，結(jié)果如圖3所示。本文算法頻率估計(jì)均方根誤差基本上不隨頻偏變化，較Rife算法有很大的改善，性能與I-Rife算法和I-Fang算法相當(dāng)。

圖3 頻率估計(jì)均方根誤差隨數(shù)字頻偏δ 變化(SNR=?22 dB)Fig.3 Frequency estimation root mean square error varies with digital frequency offset(SNR =?22 dB)

3.2 計(jì)算量分析

I-Rife算法除需要一次N點(diǎn)FFT 運(yùn)算外，還要利用DFT計(jì)算4個(gè)頻點(diǎn)的譜值。I-Rife算法共需要進(jìn)行(N/2·log2N+ 4N)次復(fù)數(shù)乘法和(N·log2N+4(N?1))次復(fù)數(shù)加法。Fang算法共需要進(jìn)行N·log2(2N)次復(fù)數(shù)乘法和2N·log2(2N)次復(fù)數(shù)加法。本文算法要計(jì)算一次N點(diǎn)FFT和2個(gè)頻點(diǎn)的譜值，共需要進(jìn)行(N/2·log2N+2N)次復(fù)數(shù)乘法和(N·log2N+2(N?1))次復(fù)數(shù)加法。改進(jìn)的Fang算法[10]要計(jì)算一次N點(diǎn)FFT和6個(gè)頻點(diǎn)的譜值，共需要進(jìn)行(N/2·log2N+6N)次復(fù)數(shù)乘法和(N·log2N+6(N?1))次復(fù)數(shù)加法。當(dāng)N取不同值時(shí)，4種算法計(jì)算量如圖4所示。

圖4 4種算法計(jì)算量隨信號(hào)長(zhǎng)度變化Fig.4 The calculation amount of the four algorithms varies with the signal length

圖4表明，F(xiàn)ang算法的計(jì)算量最大，I-Fang算法的計(jì)算量次之，I-Rife算法的計(jì)算量再次之，本文算法的計(jì)算量最小。隨著N取值的不斷增大，本文算法的計(jì)算量相比I-Rife算法和I-Fang算法更小，更有優(yōu)勢(shì)。

4 結(jié)論

本文提出了一種利用一次N點(diǎn)FFT和兩點(diǎn)細(xì)化的頻譜值|X(k0±0.5)|精確估計(jì)正弦信號(hào)頻率的方法，分析了頻偏估計(jì)的偏差和算法計(jì)算量。相比于I-Rife和I-Fang算法，本文算法不需要判別頻率修正方向，算法復(fù)雜度低，計(jì)算量小，在保證性能的同時(shí)，提高了算法的穩(wěn)定性和實(shí)用性。綜上，本文算法的整體性能優(yōu)于I-Rife算法和I-Fang算法，是一種快速有效的頻率估計(jì)方法。