戴 彥

(唐山學(xué)院 智能與信息工程學(xué)院，河北 唐山 063020)

一種改進(jìn)的FFT離散頻譜相位差加權(quán)校正算法

戴 彥

(唐山學(xué)院 智能與信息工程學(xué)院，河北 唐山 063020)

為解決用傳統(tǒng)相位差法進(jìn)行頻譜校正精度不高的問(wèn)題，提出了一種改進(jìn)的FFT離散頻譜相位差加權(quán)校正算法。該算法是通過(guò)將相位差中最大值校正誤差與次大值校正誤差按一定比例進(jìn)行加權(quán)來(lái)獲得頻譜校正結(jié)果。對(duì)此的仿真結(jié)果表明，該算法與傳統(tǒng)相位差法相比，在頻譜最大值相同的情況下，頻率校正誤差的頻率分辨率提高了近90%，提高了頻率校正精度。

頻譜校正；相位差法；加權(quán)校正算法

0 引言

快速傅里葉變化(FFT)算法是1965年由Cooely-Tukey在《計(jì)算數(shù)學(xué)》雜志上首次發(fā)表的[1]。自此，F(xiàn)FT和頻譜分析很快發(fā)展成為無(wú)線電通信、信息圖像處理、自動(dòng)控制、機(jī)械設(shè)備故障診斷和振動(dòng)分析等多種學(xué)科重要的理論基礎(chǔ)[2]。由于計(jì)算機(jī)處理樣本個(gè)數(shù)有限，F(xiàn)FT和頻譜分析也只能在有限區(qū)間內(nèi)進(jìn)行，這樣不可避免地會(huì)產(chǎn)生能量泄漏，造成譜峰值變小，精度降低[3]。目前國(guó)內(nèi)外關(guān)于對(duì)幅值譜或功率譜進(jìn)行校正的方法有四種：離散頻譜三點(diǎn)卷積幅值校正法[4]、對(duì)幅值譜進(jìn)行校正的比值法[5]、FFT譜連續(xù)細(xì)化分析的傅里葉變換法[6]、相位差法[7]。其中相位差法是利用兩段序列在FFT最大譜線處的相位差來(lái)進(jìn)行的一種離散頻譜校正方法。本文針對(duì)相位差法頻譜校正中頻率受FFT頻譜量化的影響而精度不高的問(wèn)題，提出了將最大值誤差與次大值誤差進(jìn)行加權(quán)校正的一種改進(jìn)的離散頻譜相位差校正算法。

1 加權(quán)校正算法分析

1.1 相位差校正法

設(shè)有一周期信x(t)=Acos (2πf0t+θ)，其中f0，θ，A為信號(hào)的真實(shí)頻率、相位和幅值，作傅立葉變換，結(jié)果為

(1)

式中f為FFT后的頻率變量。加窗后的諧波信號(hào)x(t)wT(t)的傅立葉變換可根據(jù)卷積定理表示為

(2)

設(shè)f1為信號(hào)峰值處的頻率，由此可得，加窗后在信號(hào)峰值處的相位為

θ1=θ-πTΔf，

(3)

其中Δf=f1-f0。

將連續(xù)信號(hào)x(t)向前平移u1T得x0(t)，窗長(zhǎng)u2T的對(duì)稱窗函數(shù)在時(shí)間上平移u2T/2得到wT0(t)。對(duì)x0(t)wT0(t)作FFT(或DFT)分析，設(shè)f2為信號(hào)峰值處的頻率，在信號(hào)峰值處的相位為

θ2=θ-πu2T(f2-f0)+2π(f1-Δf)u1T，

(4)

其中f0=f1-Δf。

相位差由式(3)減去式(4)得到式(5)，

θ=πu2T(f2-f1)-2πf1u1T+(u2+2u1-1)πTΔf。

(5)

頻率修正量為

(6)

式中u2+2u1-1≠0。

1.2 加權(quán)校正算法可行性分析

當(dāng)時(shí)域平移相對(duì)量μ1分別為0.3，0.5，0.7，矩形窗長(zhǎng)相對(duì)量μ2分別為0.3，0.5，0.7，1時(shí)，最大值校正誤差如圖1，圖2所示。

圖1 時(shí)域平移對(duì)校正誤差的影響

圖2 窗長(zhǎng)相對(duì)量對(duì)校正誤差的影響

由圖1和圖2可知，當(dāng)u1=0.5，u2=1時(shí)，校正誤差接近為0，此時(shí)校正精度最高。因此，式(6)中u1，u2，T均為固定值，而f1，f2在信號(hào)頻率一定的情況下也為常量，這樣Δθ就只與Δf相關(guān)了。

假設(shè)最大值在次大值的左方，對(duì)第一段采樣信號(hào)進(jìn)行FFT處理后，設(shè)所取次大值對(duì)應(yīng)的頻率為f1′，顯然f1′=f1+fs/N，頻率校正量Δf′=f1′-f0=f1′-f1+Δf。經(jīng)過(guò)時(shí)域平移、加窗后再進(jìn)行FFT處理，設(shè)所取最大值頻率為f2′，則f2′=f2+fs/N，fs是采樣頻率。與使用頻譜最大值校正時(shí)的相位差公式不同，使用次大值校正時(shí)：

Δφ′=φ1′-φ2′=

πμ2T(f2′-f1′)-2πμ1Tf1′+(2μ1+μ2-1)πTΔf′=

Δφ+(μ2-1)π。

(7)

由式(7)可以推出：

(8)

由次大值相位差得到的頻譜校正量Δf′與由最大值相位差得到的頻譜校正量Δf、頻率分辨率fs/N及參數(shù)u1，u2相關(guān)，而頻率分辨率fs/N，u1，u2為固定值，因此次大值相位差的頻譜校正量Δf′只與Δf相關(guān)。其中，N為采樣點(diǎn)個(gè)數(shù)。

在Matlab環(huán)境下生成信號(hào)x(t)=2cos (2πf0n)，采樣頻率fs為102.4 kHz，采樣點(diǎn)數(shù)N=1 024，頻率分辨率為100 Hz，時(shí)域平移量為512點(diǎn)，選用矩形窗，f0為實(shí)際頻率，在15～15.1 kHz范圍內(nèi)，每隔2.5 Hz，使用傳統(tǒng)相位差法分別進(jìn)行無(wú)噪和加高斯白噪聲的仿真分析，信噪比是50 dB。校正結(jié)果如圖3，圖4所示。

圖3 無(wú)噪聲校正結(jié)果

圖4 有噪聲校正結(jié)果

分析圖3，圖4可知，使用次大值得到的頻率校正誤差與使用最大值得到的頻率校正誤差符號(hào)相反，最大值和次大值的校正誤差在零點(diǎn)處相同，最大值校正誤差絕對(duì)值比次大值校正誤差絕對(duì)值要小。因此，對(duì)最大值及次大值的校正結(jié)果進(jìn)行加權(quán)處理，可以得到更加精確的校正結(jié)果。

1.3 加權(quán)校正公式的確定

根據(jù)比值校正法的原理，將頻譜次大值與頻譜最大值的比值作為頻率校正中間量，并參照?qǐng)D3，將無(wú)噪聲情況下的次大值誤差絕對(duì)值與最大值誤差絕對(duì)值相比得到一個(gè)變量，這個(gè)變量的曲線如圖5所示。

圖5 加權(quán)值確定

在比值法中頻譜最大值與頻譜次大值的比值與頻率校正量Δk相關(guān)，頻率校正量Δk的取值范圍為[-0.5，0.5]，在一個(gè)頻率分辨區(qū)間內(nèi)，即在圖4中(15 050，0)點(diǎn)之前，隨頻率的增大Δk從0增大到0.5，并呈線性；(15 050，0)點(diǎn)之后，隨頻率的增大Δk從-0.5增大到0，也呈線性。由于(15 050，0)點(diǎn)之后的Δk與其之前的Δk對(duì)稱，只需取絕對(duì)值即可以將|Δk|作為一個(gè)變量。注意到圖5中無(wú)噪聲情況下的次大值誤差絕對(duì)值與最大值誤差絕對(duì)值比值曲線的規(guī)律類似于1/x函數(shù)，結(jié)合Δk自身的規(guī)律，現(xiàn)構(gòu)造一個(gè)變量Δ=(1-|Δk|)/|Δk|，Δ的分布規(guī)律如圖5所示?？梢?jiàn)Δ曲線基本上與次大值校正誤差絕對(duì)值和最大值校正誤差絕對(duì)值的比值吻合，所以確定Δ為加權(quán)校正量，加權(quán)校正公式為[4]：

(9)

圖5中(1-x)/x擬合曲線的橫坐標(biāo)表示實(shí)際頻率的譜線號(hào)與頻譜最大值所在處譜線號(hào)之差的絕對(duì)值。

2 仿真分析

將最大值校正誤差與次大值校正誤差進(jìn)行加權(quán)，得到的校正結(jié)果如圖6，圖7所示。

圖6 無(wú)噪聲加權(quán)校正結(jié)果

圖7 有噪聲加權(quán)校正結(jié)果

有噪聲情況下得到的加權(quán)校正數(shù)據(jù)為15次運(yùn)算后的平均結(jié)果，圖8是這些數(shù)據(jù)對(duì)應(yīng)的方差。

圖8 頻率校正結(jié)果方差分布圖

由圖6，圖7，圖8，可得：

(1)無(wú)噪聲情況下，加權(quán)校正后的頻率校正誤差幾乎為零；有噪聲情況下，加權(quán)校正誤差比單純使用最大值或次大值校正得到的誤差小很多。

(2)mfs/N在左右5%的范圍內(nèi)(即頻譜線附近)加權(quán)后有所改善，但效果仍不理想，主要表現(xiàn)是校正結(jié)果誤差偏大且方差也較大。這是由于實(shí)際頻譜最大值很接近譜線，使得頻譜譜線最大值處能量很大，相應(yīng)地造成頻譜譜線次大值處的能量很小。由于噪聲的干擾，使得主瓣之外的譜線能量(大部分情況都是第一旁瓣內(nèi)靠近譜線最大值的那條譜線能量)高于主瓣內(nèi)譜線次大值處的能量，因此易造成錯(cuò)判次大值的情況。由此可見(jiàn)該方法的使用范圍有一定的局限性，即mfs/N在左右5%的范圍內(nèi)校正不適用。

3 結(jié)論

本文提出了基于FFT相位的頻譜分析原理，從理論上說(shuō)明了使用頻譜次大值進(jìn)行相位差校正的可行性，同時(shí)通過(guò)對(duì)仿真結(jié)果的分析，找到了一種使用頻譜最大值、頻譜次大值進(jìn)行相位差加權(quán)聯(lián)合校正的方法。在采樣頻率為102.4 kHz、采樣點(diǎn)數(shù)為1 024點(diǎn)、信噪比為50 dB情況下，使用該方法得到的頻率校正誤差在頻率分辨率的0.012%以內(nèi)，而在頻譜最大值相同的情況下使用傳統(tǒng)相位差法得到的頻率校正誤差在頻率分辨率的0.115%以內(nèi)。因此，使用該加權(quán)校正算法，在只增加25%運(yùn)算量的前提下，頻率分辨率提高了近90%，大大提高了頻率校正的精度。

[1] Cooley J W， Tukey J W. An algorithm for the ma-chine calculation of complex Fourier series[J]. Mathe-matics of Computation，1965，19(90)：297-301.[2] 楊志堅(jiān)，丁康.調(diào)制FFT及其在離散頻譜校正技術(shù)中的應(yīng)用[J].振動(dòng)工程學(xué)報(bào)，2009，22(6)：632-637.

[3] 丁康，鐘順聰.通用的離散頻譜相位差校正方法[J].電子學(xué)報(bào)，2003，31(1)：142-145.

[4] 丁康，謝明. Error analysis for amplitude correction method using convolution of three points in discrete spectrum[J].振動(dòng)工程學(xué)報(bào)，1996，9(1)：92-98.

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

[6] 丁康，朱文英，楊志堅(jiān)，等.FFT+FT離散頻譜校正法參數(shù)估計(jì)精度[J].機(jī)械工程學(xué)報(bào)，2010，46(7)：68-73.

[7] 譚思煒，任志良，孫常存.全相位FFT相位差頻譜校正法改進(jìn)[J].系統(tǒng)工程與電子技術(shù)，2013，35(1)：34-39.

[8] 侯慶文，陳先中.改進(jìn)的FMCW信號(hào)加權(quán)補(bǔ)償校正相位差法[J].儀器儀表學(xué)報(bào)，2010，31(4)：721-726.

(責(zé)任編校：李秀榮)

An Improved FFT Weighted Correction Algorithm for the Discrete Spectrum Phase Difference

DAI Yan

(School of Intelligence and Information Engineering，Tangshan College， Tangshan 063000， China)

The author of this paper proposes an improved FFT weighted correction algorithm for the discrete spectrum phase difference to overcome the shortcoming of inaccuracy of the traditional spectrum correction for phase difference. By this algorithm， the spectral correction can be achieved through weighting the biggest calibration error and the second biggest calibration error by a certain percentage. The simulation results show that in the case of the maximum spectrum， the frequency resolution of the frequency calibration error is increased by 90 percent with the proposed algorithm.

spectrum correction； phase difference； weighted correction algorithm

TN911.7

A

1672-349X(2015)06-0031-03

10.16160/j.cnki.tsxyxb.2015.06.012