噪聲背景下頻率估計(jì)方差的濾波方法

噪聲背景下頻率估計(jì)方差的濾波方法

fft幅度比值的頻率插值方法近年來，采用fft估計(jì)正交信號(hào)頻率的方法在振動(dòng)信號(hào)分析領(lǐng)域引起了人們的注意。由于FFT得到的是離散頻譜,若信號(hào)觀測(cè)時(shí)間為T,則譜線之間隔為Δf=1/T。因此直接利用FFT頻譜最大值對(duì)應(yīng)的頻率估計(jì)正弦信號(hào)的頻率精度受觀測(cè)時(shí)間長(zhǎng)度的限制,其誤差范圍為±Δf/2。增加觀測(cè)時(shí)間不但因FFT點(diǎn)數(shù)太多帶來實(shí)現(xiàn)上的困難,而且在條件受限場(chǎng)合不可能隨意增加觀測(cè)時(shí)間。為了突破頻率估計(jì)精度受到FFT頻率分辨率的限制,可以采用頻率細(xì)化的方法,但需要增加較多的計(jì)算量。當(dāng)信號(hào)頻率不是FFT的頻率分辨率Δf的整數(shù)倍時(shí),由于FFT的“柵欄”效應(yīng)引起頻譜泄漏,此時(shí)信號(hào)的實(shí)際頻率位于FFT主瓣內(nèi)兩條最大譜線之間,可以借助第二譜線與最大譜線的幅度比值估計(jì)信號(hào)的實(shí)際頻率在兩條譜線之間的位置(即基于FFT幅度比值的頻率插值方法)。為了抑制旁瓣通常在FFT之前對(duì)采樣數(shù)據(jù)進(jìn)行加窗處理,加窗使主瓣變寬,主瓣之內(nèi)出現(xiàn)三條以上譜線。能量重心法就是利用了主瓣內(nèi)的多條譜線來提高頻率估計(jì)精度。不加窗情況下,當(dāng)信號(hào)實(shí)際頻率接近FFT最大譜線對(duì)應(yīng)的頻率時(shí),由于確定第二大譜線容易發(fā)生錯(cuò)誤從而造成頻率插值方向錯(cuò)誤,因此限制了利用FFT幅度比值進(jìn)行頻率插值的精度。利用FFT主瓣內(nèi)第二大值與最大值之比的實(shí)部代替幅值之比,可以消除插值方向的錯(cuò)誤造成的頻率估計(jì)誤差。FFT最大值處的相位與信號(hào)的實(shí)際頻率和FFT幅度最大譜線對(duì)應(yīng)的頻率之間的偏差有關(guān),但通常信號(hào)的初相未知,因此不能直接利用FFT最大譜線的相位估計(jì)頻率。采用分段FFT的方法可以消除初相的影響,從而實(shí)現(xiàn)利用FFT的相位提高頻率估計(jì)精度。文獻(xiàn)中對(duì)上述方法的原理進(jìn)行了詳細(xì)的討論和分析,但缺少對(duì)各種方法的頻率估計(jì)性能的全面比較,特別是在噪聲背景中每種方法的頻率估計(jì)方差的對(duì)比以及與理論方差下限的比較。本文通過理論分析和計(jì)算機(jī)MonteCarlo模擬實(shí)驗(yàn)對(duì)常用的基于FFT的頻率估計(jì)方法的方差進(jìn)行分析和比較,以確定各種方法在不同信噪比及不同數(shù)據(jù)長(zhǎng)度條件下的頻率估計(jì)精度。1福色散法的頻率估算精度分析1.1復(fù)疊濾波器信號(hào)方差下限的確定在0～T時(shí)間內(nèi)對(duì)加性高斯白噪聲(AWGN)背景中單一頻率正弦信號(hào)按Δt=T/N進(jìn)行采樣,得到采樣序列s(n)=acos(2πf0Τn/Ν+θ0)+z(n)?n=0,1,2,?,Ν-1(1)式中z(n)為零均值高斯白噪聲,方差為σ2,采樣序列的信噪比為SNR=a2/(2σ2)。根據(jù)參數(shù)估計(jì)理論,在給定數(shù)據(jù)長(zhǎng)度和信噪比前提下,AWGN背景中信號(hào)參數(shù)的任意無偏估計(jì)方差不會(huì)小于某一確定值,即方差下限。方差下限有不同的確定方法,常用的為CramerRao下限(CR下限)。復(fù)正弦信號(hào)參數(shù)估計(jì)的CramerRao方差下限在文中給出,在信號(hào)頻率不接近零或二分之一采樣頻率、初始相位未知情況下,實(shí)正弦信號(hào)頻率估計(jì)方差下限為相同條件下復(fù)正弦信號(hào)的方差下限的2倍σ2CRB=3Τ2Νπ2?SΝR(2)圖1所示為N分別為64,128,256,512,1024及2048,信噪比從-10dB到40dB變化時(shí),CramerRao標(biāo)準(zhǔn)差下限σCRB的變化情況,圖中給出的是σCRB與FFT的頻率分辨率Δf的比值。當(dāng)N=2048,信噪比為40dB時(shí),標(biāo)準(zhǔn)差下限約為10-4Δf;當(dāng)N=256,信噪比為10dB時(shí)標(biāo)準(zhǔn)差下限約為10-2Δf。1.2fft幅度估計(jì)算法的頻率估計(jì)s(n)的N點(diǎn)FFT記為S(k),鑒于實(shí)序列FFT的對(duì)稱性,只考慮離散頻譜的前N/2點(diǎn),在不加窗(即加矩形窗)情況下有S(k)=asin[π(k-f0Τ)]2sin[π(k-f0Τ)/Ν]?ejθ0Ν-1Ν(k-f0Τ)π]+Ζ(k)?k=0,1,2,?,Ν/2-1(3)S(k)幅度最大值處的離散頻率索引值記做m,m=int[f0T],int[x]表示取最接近x的整數(shù)。Z(k)為z(n)的FFT。利用m可對(duì)信號(hào)的頻率做粗估計(jì)?fc=mΔf。其中Δf=1/T為FFT的頻率分辨率,即相鄰譜線之間的間隔。在不加窗的情況下,當(dāng)信號(hào)頻率不是正好為Δf的整數(shù)倍時(shí),S(k)在主瓣之內(nèi)有兩條譜線。最大譜線的幅度為S1=|S(m)|,次大譜線的位置記為m2,m2=m±1,次大譜線的幅度為S2=|S(m2)|,則根據(jù)α=S2/S1,可得到信號(hào)的實(shí)際頻率與粗測(cè)頻率?fc之間的相對(duì)偏差δ=(f0-mΔf)/Δf的估計(jì)?δ=±α1+α(4)式(4)中的符號(hào)根據(jù)次大值是在最大值的左側(cè)還是右側(cè)來決定。利用不加窗頻譜次大譜線與最大譜線的幅度之比進(jìn)行頻率插值的優(yōu)點(diǎn)是插值公式簡(jiǎn)單。不考慮噪聲影響,FFT主瓣內(nèi)的次大譜線的幅度永遠(yuǎn)大于旁瓣幅度,因此插值不會(huì)出現(xiàn)方向錯(cuò)誤。但是在有噪聲的情況下,當(dāng)δ的絕對(duì)值較小時(shí),可能出現(xiàn)位于FFT頻譜最大值另一側(cè)第一旁瓣的幅度超過主瓣內(nèi)次大值的情況,從而造成頻率插值方向相反,引起較大的頻率估計(jì)誤差。圖2(a)所示為N=1024,信噪比為18dB時(shí),第二大和第三大譜線的幅度隨δ變化的情況,可見|δ|較小時(shí)很容易出現(xiàn)插值方向錯(cuò)誤。這種頻率估計(jì)方法的誤差由兩部分組成,|δ|較小時(shí)誤差主要由插值方向錯(cuò)誤所致,|δ|較大時(shí)頻率誤差主要由FFT幅度估計(jì)誤差所造成。Rife-Jane方法的頻率估計(jì)方差可表示為σ2f=(1-|δ|)2[(1-|δ|)2+δ2]Τ2Ν?SΝRsinc2(δ)+2δ2Τ2erfc[δ|sin(πδ)|π(1-δ2)√Ν?SΝR](5)式中erfc()為補(bǔ)誤差函數(shù)。|δ|較小時(shí),不加窗FFT幅度比值法的頻率估計(jì)方差遠(yuǎn)高于式(2)表示的CramerRao下限;|δ|接近0.5時(shí),上式第二項(xiàng)的影響可以忽略,頻率估計(jì)方差可近似為σ2f≈π2/(32T2N·SNR)≈1.015σ2CRB,非常接近CramerRao下限。為了抑制旁瓣,通常在FFT之前對(duì)采樣數(shù)據(jù)進(jìn)行加(非矩形)窗處理。加窗使主瓣變寬,主瓣內(nèi)出現(xiàn)多條譜線,同時(shí)也使得分別位于最大值兩側(cè)的第二大和第三大譜線更容易區(qū)分(圖2(b)所示為加Hanning窗后第二大和第三大譜線的幅度隨δ變化的情況),因此基本避免了頻率插值方向錯(cuò)誤。1.3qoinn方法Quinn提出了利用FFT主瓣內(nèi)次大譜線與最大譜線FFT系數(shù)復(fù)數(shù)值之比的實(shí)部進(jìn)行頻率插值的方法。設(shè)FFT的最大值處的離散頻率為m,則m-1和m+1分別位于最大值的兩側(cè),且其中一個(gè)為次大值。定義α1=R{S(m-1)/S(m)}及α2=R{S(m+1)/S(m)},式中R(x)表示x取的實(shí)部,分別計(jì)算?δ1=α1/(1-α1)和?δ2=-α1/(1-α2),則頻率插值可表示為?δ={?δ1,?δ1＞0??δ2＞0?δ2?其它(6)研究S(k)的相位,并注意到式中幅度項(xiàng)符號(hào)的變化,用分別表示S(k)在幅度最大值處、主瓣內(nèi)第二大值處及主瓣另一側(cè)第一旁瓣的相位,在不考慮噪聲情況下,有。因此,當(dāng)信號(hào)實(shí)際頻率大于mΔf,即次大值位于最大值右側(cè)時(shí),α1<0,α2>0,于是?δ1,?δ2都為正;反之,次大值位于最大值左側(cè)時(shí),α1>0,α2<0,于是?δ1和?δ2都為負(fù),因此式(6)可得到正確的頻率插值。與Rife-Jane方法相比,Quinn方法的突出優(yōu)點(diǎn)是當(dāng)|δ|很小時(shí),由于|S(m-1)|/|S(m)|與|S(m+1)|/|S(m)|非常接近,所以Rife-Jane方法容易受噪聲干擾出現(xiàn)插值方向錯(cuò)誤,而S(m-1)與S(m+1)的相位因?yàn)橄嗖?80°(見圖3),所以不容易受噪聲干擾相混,Quinn方法取S(m-1)/S(m)與S(m+1)/S(m)的實(shí)部,正是利用了相位信息來判斷插值方向,從而避免了Rife-Jane方法在|δ|較小時(shí)頻率估計(jì)誤差激增的問題。下面分析Quinn方法的頻率估計(jì)方差。利用α1和α2估計(jì)頻率的方差相同,下面僅分析其中之一并略去下標(biāo)。由于正弦信號(hào)的FFT在主瓣內(nèi)的兩條譜線相位相差π,S(k)在主瓣內(nèi)最大值和次大值處的復(fù)數(shù)值可分別表示為S(m)=A1ejue001φ+Z1,S(m2)=-A2ejue001φ+Z2,其中A1和A2分別為無噪聲時(shí)S1和S2的幅度,且有A1=|(S(m)|≈Na|sin(πδ)|/(2πδ)及A2=|(S(m2)|≈Na|sin(πδ)|/[2π(1-|δ|)]。則根據(jù)Quinn估計(jì)公式,有α=R{-A2ejφ+Ζ2A1ejφ+Ζ1}(7)將Z1和Z2分別表示為Z1=X1+jY1和Z2=X2+jY2,代入上式,得α=-R{A2cosφ-X2+j[A2sinφ-Y2]A1cosφ+X1+j[A1sinφ+Y1]}=-A1A2+A2U1-A1U2-X1X2-Y1Y2A21+|Ζ1|2+2A1U1(8)式中?？紤]到主瓣內(nèi)FFT系數(shù)的信噪比一般較高,|Z1|/A1及|Z2|/A2大于或接近1的概率很小,分析估計(jì)誤差時(shí)忽略這種小概率情況而假設(shè)|Z1|/A1?1及|Z2|/A2?1。上式可近似為α≈-A2A11+U1/A1-U2/A21+2U1/A1(9)根據(jù)前述條件,可將上式展開成泰勒級(jí)數(shù)并略去高次項(xiàng),α可近似表示為α≈-A2A1(1+U1A1-U2A2)(1-2U1A1)≈-A2A1(1-U2A2-U1A1)(10)代入Quinn的頻率估計(jì)公式,整理得?δ=α1-α=-A2(1-U2/A2-U1/A1)A1+A2-U2-U1A2/A1(11)同理將上式展開級(jí)數(shù)并略去高次項(xiàng),整理得?δ≈-A2A1+A2-A1U2+A2U1(A1+A2)2(12)所以相對(duì)頻率偏差估計(jì)的方差可表示為var(?δ)=A21var(U2)+A22var(U1)(A1+A2)4(13)根據(jù)定義及X2Y2FFT的性質(zhì)有σX2=σY2=σ2Z/2=Nσ2/2,代入上式得var(δ)=A12+A22(A1+A2)4Νσ22(14)代入A1,A2,得var(δ^)=(1-|δ|)2[(1-|δ|)2+δ2]Ν?SΝR?sinc2(δ)(15)因此Quinn方法的頻率估計(jì)方差為σf2=var(δ^)Τ2=(1-|δ|)2[(1-|δ|)2+δ2]Τ2Ν?SΝR?sinc2(δ)(16)對(duì)比式(5)和式(16)可見,Quinn方法的頻率估計(jì)方差與Rife-Jane方法的頻率估計(jì)方差的第一項(xiàng)相同。2相位差法頻率估算的精度分析2.1fft相位差估計(jì)頻率的估計(jì)由式(3),在FFT幅度最大值處的相位可近似為ue001φm≈θ0-(m-f0T)π=θ0+δπ,可見FFT幅度最大譜線的相位中含有相對(duì)頻率偏差δ的信息。但由于通常初相θ0未知,因此并不能直接利用估計(jì)δ。為了消除θ0的影響,可以將式(1)表示的采樣序列分為兩段長(zhǎng)度各為N/2的序列并分別進(jìn)行FFT。兩個(gè)FFT在最大值處的相位分別為,則利用可以得到的估計(jì)δ^2=Δφ2πΔf2=ΔφπΤ(17)式中Δf2=2/T為N/2點(diǎn)FFT的頻率分辨率。為了區(qū)別幅度插值法的相對(duì)頻偏,此處的相對(duì)頻率偏差記為δ2,δ2=(f0-nΔf2)/Δf2,n為N/2點(diǎn)FFT幅度最大值對(duì)應(yīng)的離散頻率。兩段FFT相位差法頻率估計(jì)方差為σf2=4Τ2Νπ2?SΝR?sinc2(δ2)(18)可見分段FFT相位差估計(jì)頻率的方差不僅與信噪比及FFT長(zhǎng)度有關(guān),還與δ2的值有關(guān)。當(dāng)δ2=0時(shí),σf2=(4/3)σCRB2,略高于CR下限。而當(dāng)δ2=±0.5時(shí),頻率估計(jì)方差為σf2=(π2/3)σCRB2≈3.3σCRB2,約為δ2=0時(shí)的2.5倍。為了進(jìn)一步降低頻率估計(jì)方差,可以利用第二大譜線的相位與最大譜線的相位的加權(quán)平均值來進(jìn)行頻率估計(jì)。兩段FFT幅度第二大值處的相位差記為Δue001φ′,則利用Δue001φ′估計(jì)δ2為δ^′2=Δφ′2πΔf2(19)由文獻(xiàn)可知var(δ^′2)/var(δ^2)=var(Δφ′)/var(Δφ)=(1-|δ2|)2/δ22。根據(jù)維納濾波理論,對(duì)δ^2和δ^′2進(jìn)行加權(quán)平均,得到δ2的估計(jì)值為δ^ˉ2=(1-|δ2|)2δ^2+δ22δ^′2(1-|δ2|)2+δ22(20)實(shí)際估計(jì)時(shí)可用δ^2來代替上式中的δ2。按式(20)估計(jì)頻率的方差為σf2=4Τ2Νπ2SΝR[δ22+(1-|δ2|)2]2?[(1-|δ2|)4sinc2(δ2)+δ24sinc2(1-|δ2|)](21)當(dāng)δ2=±0.5時(shí),經(jīng)濾波后σf2為濾波前的一半。利用FFT相位差估計(jì)頻率應(yīng)注意相位差的取值范圍,因?yàn)楸旧淼娜≈捣秶鶠?π～π,所以直接按計(jì)算相位差的結(jié)果范圍將在-2π～2π,為了得到正確的頻率估計(jì)值應(yīng)先對(duì)進(jìn)行折疊(即將變換到主值范圍)。但若直接在±π處設(shè)一個(gè)門限進(jìn)行折疊,在有噪聲時(shí)會(huì)出現(xiàn)錯(cuò)誤。無噪聲時(shí),當(dāng)δ2接近0.5時(shí),值接近(小于)π,但受噪聲影響可能超過π。簡(jiǎn)單的折疊方法將大于π的值都轉(zhuǎn)換成負(fù)值,從而造成頻率插值方向相反,帶來很大的頻率估計(jì)誤差。δ2接近-0.5時(shí)也會(huì)出現(xiàn)同樣問題(利用Δue001φ=arg{S2S*1}來計(jì)算S2與S1相位差,雖然直接將Δue001φ轉(zhuǎn)換到-π～π范圍內(nèi),但仍然存在上述問題)。為了盡量避免上述情況的出現(xiàn),當(dāng)接近±π時(shí)可利用次大譜線位置n2與最大譜線位置n之間的關(guān)系來判斷是否需要對(duì)進(jìn)行折疊。因?yàn)槿鬾2>n,則f0>nΔf2,應(yīng)有Δue001φ>0,若n2<n則f0<nΔf2,應(yīng)有。因此可按下面方法對(duì)進(jìn)行折疊:式中β>0,實(shí)驗(yàn)表明取β=π/3效果較好。N=512時(shí),在信噪比大于0的條件下,基本不出現(xiàn)頻率插值方向錯(cuò)誤。2.2fft的頻率和長(zhǎng)度文獻(xiàn)提出了利用兩個(gè)不同長(zhǎng)度的FFT在最大值處的相位差估計(jì)頻率的方法。設(shè)兩個(gè)FFT的長(zhǎng)度分別為N和M(M<N),與2.1節(jié)介紹的分段FFT相位差法不同的是,這里兩個(gè)序列是有重疊的,第二個(gè)序列是第一個(gè)序列的一部分。兩個(gè)FFT的幅度最大值分別出現(xiàn)在離散頻率m和n處,對(duì)應(yīng)的相位分別記為ue001φ2和ue001φ1,則頻率估計(jì)為f^=Ν(Ν-Μ)Τ[Ν-1Νm-Μ-1Μn-φ2-φ1π](22)由于兩個(gè)FFT的長(zhǎng)度不同,m和n不在同一個(gè)位置發(fā)生變化,很難消除ue001φ2-ue001φ1的模糊。為了避免相位差出現(xiàn)模糊,該方法只適用于|δ|較小的情況,而Rife-Jane方法在|δ|較小時(shí)誤差較大,因此這兩種方法可以互補(bǔ)應(yīng)用。3分段fft相位差法模擬結(jié)果為了對(duì)各種估計(jì)方法的精度進(jìn)行對(duì)比并驗(yàn)證前面給出的方差計(jì)算公式,進(jìn)行了計(jì)算機(jī)MonteCarlo模擬試驗(yàn)。為了方便與FFT本身的頻率分辨率Δf對(duì)比,下面在按公式計(jì)算和計(jì)算機(jī)模擬過程中均采用頻率估計(jì)的標(biāo)準(zhǔn)差σf來進(jìn)行分析,同時(shí)為了便于與估計(jì)誤差的理論下限進(jìn)行比較,下面給出的結(jié)果均為歸一化頻率估計(jì)標(biāo)準(zhǔn)差,即頻率估計(jì)標(biāo)準(zhǔn)差與相應(yīng)的標(biāo)準(zhǔn)差下限σCRB的比值。圖4中虛線所示為根據(jù)式(16),式(18)和式(21)計(jì)算的,用Quinn插值方法和分段FFT相位差法未進(jìn)行平均和加權(quán)平均之后的歸一化頻率估計(jì)標(biāo)準(zhǔn)差在一個(gè)FFT分辨單元內(nèi)變化情況的對(duì)比。因?yàn)榉侄蜦FT相位差法的一個(gè)FFT分辨單元相當(dāng)于Quinn方法的兩個(gè)FFT分辨單元(Δf2=2Δf),所以圖中δ2(分段FFT相位差法的相對(duì)頻偏)從0～0.5,而δ(Quinn方法的相對(duì)頻偏)從-0.5～0.5。可見Quinn方法在δ接近±0.5時(shí),標(biāo)準(zhǔn)差接近理論下限,但δ接近零時(shí),標(biāo)準(zhǔn)差達(dá)到下限的1.8倍。分段FFT相位差法經(jīng)兩點(diǎn)加權(quán)平均后,標(biāo)準(zhǔn)差隨δ2變化很小,約為1.2δCRB。Rife-Jane插值法的歸一化估計(jì)誤差與數(shù)據(jù)長(zhǎng)度及信噪比有關(guān)。圖4實(shí)線所示為加矩形窗時(shí)Quinn插值方法和分段FFT相位差法的歸一化頻率估計(jì)標(biāo)準(zhǔn)差模擬實(shí)驗(yàn)結(jié)果。模擬是在N=512,SNR=12dB,f0T=128+δ,初始相位隨機(jī)變化條件下進(jìn)行的,模擬次數(shù)為1000?？梢娔M結(jié)果與本文給出的公式計(jì)算結(jié)果很好地吻合。圖5所示為N=512,信噪比為0～30dB,Rife-Jane插值法加矩形窗時(shí)的歸一化估計(jì)誤差隨δ的情況(圖中所示為按式(5)計(jì)算結(jié)果。模擬結(jié)果與計(jì)算結(jié)果接近,圖中未示出)。可見當(dāng)信噪比較低時(shí),Rife-Jane方法的標(biāo)準(zhǔn)差在很大范圍內(nèi)都遠(yuǎn)高于σCRB,即使信噪比較高時(shí),該方法在|δ|<0.1時(shí)也很難得到滿意的精度,因此這種方法只能用在|δ|較大的情況。圖6為與圖4相同條件下Rife-Jane方法和重疊FFT相位差法的規(guī)一化標(biāo)準(zhǔn)差模擬結(jié)果與Quinn方法和分段FFT相位差法的對(duì)比。模擬結(jié)果表明:不加窗(即矩形窗)且|δ|較小時(shí),Rife-Jane方法誤差較大;|δ|較大(接近0.5)時(shí),Rife-Jane方法與Quinn方法誤差相同。重疊FFT相位差法在δ接近±0.5時(shí)出現(xiàn)較大誤差,|δ|較小時(shí)誤差略高于分段FFT相位差法。圖7所示為加Hanning窗時(shí)的規(guī)一化標(biāo)準(zhǔn)差的模擬結(jié)果,模擬條件與圖4和圖6相同。因?yàn)镕FT具有提高白噪聲背景中正弦信號(hào)信噪比的作用,其信噪比增益與數(shù)據(jù)長(zhǎng)度有關(guān),加(非矩形)窗使得有效數(shù)據(jù)長(zhǎng)度縮短,造成FFT頻譜的信噪比增益損失約3～5dB。因此與不加窗相比,加Hanning窗時(shí)的頻率估計(jì)誤差均略有增加(除Rife-Jane方法外,因?yàn)榧哟氨苊饬薘ife-Jane