摘 要:信號(hào)中含有噪聲或非整周期截?cái)鄷r(shí)發(fā)生的頻譜泄漏是導(dǎo)致正弦信號(hào)頻率估計(jì)精度不高的主要原因。針對(duì)這一問(wèn)題，從擴(kuò)展信號(hào)頻譜表征方式出發(fā)，將經(jīng)典幅值譜擴(kuò)展至不受頻譜泄漏制約、表現(xiàn)力更強(qiáng)、可讀性更好的二維幅值譜。與經(jīng)典幅值譜相比，二維幅值譜除包含信號(hào)的頻率個(gè)數(shù)、幅值信息外還包含了易于獲取的周期信息，且具有一定的抗噪性，在信噪比低至-10 dB時(shí)仍有較好表現(xiàn)力。提出一種估計(jì)方法，先從二維幅值譜中估計(jì)出信號(hào)的周期T，然后根據(jù)信號(hào)采樣頻率、信號(hào)頻率、以及信號(hào)周期T之間的定量關(guān)系完成信號(hào)的頻率估計(jì)。實(shí)驗(yàn)結(jié)果證明了該方法的有效性?；诙S幅值譜的正弦信號(hào)頻率估計(jì)方法為正弦信號(hào)的頻譜估計(jì)提供了一種新思路。

關(guān)鍵詞:DFT; 二維幅值譜; 頻率估計(jì); 正弦信號(hào)

中圖分類號(hào):TN911 文獻(xiàn)標(biāo)識(shí)碼:A

文章編號(hào):1004-373X(2010)13-0103-04

Sinusoid Signal Frequency Estimation Based on Two-dimensional Amplitude Spectrum

HAN Feng， TIAN Min， XU Gang

(Mechanical College， Inner Mongolia University of Technology， Huhhot 010051， China)

Abstract: The spectrum leakage producing as the signal has noise or non-integer period truncation is a principal reason of reducing the accuracy of the sinusoid signal frequency estimation. The typical amplitude spectrum of signals is expanded from one dimension to two dimension for solving that problem. The 2D amplitude spectrum shows more useful information of signals and is not limited by the signal integer period truncation. It contains the number of signal frequencies， amplitude information and easily obtained periodic information. Because of that， more useful information can be obtained from 2D spectrum than 1D′s. The 2D amplitude spectrum still has noise cancelling performance when SNR decreases to -10 dB. A new method to accomplish the frequency estimation of the signal according to the quantitative relation among sampling frequency， signal frequency and signal period T which is obtained by the estimation of 2D amplitude epectrum is proposed. It has an agreeable accuracy and good robustness.

Keywords: DFT; two-dimensional amplitude spectrum; frequency estimation sine signal

0 引 言

含噪聲正弦信號(hào)的頻率估計(jì)是信號(hào)處理領(lǐng)域研究的重要課題之一，被廣泛應(yīng)用于雷達(dá)、通信、聲納、電子對(duì)抗等多個(gè)領(lǐng)域。目前正弦信號(hào)頻率估計(jì)方法眾多[1-5]，各種現(xiàn)代譜估計(jì)方法可實(shí)現(xiàn)正弦信號(hào)頻率的精確估計(jì)，但這些方法算法復(fù)雜、計(jì)算量大[6]?；贒FT的經(jīng)典幅值譜(一維譜)分析方法，由于可利用FFT而具有運(yùn)算速度快、對(duì)正弦信號(hào)具有顯著的信噪比增益、算法參數(shù)不敏感等優(yōu)點(diǎn)而被廣泛使用[7]。但當(dāng)采樣頻率為DFT頻率分辨率[8]的非整數(shù)倍，即信號(hào)截?cái)嚅L(zhǎng)度不是信號(hào)周期整數(shù)倍時(shí)，信號(hào)頻譜就發(fā)生泄漏[9]。即使無(wú)噪聲影響，信號(hào)真實(shí)頻率仍落于主瓣內(nèi)兩根離散的譜線之間，導(dǎo)致頻率估計(jì)精度無(wú)法滿足要求。低信噪比情況下，頻率估計(jì)精度更低。基于FFT的各種插值算法改善了估計(jì)精度，但以增大計(jì)算量或復(fù)雜算法為代價(jià)，因此尋求一種算法簡(jiǎn)單、估計(jì)精度高，且具有一定抗噪聲能力的正弦信號(hào)頻率估計(jì)方法成為人們的研究目標(biāo)。本文從拓展信號(hào)幅值譜的表征方式出發(fā)，以Matlab為平臺(tái)，在DFT的基礎(chǔ)上將一維譜擴(kuò)展到突顯信號(hào)更多信息的二維幅值譜(簡(jiǎn)稱二維譜)，并在其基礎(chǔ)上推導(dǎo)出正弦信號(hào)頻率與信號(hào)采樣頻率、信號(hào)一個(gè)周期內(nèi)樣本點(diǎn)數(shù)以及頻率直線斜率之間的定量關(guān)系，進(jìn)而實(shí)現(xiàn)正弦信號(hào)的頻率估計(jì)。實(shí)驗(yàn)結(jié)果驗(yàn)證了該方法具有魯棒性好、估計(jì)精度高的特點(diǎn)。

1 二維譜的理論依據(jù)及解讀

1.1 二維譜的理論依據(jù)

設(shè)時(shí)域離散信號(hào)為x(n)，頻域?yàn)閄(k)，信號(hào)時(shí)域截?cái)嚅L(zhǎng)度為N，當(dāng)N為定值、k為變量時(shí)，可將X(k)看成關(guān)于k的函數(shù)，則根據(jù)有限長(zhǎng)序列的DFT公式:



X(k)=∑N-1n=0x(n)e-j2πkn/N，k=0，1，2，…，N-1 (1)

可做出x(n)在某個(gè)頻率點(diǎn)k的一維譜。若將DFT公式中信號(hào)時(shí)域截?cái)嚅L(zhǎng)度N、頻率點(diǎn)k均看作變量，則信號(hào)幅值與N，k均有關(guān)，可將信號(hào)幅值表達(dá)為X(k，N)，那么有限長(zhǎng)序列的DFT公式就可看作以信號(hào)時(shí)域截?cái)嚅L(zhǎng)度N和頻率點(diǎn)k為自變量、以信號(hào)幅值X(k，N)為因變量的二元函數(shù)，如式(2):

X(k，N)=∑N-1n=0x(n)e-j2πkn/N，k=0，1，2，…，N-1 (2)

式(2)確定了一個(gè)以信號(hào)長(zhǎng)度和頻率為自變量的二維信號(hào)譜。

對(duì)有限長(zhǎng)信號(hào)序列，非整周期截?cái)鄷r(shí)一維譜發(fā)生頻譜泄漏而使其有用信息可能淹沒于泄漏的能量之中，致使信號(hào)分析時(shí)無(wú)法解讀出足夠的有用信息;而對(duì)于二維譜(二維譜亦為雙邊譜，文中二維譜均指頻率點(diǎn)k取正值的右半邊譜)，由于將信號(hào)長(zhǎng)度亦作為變量使采樣信號(hào)中必包含信號(hào)的整周期點(diǎn)而使二維譜的可讀性受能量泄漏影響較一維譜小，故而二維譜突顯了信號(hào)更多有用信息，具有良好的可讀性。

1.2 二維譜的解讀

當(dāng)用來(lái)作二維譜的信號(hào)長(zhǎng)度滿足條件N≥8T時(shí)，二維譜頻率分辨率好[10]，展示出的信號(hào)周期(頻率)信息完整，具有較佳表現(xiàn)力，便于做譜分析來(lái)完成信號(hào)的頻率估計(jì)，故稱之為二維譜較佳表現(xiàn)力條件。與一維譜相比較，二維譜具有以下特點(diǎn):

(1) 直觀反映出信號(hào)周期及幅值信息

設(shè)T為信號(hào)一個(gè)周期里的樣本點(diǎn)數(shù)，N為信號(hào)的樣本點(diǎn)數(shù)，則T=100，N=1 000的正弦信號(hào)在k取1~10時(shí)的二維譜如圖1所示，可解讀為:一個(gè)直線走向的“山脊”表示信號(hào)中含有一個(gè)頻率成分，信號(hào)的頻率信息包含在直線狀的“山脊”中?！吧郊埂庇扇舾伞吧巾敗苯M成。實(shí)驗(yàn)觀察發(fā)現(xiàn):一個(gè)“山頂”對(duì)應(yīng)信號(hào)的一個(gè)周期，若作二維譜的信號(hào)長(zhǎng)度不是信號(hào)周期的整數(shù)倍時(shí)，則二維譜中會(huì)出現(xiàn)不完整的“山頂”。對(duì)“山頂”做定量分析有:二維譜中的“山頂”個(gè)數(shù)是作二維譜的信號(hào)中實(shí)際包含的信號(hào)周期個(gè)數(shù)的一半，若作二維譜的信號(hào)長(zhǎng)度不為T的整數(shù)倍，則“山頂”對(duì)應(yīng)的周期數(shù)為非整數(shù)。“山頂”高度對(duì)應(yīng)信號(hào)的幅值大小。

二維譜中表征信號(hào)周期及幅值信息的“山頂”對(duì)應(yīng)存在著尖點(diǎn)。對(duì)于確定的有限長(zhǎng)信號(hào)序列，在N方向，尖點(diǎn)隨著信號(hào)長(zhǎng)度N的變化遵從信號(hào)周期而周期性地出現(xiàn);沿k方向，尖點(diǎn)隨著信號(hào)長(zhǎng)度N的變化有規(guī)律地平移，尖點(diǎn)的這些特征包含了信號(hào)的周期信息，尖點(diǎn)高度表征了信號(hào)的幅值信息。暫稱這些尖點(diǎn)為信息點(diǎn)。有限長(zhǎng)信號(hào)序列的信息點(diǎn)位于一條直線上。根據(jù)上文對(duì)信息點(diǎn)的分析可知，信息點(diǎn)包含了精確的信號(hào)周期信息。這是基于二維譜實(shí)現(xiàn)正弦信號(hào)頻率估計(jì)的重要依據(jù)。

圖1 正弦信號(hào)的二維幅值譜

(2) 受制約條件少、頻率分辨率高

工程實(shí)際中，信號(hào)的周期信息往往是未知的，信號(hào)采集時(shí)很難做到整周期截取。對(duì)于有限長(zhǎng)信號(hào)，無(wú)干擾整周期截?cái)鄷r(shí)一維譜具有良好的頻率分辨率，此處指從譜中能夠判斷出信號(hào)所含的頻率個(gè)數(shù)[10]，非整周期截?cái)嗷蛐盘?hào)含有噪聲干擾時(shí)時(shí)頻譜泄漏導(dǎo)致一維譜頻率分辨率急劇下降。作二維幅值譜的信號(hào)長(zhǎng)度只要滿足二維譜較佳表現(xiàn)力條件，就可得到頻率分辨率高的二維譜，信號(hào)周期信息未知時(shí)，延長(zhǎng)信號(hào)采集時(shí)間即可滿足上述條件。對(duì)于信噪比低至-10 dB的含噪聲正弦信號(hào)，在滿足二維譜較佳表現(xiàn)力條件時(shí)，從其二維譜仍能清晰地判斷信號(hào)能量集中的位置，且易知該信號(hào)僅包含一個(gè)頻率成分，如圖2所示。

(3) 可解讀出更多信息

與一維譜相比較，二維譜不僅包含了頻率、幅值信息，還包含了幅值隨信號(hào)長(zhǎng)度N的變換規(guī)律以及幅值隨頻率點(diǎn)k的變化規(guī)律。從不同角度觀察二維譜可獲得N，k及幅值各物理量之間的相互變化規(guī)律。

① 沿信號(hào)長(zhǎng)度N方向縱切二維譜，即頻率點(diǎn)k一定時(shí)，隨著N的連續(xù)取值可得信號(hào)在某個(gè)頻率點(diǎn)對(duì)應(yīng)的幅值與信號(hào)長(zhǎng)度N之間的變化規(guī)律遵從|sinc|函數(shù)，如圖3(a)所示。頻率不同，對(duì)應(yīng)的|sinc|函數(shù)大小及出現(xiàn)周期亦不同。對(duì)于單頻信號(hào)，在采樣信號(hào)的長(zhǎng)度范圍內(nèi)N每變化一個(gè)周期里樣本點(diǎn)數(shù)即T對(duì)應(yīng)的長(zhǎng)度時(shí)，|sinc|函數(shù)沿N軸做T長(zhǎng)度的平移。

② 沿頻率點(diǎn)k方向縱切二維譜，即信號(hào)長(zhǎng)度N為定值時(shí)，即為經(jīng)典幅值譜，如圖3(b)所示。幅值與頻率點(diǎn)k的變化規(guī)律不再贅述。

2 基于二維譜的正弦信號(hào)頻率估計(jì)

根據(jù)二維譜理論依據(jù)，對(duì)信號(hào)序列做DFT，由變換后得到的數(shù)據(jù)做二維譜，觀察二維譜并根據(jù)二維譜規(guī)律特征建立模板將信息點(diǎn)所在的k，N位置找出。信息點(diǎn)的位置坐標(biāo)(k，N)二維數(shù)據(jù)在二維譜的k-N平面中表征信號(hào)的周期信息。在以頻率點(diǎn)k為橫軸、以信號(hào)長(zhǎng)度N為縱軸的二維平面(以下簡(jiǎn)稱k-N二維平面)中根據(jù)信息點(diǎn)的(k，N)坐標(biāo)描出信息點(diǎn)，可獲得一條由信息點(diǎn)組成的直線，稱之為信號(hào)的頻率直線，如圖4所示，其物理意義如圖5所示。實(shí)驗(yàn)發(fā)現(xiàn)，對(duì)有限長(zhǎng)信號(hào)做DFT時(shí)，隨著信號(hào)序列采樣位置的變化，信息點(diǎn)的位置有規(guī) 律地發(fā)生變化，但信息點(diǎn)組成的頻率直線的斜率不發(fā)生變化，即頻率直線所包含的信號(hào)周期信息不發(fā)生變化，這是采用分析頻率直線進(jìn)而實(shí)現(xiàn)正弦信號(hào)頻率估計(jì)的方法可靠性所在。下面依據(jù)此結(jié)論進(jìn)行正弦信號(hào)頻率估計(jì)原理的推導(dǎo)。

圖2 時(shí)域波形及其一維、二維幅值譜分析

圖3 幅值隨信號(hào)長(zhǎng)度的變化規(guī)律及不同截?cái)喾绞较路惦S頻率k點(diǎn)的變化規(guī)律

圖4 正弦信號(hào)信息點(diǎn)示意圖

圖5 頻率直線斜率的物理意義

設(shè)信號(hào)采集設(shè)備的采樣頻率為fs，信號(hào)采集時(shí)間為t，采集到的信號(hào)長(zhǎng)度為M，那么有:

M=fst (3)

若設(shè)信號(hào)實(shí)際頻率為f0，信號(hào)一個(gè)周期內(nèi)樣本點(diǎn)數(shù)為T，長(zhǎng)度為M的信號(hào)中實(shí)際包含的周期數(shù)(可為非整數(shù))為n，則有:

M=Tn (4)

t=n/f0 (5)

由式(3)~式(5)可推導(dǎo)出:

f0=fs/T (6)

式(6)即為正弦信號(hào)頻率估計(jì)計(jì)算公式。可知，T的估計(jì)是實(shí)現(xiàn)正弦信號(hào)頻率估計(jì)的關(guān)鍵。

在k-N二維平面中，N′表示頻率直線起止信息點(diǎn)之間的信號(hào)長(zhǎng)度，Δk表示起止信息點(diǎn)所占頻率點(diǎn)k的范圍差。則正弦信號(hào)一個(gè)周期里的樣本點(diǎn)數(shù)T可表示為:

T=N′/Δk (7)

設(shè)正弦信號(hào)的頻率直線為y=Ax+B，如圖5所示，易知其斜率A為N′/Δk，那么頻率直線斜率A的物理意義為:正弦信號(hào)一個(gè)周期里的樣本點(diǎn)數(shù)。則有:

A=N′/Δk=T (8)

故正弦信號(hào)頻率估計(jì)公式亦可表示為:

f0=fs/A (9)

以上分析可知，要實(shí)現(xiàn)正弦信號(hào)的頻率估計(jì)，只需求出頻率直線的斜率然后將其代入式(9)即可。

3 實(shí)驗(yàn)及結(jié)論分析

表1為在滿足二維譜較佳表現(xiàn)力條件下，對(duì)任意長(zhǎng)度(相對(duì)于一維譜要不發(fā)生頻譜泄漏而須滿足信號(hào)整周期截?cái)喽?、不同信噪比的正弦信號(hào)進(jìn)行頻率估計(jì)的實(shí)驗(yàn)結(jié)果。其中，N為做DFT的信號(hào)長(zhǎng)度;T1為正弦信號(hào)一個(gè)周期里樣本點(diǎn)數(shù)估計(jì)值;f1為正弦信號(hào)頻率估計(jì)值;εf為頻率估計(jì)的相對(duì)誤差。不含噪聲正弦信號(hào)的T=100、幅值為單位1、信號(hào)采樣頻率fs=100 Hz、信號(hào)頻率f0=1 Hz。頻率相對(duì)誤差為εf=f0-f1/f0。具體頻率估計(jì)精度及相對(duì)誤差如表1所示。由表1可知，對(duì)任意信號(hào)長(zhǎng)度不為周期整數(shù)

的正弦信號(hào)進(jìn)行基于二維譜的頻率估計(jì)時(shí)，信號(hào)長(zhǎng)度

N對(duì)頻率估計(jì)精度無(wú)影響。頻率估計(jì)精度隨信號(hào)信噪比的增大而升高，當(dāng)信噪比為-10 dB時(shí)，通過(guò)二維譜仍可準(zhǔn)確地估計(jì)出信號(hào)頻率，相對(duì)頻率誤差為0.026，而當(dāng)信噪比不斷增大時(shí)，相對(duì)頻率誤差不斷減小，信噪比SNR≥6 dB時(shí)，相對(duì)頻率誤差降至0.006 6，并趨于恒定。這個(gè)恒定誤差是由生成時(shí)域離散信號(hào)仿真軟件的系統(tǒng)誤差引起的，若不考慮上述因素，理論上本頻率估計(jì)方法的相對(duì)誤差可趨于零。

表1 不同信噪比下正弦信號(hào)頻率估計(jì)及精度(T=100，幅值為1，fs=100 Hz，f0=1 Hz)

SNR /dB-10-8-6-4-2024681012

N950940930920910890906908912932946984

T1102.666 7102.333 3101.333 3101.666 7101.666 7101.000 0101.000 0101.000 0100.666 7100.666 7100.666 7100.666 7

f1 /Hz0.974 00.977 20.986 80.986 80.983 60.990 10.990 10.990 10.993 40.993 40.993 40.993 4

εf0.026 00.022 80.013 20.013 20.016 40.009 90.009 90.009 90.006 60.006 60.006 60.006 6

4 結(jié) 語(yǔ)

通過(guò)擴(kuò)展傳統(tǒng)DFT計(jì)算公式中信號(hào)長(zhǎng)度N亦為變量，給出了能夠突顯出信號(hào)更多有用信息、具有一定抗噪聲能力、頻率分辨率高的二維譜。在信號(hào)長(zhǎng)度滿足N≥8T時(shí)二維譜具有較佳表現(xiàn)力和準(zhǔn)確度。二維譜的應(yīng)用可推廣至多頻信號(hào)的分析處理中，為信號(hào)的譜分析提供了一種新的頻譜表征方式。本文給出的是一種以適度增大數(shù)據(jù)量為代價(jià)獲得二維幅值譜，從而達(dá)到提高信號(hào)頻率估計(jì)精度和穩(wěn)定性的方法，為正弦信號(hào)的頻率估計(jì)方法提供了一種新的參考思路。

參考文獻(xiàn)

[1]祝俊，唐斌，陳真.快速高精度實(shí)正弦信號(hào)頻率估計(jì)算法[J].電子測(cè)量與儀器學(xué)報(bào)，2008，22(6):65-69.

[2]QUINN B G. Estimating frequency by interpolation using Fourier coefficients[J]. IEEE Transaction on Signal Processing， 1994， 42(5): 1264-1268.

[3]RIFE D C， VINCENT G A. Use of the discrete Fourier transform in the measurement of frequencies and levels of tones[J]. Bell. Syst .Tech.， 1970， 49(2): 197-228.

[4]FUNGA H W， ALEX C， KOTB K H， et al. Parameter estimation of a real signal tone from short data records[J]. Signal Processing， 2004， 84: 601-617.

[5]ZAKHAROV Y V， TOZER T C. Frequency estimator with dichotomous search of periodogram peak[J]. Electronics Letters， 1999， 35(19): 1608-1609.

[6]張松.基于FFT的正弦信號(hào)頻率估算新方法[J].大理學(xué)院學(xué)報(bào)，2009，8(8):36-39.

[7]黃玉春，黃載祿，黃本雄，等.基于FFT滑動(dòng)平均極大似然法的正弦信號(hào)頻率估計(jì)[J].電子與信息學(xué)報(bào)，2008，30(4):831-835.

[8]胡廣書.數(shù)字信號(hào)處理——理論、算法與實(shí)現(xiàn)[M].北京:清華大學(xué)出版社，1997.

[9]劉春艷，韓峰，梅秀莊.一種基于最小能量泄漏的信號(hào)周期估計(jì)方法[J].現(xiàn)代電子技術(shù)，2009，32(7):4-10.

[10]王剛，王艷芬，張曉光，等.關(guān)于離散傅里葉變換頻率分辨率的討論[J].電氣電子教學(xué)學(xué)報(bào)，2006，28(6):18-20.