張 威， 李 松, 劉進(jìn)忙, 聶 瑩

(1.空軍工程大學(xué) 防空反導(dǎo)學(xué)院，陜西 西安 710051；2.中國人民解放軍 93861部隊(duì)，陜西 三原 713800)

基于頻差修正的余弦信號(hào)頻率估計(jì)算法*

張 威1， 李 松1, 劉進(jìn)忙1, 聶 瑩2

(1.空軍工程大學(xué) 防空反導(dǎo)學(xué)院，陜西 西安 710051；2.中國人民解放軍 93861部隊(duì)，陜西 三原 713800)

針對(duì)余弦振動(dòng)信號(hào)的頻率高精度估計(jì)需求，提出了一種基于頻差修正的頻率估計(jì)算法。對(duì)連續(xù)時(shí)間信號(hào)進(jìn)行采樣后，使用Candan算法估計(jì)出頻差，運(yùn)用頻差對(duì)信號(hào)的頻率進(jìn)行修正。對(duì)修正后的信號(hào)使用Liang算法再次進(jìn)行頻偏估計(jì)。最后將2次估計(jì)得到的頻率值相加求得最終估計(jì)頻率。通過頻差修正，避免了Candan算法因插值方向錯(cuò)誤和Liang算法自身特點(diǎn)導(dǎo)致估計(jì)精度降低的問題，雖然增加了計(jì)算量，但并不影響信號(hào)實(shí)時(shí)處理。仿真結(jié)果表明：在相對(duì)頻偏為任意值的情況下，改進(jìn)算法的頻率估計(jì)均方誤差接近克拉美羅下限(CRLB)，性能優(yōu)于現(xiàn)有頻率估計(jì)算法。

余弦信號(hào)； 頻差修正； 頻率估計(jì)； 離散傅立葉變換

0 引 言

在眾多行業(yè)中，通常從實(shí)際測量的振動(dòng)信號(hào)中提取各種特征[1～3]，用以參數(shù)檢測、質(zhì)量評(píng)價(jià)、狀態(tài)監(jiān)視和故障診斷等。隨著各種數(shù)字信號(hào)處理方法的不斷出現(xiàn)和發(fā)展，振動(dòng)信號(hào)數(shù)字處理方法一直是近年來的研究熱點(diǎn)[4，5]。在真實(shí)環(huán)境下，因振蕩、旋轉(zhuǎn)、脈動(dòng)等而產(chǎn)生余弦振動(dòng)信號(hào)[6]的情況非常普遍，本文正是在此類情況下研究高斯噪聲背景下的余弦信號(hào)頻率估計(jì)算法。文獻(xiàn)[7]中提出的最大似然估計(jì)法估計(jì)誤差達(dá)到克拉美羅下限(Cramer-Rao lower bound,CRLB)，是最優(yōu)估計(jì)算法，但是算法復(fù)雜，運(yùn)算量大，無法實(shí)時(shí)處理。 基于離散傅立葉變換(discrete Fourier transform,DFT)的頻率估計(jì)算法[8]具有運(yùn)算速度快、算法參數(shù)不敏感等特點(diǎn)，得到了國內(nèi)外學(xué)者的廣泛關(guān)注[9～18]。

基于DFT的余弦信號(hào)頻率估計(jì)算法主要分2步：1)利用快速傅氏變換(fast Fourier transform，F(xiàn)FT)對(duì)信號(hào)進(jìn)行粗估計(jì)；2)運(yùn)用插值方法進(jìn)行精估計(jì)。目前,此類算法的研究主要在步驟(2)。文獻(xiàn)[13]提出的Jacobsen算法，利用FFT頻譜中最大的三根譜線對(duì)頻率進(jìn)行精估計(jì)，但估計(jì)精度不高。文獻(xiàn)[14]中Candan針對(duì)Jacobsen算法系數(shù)進(jìn)行了糾正，估計(jì)精度有所提高,但在低信噪比時(shí)可能出現(xiàn)插值方向錯(cuò)誤，不能滿足高精度的特定場合。文獻(xiàn)[16]中Candan對(duì)文獻(xiàn)[14]的算法進(jìn)行了改進(jìn)，消除了估計(jì)偏差，但在低信噪比時(shí)頻率估計(jì)均方誤差明顯高于CRLB。文獻(xiàn)[17]中Fang L Y對(duì)信號(hào)補(bǔ)零后進(jìn)行2N點(diǎn)的DFT，推導(dǎo)了一種有近似運(yùn)算的雙根譜線頻率估計(jì)算法，比前面提到的算法精度高。文獻(xiàn)[18]中Liang X H也對(duì)信號(hào)補(bǔ)零后進(jìn)行2N點(diǎn)的離散傅里葉變換DFT，推導(dǎo)了一種新的無近似運(yùn)算的三根譜線頻率估計(jì)算法；同時(shí)利用DFT運(yùn)算的周期性和補(bǔ)零算法的特殊取值，實(shí)現(xiàn)了1次N點(diǎn)FFT運(yùn)算和5次單點(diǎn)DFT運(yùn)算完成估計(jì)，大大降低了運(yùn)算復(fù)雜度，同時(shí)提高了算法的估計(jì)精度和抗噪性。對(duì)信號(hào)補(bǔ)零后再進(jìn)行頻率估計(jì)的算法中存在一個(gè)問題：在頻率偏差較小時(shí)，估計(jì)方差接近CRLB；但當(dāng)頻率偏差較大時(shí)，估計(jì)方差遠(yuǎn)離CRLB。文獻(xiàn)[19]針對(duì)Rife方法不能在整個(gè)頻偏范圍內(nèi)保持較高精度的問題，提出了一種頻移修正的改進(jìn)方案。

本文根據(jù)文獻(xiàn)[14]和文獻(xiàn)[18]算法，結(jié)合文獻(xiàn)[19]的思想提出了一種基于頻差修正的改進(jìn)算法：將采樣后的信號(hào)通過希爾伯特變換變成復(fù)信號(hào)，運(yùn)用Candan算法估計(jì)出頻率偏差，再用此偏差對(duì)原始信號(hào)進(jìn)行修正。然后使用Liang算法對(duì)修正后的信號(hào)進(jìn)行頻偏估計(jì)。將兩次估計(jì)得到的頻率值相加求得信號(hào)最終估計(jì)頻率。該算法適當(dāng)增加了計(jì)算復(fù)雜度，但估計(jì)精度更高，更加接近于CRLB。

1 信號(hào)模型

在沒有噪聲的情況下，任意單頻余弦信號(hào)表示為

x(t)=cos(2πfdt+φd)

(1)

式中fd和φd分別為余弦信號(hào)的頻率和初始相位。

經(jīng)過采樣頻率為fs、采樣點(diǎn)數(shù)為N的采樣后，通過希爾伯特變換，在高斯噪聲背景下的復(fù)信號(hào)表示為

(2)

式中A為復(fù)信號(hào)的振幅；w(n)為均值為零、方差為σ2的復(fù)高斯白噪聲。

2 頻率估計(jì)算法

一般情況下，N點(diǎn)采樣后的信號(hào)真實(shí)頻率為f=(km+δ)fs/N，其中,δ<|0.5|，km為正整數(shù)?；贒FT的頻率估計(jì)現(xiàn)有算法中，在進(jìn)行頻率粗估計(jì)時(shí)均使用FFT找到頻譜峰值位置，算法的差異主要體現(xiàn)在偏差估計(jì)階段，表1中列出了近幾年來部分算法的偏差求解表達(dá)式。

3 基于頻差修正的改進(jìn)算法

Candan算法[14]計(jì)算簡單、精度較高，但在推導(dǎo)過程中只考慮了大信噪比的情況，忽略了噪聲對(duì)信號(hào)的影響；當(dāng)信噪比較低時(shí)，容易出現(xiàn)插值方向錯(cuò)誤，導(dǎo)致誤差較大。Liang算法[18]對(duì)信號(hào)作2N點(diǎn)DFT運(yùn)算，推導(dǎo)過程無近似運(yùn)算，達(dá)到理論精度；但當(dāng)偏差δ較大時(shí)，算法精度受限。本文算法綜合以上兩種算法的優(yōu)點(diǎn)，主要思路如下：對(duì)信號(hào)進(jìn)行N點(diǎn)采樣，使用Candan算法估計(jì)出信號(hào)頻差，用該頻差對(duì)采樣信號(hào)進(jìn)行頻率修正。對(duì)修正后的信號(hào)用Liang算法進(jìn)行頻偏估計(jì)。將2次估計(jì)得到的頻率值相加求得信號(hào)最終估計(jì)頻率。操作過程如圖1所示，具體步驟如下：

表1 基于DFT的正弦信號(hào)頻率估計(jì)算法

對(duì)信號(hào)x(n)作N點(diǎn)DFT，則譜峰位置km及左、右兩點(diǎn)處的DFT值為

X(km-1)=Aejφdf(δ+1)+W(km-1)

(3)

X(km)=Aejφdf(δ)+W(km)

(4)

X(km+1)=Aejφdf(δ-1)+W(km+1)

(5)

式中W(km)為高斯白噪聲w(n)的DFT，函數(shù)f(·)定義為

對(duì)f(δ)進(jìn)行泰勒級(jí)數(shù)展開，保留δ一次項(xiàng)，忽略其高次項(xiàng)得

(7)

在高信噪比時(shí)，可以用信號(hào)的DFT值代替Af(·)，代入式(7)得

(8)

2)對(duì)原始信號(hào)x(n)進(jìn)行頻率修正，得

(9)

式中w1(n)為高斯白噪聲w(n)經(jīng)頻率修正后的表達(dá)式。

3)在頻偏估計(jì)階段，使用Liang算法對(duì)信號(hào)x1(n)補(bǔ)零至2N點(diǎn)，作2N點(diǎn)DFT運(yùn)算，搜索到譜峰位置k2，估計(jì)出頻偏2。忽略w1(n)的影響，其推導(dǎo)過程如下:

假設(shè)補(bǔ)零后的信號(hào)為

x′(n):{x1(0),x1(1),x1(2),…,x1(N-1),0,…,0}

(10)

對(duì)x′(n)進(jìn)行DFT得

(11)

式中 ξ=fdN/fs。

簡記YΔ=|Ω(k2+Δ)|，Δ=-1,0,+1，則譜峰位置k2及其左、右兩點(diǎn)的DFT值為

(12)

(13)

(14)

由于，|δ|<0.5，N>0，所以，式(12)～式(14)可以去掉模值符號(hào)。由式(12)和式(14)進(jìn)行簡化運(yùn)算后得

(15)

式(13)代入式(15)可得

(16)

則偏差估計(jì)表達(dá)式為

(17)

4)最終頻率估計(jì)表達(dá)式為

(18)

圖1 改進(jìn)算法框圖

4 性能和復(fù)雜度分析

仿真參數(shù)設(shè)置如下：A=1；φd=0；fs=1 (歸一化)；N=256；fd=(N/4+δ)fs/N,|δ|<0.5；所加噪聲為高斯白噪聲，其均值為0;方差為σ2。MonteCarlo仿真次數(shù)為20 000次模擬，計(jì)算的均方誤差(meansquareerror,MSE)歸一化處理。

CRLB[7]為

CRLB=3f2sσ2/(2π2N(N2-1)A2)

(19)

信噪比定義為

SNR=10lg{A2/(2σ2)}

(20)

在圖2(a)和(b)中，分別畫出了SNR=10dB和SNR=5dB下不同方法估計(jì)頻偏的MSE。從圖中可以看出：本文算法在0≤δ≤0.49整個(gè)范圍內(nèi)，性能比較穩(wěn)定；在0.418≤δ≤0.49時(shí)，本文算法的MSE稍高于其他算法，但通過局部放大圖可以看出，性能幾乎未下降。

在圖2(c)中，模擬真實(shí)情況下大多數(shù)信號(hào)的信噪比范圍，分別畫出了不同方法估計(jì)頻偏的MSE。從仿真圖可以看出：本文算法在整個(gè)信噪比范圍內(nèi)，優(yōu)于其他算法，更接近于CRLB；當(dāng)SNR=-7dB，本文算法的MSE略高于Candan和Liang算法，但通過局部放大可以發(fā)現(xiàn)MSE均在10-5量級(jí)上。

圖2 不同δ時(shí)的MSE對(duì)比

文獻(xiàn)[18]針對(duì)運(yùn)算量大的問題,根據(jù)補(bǔ)零后的信號(hào)本身特點(diǎn)和DFT的周期性，提出了采用1次N點(diǎn)FFT和5次單點(diǎn)DFT來優(yōu)化算法。所有算法都采用FFT實(shí)現(xiàn)，其復(fù)雜度對(duì)比如表2所示。

5 結(jié) 論

在Candan算法和Liang算法的基礎(chǔ)上，綜合2種算法的優(yōu)缺點(diǎn)，本文提出了一種改進(jìn)的頻率估計(jì)算法:算法使用復(fù)雜度低的Candan算法估計(jì)出頻差，對(duì)原始信號(hào)先進(jìn)行一步修正；針對(duì)修正后的信號(hào)，使用Liang算法估計(jì)出剩余頻偏。相對(duì)于其他算法，本文算法犧牲了運(yùn)算復(fù)雜度而提高了估計(jì)精度，但是并不影響信號(hào)實(shí)時(shí)處理以及硬件復(fù)雜度。仿真結(jié)果表明:在任意頻偏下，改進(jìn)算法的頻率估計(jì)MSE更接近于CRLB。針對(duì)精度要求低的應(yīng)用場合，為了減少計(jì)算量，采用Liang算法也能取得好的性能。但是對(duì)于那些高精度和高穩(wěn)定度的場合，采用本文算法更加合適。本文算法未考慮加窗情況下的頻率估計(jì)精度問題，下一步應(yīng)該關(guān)注此方向。

表2 算法的復(fù)雜度對(duì)比

[1] 胡 燦，王旭峰，孟祥營，等.塔里木油田鉆機(jī)平臺(tái)振動(dòng)測試與振動(dòng)激勵(lì)源分析[J].傳感器與微系統(tǒng)，2014,33(8):56-58.

[2] 孫小超，周 文，李明濤，等.基于PSD的轉(zhuǎn)軸振動(dòng)位移測量與分析[J].傳感器與微系統(tǒng)，2014,33(9):135-137.

[3] 王威濤，姜印平.基于霍爾效應(yīng)的機(jī)械微振動(dòng)量檢測系統(tǒng)研究[J].傳感器與微系統(tǒng)，2013,32(10):34-36.

[4] 丁 康，陳健林，蘇向榮.平穩(wěn)和非平穩(wěn)振動(dòng)信號(hào)的若干處理方法及發(fā)展[J].振動(dòng)工程學(xué)報(bào)，2003,16(1):1-10.

[5] 肖棟林，吳 賓，向 函.基于雙通道異步采樣的振動(dòng)傳感器測量方法[J].傳感器與微系統(tǒng)，2016,35(1):146-148.

[6] 張步云，陳懷海，賀旭東.多輸入多輸出正弦掃頻試驗(yàn)控制新方法[J].振動(dòng)與沖擊，2015,34(8):198-202.

[7] Rife D C,Boorstyn R R.Single-tone parameter estimation from discrete-time observations[J].IEEE Transactions on Information Theory,1974,20(5):591-598.

[8] Rife D C,Vincent G A.Use of the discrete Fourier transform in the measurement of frequencies and levels of tones[J].Bell Sy-sten Technical Journal,1970,49(2):197-228.

[9] Quinn B G.Estimating frequency by interpolation using Fourier coefficients[J].IEEE Transactions on Signal Process,1994,42(5):1264-1268.

[10] Quinn B G.Estimation of frequency,amplitude, and phase from the DFT of a time series[J].IEEE Transactions on Signal Process,1997,45(3):814-817.

[11] Macleod M D.Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones[J].IEEE Transactions on Signal Processing,1998,46(1):141-148．

[12] 張英龍，劉 渝，王旭東.基于頻偏校正的正弦波頻率估計(jì)算法[J].南京航空航天大學(xué)學(xué)報(bào)，2007,39(5):597-600.

[13] Jacobsen E,Kootsookos P.Fast,accurate frequency estimator-s[J].IEEE Signal Processing Magazine,2007,24(3):123-125.

[14] Candan C.A method for fine resolution frequency estimation from three DFT samples[J].IEEE Signal Processing Letters,2011,18(6):351-354．

[15] Yang C,Wei G.A noniterative frequency estimator with rational combination of three spectrum lines[J].IEEE Transactions on Signal Processing,2011,59(10):5065-5070.

[16] Candan C.Analysis and further improvement of fine resolution frequency estimation method from three DFT samples[J].IEEE Signal Processing Letters,2013,20(9):913-916．

[17] Fang L Y,Duan D L,Yang L Q.A new DFT-based frequency estimator for single-tone complex sinusoidal signals[C]∥Procee-dings of the 2012 IEEE Military Communications Conference,MILCOM 2012:Piscataway:IEEE,2012:1-6.

[18] Liang X H,Liu A J,Pan X F,et al.A new and accurate estimator with analytical expression for frequency estimation[J].IEEE Communications Letters,2016,20 (1):105-108.

[19] 鄧振淼，劉 渝，王志忠.正弦波頻率估計(jì)的修正Rife算法[J].數(shù)據(jù)采集與處理，2006,21(4):473-477.

Frequency estimation algorithm for cosine signal using frequency-offset correction*

ZHANG Wei1, LI Song1, LIU Jin-mang1, NIE Ying2

(1.Air and Missile Defense College,Air Force Engineering University,Xi’an 710051,China; 2.Unit 93861 of the PLA,Sanyuan 713800,China)

To meet the demand of high precision frequency estimation of cosine vibration signal,a new algorithm based on frequency-offset correction is proposed.The algorithm estimates the frequency-offset,with which the frequency of signal can be revised,by Candan algorithm after sampling the continuous time signal.For the signal revised, use Liang algorithm to estimate the frequency-offset again.The ultimate estimated frequency is given by adding the two estimated frequency together.By frequency-offset revision,the proposed algorithm avoids the decrease of precision caused by interpolation error of Candan algorithm and characteristics of Liang algorithm.The new algorithm guarantees the real-time signal processing,although it increases the amount of computation.Numerical results demonstrate that the mean square error of frequency estimation of the improved algorithm is more closer to the Cramer-Rao lower bound(CRLB)bound character is prior to the existing frequency estimation algorithm.

cosine signal; frequency-offset correction; frequency estimation; discrete Fourier transform(DFT)

10.13873/J.1000—9787(2017)09—0121—04

2016—08—03

國家自然科學(xué)基金資助項(xiàng)目(61201287)

TN 911.7

A

1000—9787(2017)09—0121—04

張 威(1990-)，男，碩士研究生，主要研究方向?yàn)閿?shù)字信號(hào)處理。

李 松，通訊作者，E-mail：z199007w@163.com。