對(duì)兩個(gè)頻率相近成分作頻譜校正的非迭代形式研究

張 強(qiáng)，張 頻，陳奎孚

(1.上海師范大學(xué) 建筑工程學(xué)院，上海 201418；2.江西農(nóng)業(yè)大學(xué) 國(guó)土資源與環(huán)境學(xué)院，南昌 330045；3.中國(guó)農(nóng)業(yè)大學(xué) 理學(xué)院，北京 100083)

對(duì)于非密集頻率復(fù)雜諧波信號(hào)，基于復(fù)正弦模型(Complex Sinusoid Model：CSM)的校正技術(shù)已經(jīng)發(fā)展得幾近完美［1－6］，有關(guān)的綜述參考文獻(xiàn)［7－8］。該方法僅僅利用主瓣附近幾條的離散譜線數(shù)據(jù)，經(jīng)過(guò)簡(jiǎn)單運(yùn)算就可以得到相當(dāng)精確的參數(shù)。對(duì)于頻率密集模型，是否還存在這種不需迭代的校正技術(shù)，尚未見(jiàn)文獻(xiàn)給出明確答案。最簡(jiǎn)單的頻率密集模型是雙頻率模型(Double-Frequency Model：DFM)。針對(duì)這種模型，文獻(xiàn)［9］通過(guò)搜索兩個(gè)頻率分量在復(fù)平面上的方位角將其分離，它需要使用數(shù)值法搜索。文獻(xiàn)［9－10］直接用數(shù)值方法迭代求解非線性頻率方程，因而兩者算法復(fù)雜性和運(yùn)算量都遠(yuǎn)遠(yuǎn)超出CSM的簡(jiǎn)單校正方法。

本文將探索對(duì)于DFM，是否存在類(lèi)似CSM的顯式校正公式。本文的理論分析表明：DFM加矩形窗情形確實(shí)有顯式校正公式；但對(duì)其他窗函數(shù)，校正公式要么不存在，要么過(guò)度復(fù)雜，因而不具有實(shí)用性。

1 CSM存在簡(jiǎn)單校正公式的原因

對(duì)周期信號(hào)和平穩(wěn)信號(hào)，所使用的窗函數(shù)一般均為對(duì)稱(chēng)。仔細(xì)研究CSM校正公式的推導(dǎo)過(guò)程，可以發(fā)現(xiàn)：它之所以簡(jiǎn)單是因?yàn)槌Ｒ?jiàn)對(duì)稱(chēng)窗函數(shù)w(t)的頻譜W(ω)具有如下的兩個(gè)特性。

特性1常見(jiàn)的對(duì)稱(chēng)窗函數(shù)的頻譜W(ω)一般可以表示為如下形式：

其中：WR(ω)為實(shí)有理分式，有的甚至為實(shí)多項(xiàng)式；WP(ω)是由三角函數(shù)構(gòu)成的超越實(shí)函數(shù)。

不失一般性，本文僅考慮WP(ω+Δω)=－WP(ω)的情形。

對(duì)于單頻復(fù)解析信號(hào) xa(t)=A0exp(jφ0+jω0t)加w(t)的傅里葉變換為：

記主瓣內(nèi)相鄰的兩條FFT譜線ω1和ω2對(duì)應(yīng)的頻譜為X1和 X2，即：

由于FFT的譜線間隔Δω=ω2－ω1=2π/T，因此二者的比值為：

其中利用了特性1和特性2。正因?yàn)檫@兩個(gè)特性，X1和X2的超越部分抵消而使得X2/X1僅為ω0的有理函數(shù)，進(jìn)而才有可能建立簡(jiǎn)單的校正公式。

比如將漢寧窗 WR(ω)=(ωT/2)［1－(πT/2)2］代入式(2)，立即可以得到它的校正公式。

對(duì)于DFM，不再有類(lèi)似于式(3)的簡(jiǎn)單比例關(guān)系。但是如果能將兩個(gè)所對(duì)應(yīng)頻譜分離，則各自參數(shù)應(yīng)滿足這個(gè)比例關(guān)系。

2 DFM的頻率方程

記DFM的信號(hào)為：

其中的 Aa，φa，ωa，Ab，φb，ωb分別為兩個(gè)頻率成分的幅值、相位和頻率參數(shù)。根據(jù)傅里葉變換的線性性質(zhì)，有：

考察雙頻率主峰附近的四條連續(xù)譜線。記頻率為ω1～ω4，對(duì)應(yīng)的頻譜為 X1～X4，即：

方程組(4)有8個(gè)待定參數(shù)，但是每個(gè)方程都是復(fù)方程，虛部和實(shí)部各提供一個(gè)方程，因此有8個(gè)方程。但是本文不直接解這8個(gè)方程，而是利用特性2。根據(jù)該特性，方程組(4)的每個(gè)方程右端兩個(gè)分子為：

盡管方程組(4)無(wú)超越函數(shù)，但也未必能夠得到顯式解。比如高于5次多項(xiàng)式的方程，有所謂的阿貝耳定理，即方程的根不可能用方程系數(shù)經(jīng)過(guò)有限次四則運(yùn)算和開(kāi)方運(yùn)算表達(dá)出來(lái)。退一步，即使有顯式解，若其表達(dá)式非常復(fù)雜，則實(shí)用性也不大。如簡(jiǎn)單的三次多項(xiàng)式方程，它的顯式根很復(fù)雜。

3 DFM加矩形窗的校正

下面研究DFM加矩形窗簡(jiǎn)化情形。

這樣方程組(4)的第2和第3兩式變?yōu)椋?/p>

圖1 頻率參數(shù)之間的關(guān)系Fig.1 The relationship between the frequency variables

不失一般性，取T=2。將WR(ω)=ω代入式(5)，利用克萊姆法則可以解出：

對(duì)于方程組(4)的第1和4式，只需將X1→X2，X3→X4，3δω→δω 代入式(6)即有：

將方程(8)的第1式展開(kāi)得：

將方程(8)第2式作類(lèi)似展開(kāi)，然后與式(9)相減，并消去公因式δa－δb得到：

式中：

將式(10)的第一項(xiàng)δaδb表示出來(lái)：

然后代入式(9)的第1項(xiàng)，經(jīng)整理可得另外一個(gè)方程：

其中：

聯(lián)立式(10)和式(11)，可以解出：

其中Δ，Δ1和Δ2是中間變量，具體為：

由式(12)可以看出頻率校正量仍然僅與頻譜比值有關(guān)。但是對(duì)于DFM，離散譜線在復(fù)數(shù)平面內(nèi)不再共線［11］，因此虛部和實(shí)部都出現(xiàn)在校正公式上。

由于Δ，Δ1和 Δ2均為復(fù)數(shù)，因此式(12)中－4Δ1Δ是否大于零并不重要。在理論上δa，δb應(yīng)該全為實(shí)數(shù)，但是由于誤差存在，很難保證其虛部為零，通常取其實(shí)部作為δa，δb的近似即可。

由圖1的關(guān)系，可以得到校正后的頻率：

以上是針對(duì)時(shí)間零點(diǎn)位于窗函數(shù)中心建立的公式。對(duì)于常用軟件包的FFT，零點(diǎn)位于窗口最左側(cè)。根據(jù)傅里葉變換的時(shí)移性質(zhì)，F(xiàn)FT的奇數(shù)條譜線將附加一個(gè)π相位，而偶數(shù)條譜線不變。因而所有推導(dǎo)過(guò)程和最終結(jié)果應(yīng)將 －X1→X1，－X3→X3。

由式(12)的推導(dǎo)過(guò)程可以看出，如果將矩形窗換成更復(fù)雜函數(shù)，則無(wú)法得到平行于式(10)和(11)的二次多項(xiàng)式方程，因此校正量雖然與頻譜的比值存在確定關(guān)系，但是顯式表達(dá)式將很復(fù)雜或者根本就得不到。

4 考核結(jié)果與討論

為了驗(yàn)證式(12)和式(13)的正確性，以及評(píng)價(jià)它們的實(shí)用性，計(jì)算參數(shù)取T=2，F(xiàn)FT的長(zhǎng)度N=1 024。因此，ΔT=T/N=1/512，Δω =2π/T=π。DFM 信號(hào)的兩個(gè)成分為：成分1為強(qiáng)幅成分，參數(shù)有兩組，即Aa=10，ωa=100.05Δω，φa= π/3 和 Aa=10，ωa=100.55 Δω，φa=π/3，前者接近整周期采樣，后者接近半周期采樣；成分2為弱幅成分，Ab=1，φb=5π/4。為了細(xì)致考核校正公式的效果，ωb在一個(gè)范圍內(nèi)取值，即ωb=ωa+ρΔω，其中 ρ=0.3～3Δω。

仿真考核的誤差如圖 2 所示，其中 εω/Δω，εA，εφ分別為校正后的頻率絕對(duì)誤差、幅值相對(duì)誤差和相位誤差。

圖2 仿真考核的誤差特性Fig.2 The error characteristics of the tested example

由圖2可以看出，強(qiáng)幅成分的參數(shù)精度很高，誤差的總趨勢(shì)隨兩個(gè)頻率的分離而下降，當(dāng)ρ＞0.5之后，εω只有萬(wàn)分之幾個(gè) Δω，幅值相對(duì)誤差量級(jí)不超過(guò)10－4，相位誤差也只有零點(diǎn)幾度，且整周期采樣與否對(duì)幅值和頻率的精度影響不大(但成分1是否為整周期采樣對(duì)成分2的精度影響比較大)。

對(duì)影響誤差的各種因素比較發(fā)現(xiàn)：對(duì)誤差影響最大的因素是FFT的矩形積分法對(duì)精確積分接近程度［12］，特別是窗函數(shù)兩端權(quán)不為零，而數(shù)值積分又沒(méi)有計(jì)入這個(gè)值的影響。

在成分1的誤差曲線上標(biāo)有X的點(diǎn)，該點(diǎn)附近的誤差很大，其原因是式(12)中的Δ→0，因此數(shù)值上出現(xiàn)奇異。在校正實(shí)施中必須檢查這個(gè)值，以確定結(jié)果的可靠性。如果不考慮這一點(diǎn)，根據(jù)誤差隨頻率差而總體下降的趨勢(shì)，式(12)也適用于CSM的頻譜校正。

弱幅成分的參數(shù)識(shí)別誤差較大，但當(dāng)ρ=0.5～2(即若兩成分頻率間隔界于 0.5Δω ～2.0Δω)，則頻率誤差不超過(guò)0.1Δω，幅值相對(duì)誤差低于5%，相位誤差在幾度范圍之內(nèi)，大體可滿足工程精度的要求。

但是隨ρ繼續(xù)增大，弱幅成分的參數(shù)誤差趨勢(shì)并非繼續(xù)降低。這是因?yàn)殡Sρ增大，弱幅成分對(duì)強(qiáng)幅成分主峰附近的頻譜貢獻(xiàn)越小，因而利用強(qiáng)幅成分主峰周?chē)V線估計(jì)的弱幅信號(hào)參數(shù)的精度越來(lái)越差，尤其是奇點(diǎn)(圖中標(biāo)有X的點(diǎn))頻度也加增大。因而對(duì)于弱幅成分，如果ρ＞2，不如采用CSM校正公式更為可靠。

弱幅成分的頻率估計(jì)奇點(diǎn)往往由Δ=0所造成，但式(13)的Wp(ω2－ωa)也會(huì)導(dǎo)致幅值和相位的估計(jì)奇點(diǎn)。

5 結(jié)論

分析了單頻率模型的頻譜存在簡(jiǎn)單校正公式的原因。原因是常用窗函數(shù)的譜函數(shù)可以分解為超越函數(shù)與有理分式的乘積，前者對(duì)離散譜線間隔有周期性。利用這個(gè)特性很容易建立單頻率模型的加對(duì)稱(chēng)窗的頻譜校正公式。將該特性用于雙頻率模型(DFM)，也可以得到對(duì)應(yīng)的頻率方程，但是其關(guān)系非常復(fù)雜。除了加矩形窗外，其他窗情形涉及到高于2次的多項(xiàng)式方程，因此顯示校正公式可能不具有實(shí)用性或根本不存在。

對(duì)于DFM加矩形窗，頻率校正量滿足二次多項(xiàng)式方程，因此可以寫(xiě)出其顯式解，但是其復(fù)雜性遠(yuǎn)遠(yuǎn)超出單頻率模型的校正公式。采用含幅值差異較大的兩個(gè)頻率成分的信號(hào)，對(duì)這組顯式校正公式進(jìn)行了考核。結(jié)果表明，對(duì)于強(qiáng)幅成分，在工程精度內(nèi)可以認(rèn)為無(wú)誤差，而對(duì)弱幅成分，當(dāng)兩頻率間隔在0.5～2個(gè)頻率分辨率范圍之內(nèi)，誤差在工程精度上可以容忍。如果頻率相距比較大，還是各自采用單頻率模型為宜。

該方法的噪聲特性需要進(jìn)一步探討。

［1］沈玉娣，劉 雄，趙振毅.機(jī)械故障診斷－FORTRAN源程序匯編［M］.西安：西安交通大學(xué)出版社，1990.

［2］黃迪山.FFT相位誤差分析及實(shí)用校正方法［J］.振動(dòng)工程學(xué)報(bào)，1994，7(2)：185－189.

［3］謝 明，丁 康.頻譜分析的校正方法［J］.振動(dòng)工程學(xué)報(bào)，1994，7(2)：172－179.

［4］丁 康，謝 明.離散頻譜三點(diǎn)卷積幅值校正法的誤差分析［J］.振動(dòng)工程學(xué)報(bào)，1996，9(1)：92－98.

［5］丁 康，鐘舜聰，朱小勇.離散頻譜相位差校正方法研究［J］.振動(dòng)與沖擊，2001，20(2)：52－55.

［6］朱利民，鐘秉林，黃 仁.離散頻譜多點(diǎn)卷積幅值修正法的理論分析［J］.振動(dòng)工程學(xué)報(bào)，1999，12(1)：120－126.

［7］丁 康，張曉飛.頻譜校正理論的發(fā)展［J］.振動(dòng)工程學(xué)報(bào)，2000，13(1)：14－22.

［8］段虎明，秦樹(shù)人，李 寧.離散頻譜的校正方法綜述［J］.振動(dòng)與沖擊，2007，26(11)：138－144.

［9］謝 明，丁 康.兩個(gè)密集頻率成分重疊頻譜的校正法［J］.振動(dòng)工程學(xué)報(bào)，1999，12(1)：109－114.

［10］陳奎孚，張森文.利用三條譜線計(jì)算頻率緊鄰的兩個(gè)成分的參數(shù)［J］.振動(dòng)工程學(xué)報(bào)，2004，17(2)：153－158.

［11］謝 明，丁 康，莫克斌.頻譜校正時(shí)譜線干涉的影響及判定方法［J］.振動(dòng)工程學(xué)報(bào)，1998，11(1)：22－28.

［12］陳奎孚，焦群英，高小榕.改善FFT精度的2種算法比較［J］.中國(guó)農(nóng)業(yè)大學(xué)學(xué)報(bào)，1996，1(6)：74－79.