基于變窗函數(shù)方法的頻率檢測技術(shù)研究

文超斌，李世平，劉如峰，宋 兵

（第二炮兵工程學(xué)院，陜西 西安 710025）

1 引 言

獲取信號頻率作為信號分析的主要內(nèi)容之一，在工程應(yīng)用中占有很重要的地位。目前國內(nèi)外對多頻信號通過頻譜校正進(jìn)行精確測量的方法很多，如比值校正法，相位差校正法等[1-6]。這些方法對解析單頻信號或頻率成分間隔較大的情景都具有優(yōu)良的特性，但由于沒有考慮負(fù)頻譜的影響，所以應(yīng)用在正負(fù)頻譜發(fā)生嚴(yán)重干涉的低頻信號譜中時，其校正精度明顯下降，甚至失去校正作用。

頻率的高低與記錄的長短有關(guān)，信號頻率的高低是相對的概念，如果分析樣本內(nèi)含有足夠多的波動周期數(shù)則可視為高頻，這時就可以忽略負(fù)頻譜的影響。然而在很多種情形下，采集樣本可能不滿足這個條件，例如工程實踐中往往要進(jìn)行信號實時處理從而要求被處理信號足夠短。另外，很多情況下已采集到的數(shù)據(jù)不夠長，但再次采集成本又很高，這時信號正負(fù)頻譜發(fā)生了嚴(yán)重的干涉，常用校正方法將達(dá)不到精度的要求，甚至不能使用。因此，有必要進(jìn)一步研究建立系統(tǒng)的考慮負(fù)頻響應(yīng)的高精度多頻率測量技術(shù)，文獻(xiàn)[3-5]給出了考慮負(fù)頻譜影響的幾種方法，但是這些方法的測頻精度受信號中直流分量的影響較大。該文在考慮負(fù)頻響應(yīng)的情況下建立信號頻譜模型，多個窗函數(shù)組合應(yīng)用，總結(jié)出一套顯式測頻的新方法。仿真結(jié)果表明，這個方法比簡單比值法精度更高、應(yīng)用范圍更廣、實用性更強。

2 考慮負(fù)頻響應(yīng)的頻譜模型

帶直流偏移的正弦波模型x（t）如下：

式中：f0、A、φ0和B——信號的頻率、幅值、初相位和直流偏移量4個參數(shù)。

假定式（1）在[0，T]上截斷。采樣樣本內(nèi)包含的信號周期個數(shù)可簡記為 CiR（Cycles in Record）[5]，其值等于f0與Δf的比，其中Δf表示信號FFT變換頻譜的頻率分辨。設(shè)w（t）為定義在區(qū)間[-T/2，T/2]上的對稱窗函數(shù)，其頻域表達(dá)式為W（f），把該窗函數(shù)平移至區(qū)間[0，T]得窗函數(shù)wT（t）。用窗函數(shù)wT（t）對正弦波模型信號x（t）進(jìn)行截斷得截斷信號xT（t）。則該截斷信號的傅里葉變換頻譜函數(shù)XT（f）為頻率f的連續(xù)函數(shù)，即：

由傅里葉變換時移特性有：

所以由式（2）、式（3）、式（4）并且結(jié)合傅里葉變換頻移特性有：

3 算法推導(dǎo)

考慮式（5）截斷信號的傅里葉變換頻譜函數(shù)XT（f）在零頻率處的頻譜值為：

用任意3個對稱窗函數(shù)如上面敘述操作：分別平移后對正弦波信號模型x（t）截斷，并且分別求出其對應(yīng)的3個傅里葉變換頻譜函數(shù)XTi（f），考慮函數(shù)XTi（f）在零頻率處的頻譜值，記為XTi（0）（其中i=1，2，3）。

構(gòu)造運算：

并且記：

研究發(fā)現(xiàn)，若利用窗函數(shù)幅值恢復(fù)系數(shù)在FFT變換前對所涉及到的窗函數(shù)自身進(jìn)行補償，確保截斷信號的傅里葉變換頻譜函數(shù)對應(yīng)頻譜在零頻率處幅值相等，這樣式（7）構(gòu)造的運算：做差消去了直流分量對測頻精度的影響；相除可使式（7）右邊運算中涉及到的復(fù)雜非線性部分對消，最終得到關(guān)于要求頻率的函數(shù)。

算例1：取窗函數(shù)1，2，3進(jìn)行研究，利用各自窗函數(shù)幅值恢復(fù)系數(shù)對自身進(jìn)行補償。

將式（9）3個新的窗函數(shù)分別代入式（6），求出此3個新的窗函數(shù)分別由區(qū)間[-T/2，T/2]平移至區(qū)間[0，T]，而后對正弦波信號模型 x（t）進(jìn)行截斷，所得3個不同截斷信號的傅里葉變換頻譜函數(shù)在零頻率處的頻譜值：

其中 i=1，2，3，于是由式（9）、式（10）并且結(jié)合式（7）構(gòu)造運算得：

整理可得：

另外，對正弦波截斷信號xT（t）進(jìn)行采樣得離散信號xT（n），該離散信號進(jìn)行FFT變換可得截斷信號xT（t）對應(yīng)的離散頻譜在頻率fi處的頻譜值為：

（其中 i=0，1，2，3，…，N）可見由計算機對采樣信號做FFT變換可直接求出信號離散頻譜在零頻率處的頻譜值，而這正是前面敘述中截斷信號的傅里葉變換頻譜函數(shù)XT（f）在零頻率處的頻譜值XT（0）。于是利用式（9）中經(jīng)過補償?shù)?個窗函數(shù)分別對正弦波模型信號x（t）截斷，采樣得其對應(yīng)的離散信號，而后進(jìn)行FFT變換可得到3個不同截斷信號的傅里葉變換頻譜函數(shù)XTi（f）在零頻率處的頻譜值XTi（0）（其中i=1，2，3），進(jìn)而由式（8）可求出K。這樣由式（11）即可以解出正弦波信號的頻率f0：

當(dāng)W1（f）=Wr（f），W2（f）=ahnWhn（f），

W3（f）=abWb（f）時，同理得：

當(dāng)W1（f）=Wr（f），W2（f）=abWb（f），

W3（f）=ahnWhn（f）時，同理得：

算例2：取窗函數(shù)2，3，4進(jìn)行研究時，同樣方法可得對應(yīng)于3種窗函數(shù)不同的組合順序頻率計算公式分別為：

算例3：取窗函數(shù)1，3，4進(jìn)行研究時，同樣方法可得對應(yīng)于3種窗函數(shù)不同的組合順序頻率計算公式分別為：

4 仿真研究與結(jié)果討論

4.1 仿真參數(shù)與計算條件

按照式（1）正弦波信號模型生成原始數(shù)據(jù)，這里只考慮無噪聲干擾的情形。信號直流偏移量B=1V，幅值A(chǔ)=1V，初相位φ0=0。采樣樣本大小N=2000，采樣頻率fS=50000Hz，信號FFT變換頻譜的頻率分辨Δf=fS/N=25Hz。

為全面反映該文所述新方法測頻的效能，這里用以上總結(jié)出的各個測頻計算公式，分別在信號頻率變動化時，求出各個仿真信號的頻率歸一化絕對誤差，并繪出了相應(yīng)的曲線圖，被測低頻信號頻率f0從0.1Δf掃描到10Δf，掃描步長為0.1Δf（電網(wǎng)中的諧波頻率變化通常為45～55 Hz對應(yīng)這里為頻率f0從1.8Δf掃描到2.2Δf[7-9]）。也就是仿真信號采樣樣本內(nèi)包含的信號周期個數(shù)CiR從0.1線性增加到10，步長為 0.1。

圖1～圖3中給出了以上算法推導(dǎo)部分3種不同的算例情況，分別用各自算例3個不同公式計算頻率的歸一化絕對誤差曲線。第一個算例（圖1）取窗函數(shù) 1，2，3，分別用式（11）、式（13）、式（14）進(jìn)行計算；第二個算例（圖2）取窗函數(shù)2，3，4分別用式（15）、式（16）、式（17）進(jìn)行計算；第三個算例（圖 3）取窗函數(shù) 1，3，4 分別用式（18）、式（19）、式（20）進(jìn)行計算；圖4畫出了當(dāng)CiR在區(qū)間0.5～2.67變化，步長變?yōu)?0.02Δf，取窗函數(shù) 1，2，3（算例 1）時，分別用3 個不同式（11）、式（13）、式（14）計算得到的頻率歸一化絕對誤差曲線；而且每個算例3個不同的頻率計算公式計算出的結(jié)果都用不同的形狀表示。圖5給出了正弦波信號模型x（t）在直流偏移量B=0，其他信號參數(shù)如前述不變情況下，仍在變頻率步長0.1Δf情況下對每一個信號分別加該文分析的4種窗函數(shù)后，不同截斷信號的傅里葉變換頻譜函數(shù)在零頻率處的頻譜值變化曲線。5幅圖橫坐標(biāo)均為記錄樣本中所包含的信號周期數(shù)（CiR）。

圖1 算例1測頻誤差曲線（步長0.1Δf）

圖2 算例2測頻誤差曲線（步長0.1Δf）

圖3 算例3測頻誤差曲線（步長0.1Δf）

圖4 算例1測頻誤差曲線（步長0.02Δf）

圖5 零頻率頻譜值曲線

4.2 結(jié)果與討論

（1）對照圖1～圖3可以看出，算例2的3個頻率計算公式，算例3的頻率計算式（18）、式（20）測頻誤差太大，大部分頻率點測頻誤差都落在0.4倍的頻率分辨率以外，所以這幾個校正公式不能使用；算例1的3個頻率計算公式以及算例3的頻率計算公式（19）均能得到比較好的測頻結(jié)果。

（2）由圖1、圖3可以看出，隨著采樣樣本內(nèi)包含的信號周期個數(shù)CiR變大測頻誤差控制在0.1倍頻率分辨率的點越來越少，而且誤差比較大的點在圖5曲線過零點對應(yīng)位置成周期出現(xiàn)。原因可從圖5得出，如果對信號分別加窗函數(shù)1，2，3，4進(jìn)行FFT變換，對應(yīng)于同一個CiR，四個離散頻譜在零頻率點對應(yīng)頻譜值相差不大，就會造成求K時式（8）分母接近零，出現(xiàn)奇點造成計算誤差變大。一方面隨著采樣樣本內(nèi)包含的信號周期個數(shù)CiR變大，圖5中對應(yīng)于同一個CiR 4個離散頻譜在零頻率點對應(yīng)頻譜值逐漸衰減趨于同樣大小，從而引起這一奇點效應(yīng)；另一方面隨著圖5零頻率點對應(yīng)頻譜值的過零點呈周期出現(xiàn)，同樣也會引起這一奇點效應(yīng)，造成計算誤差也跟著呈周期變大?？梢姲l(fā)生奇點效應(yīng)是引起這一算法測頻精度下降最主要的原因，但由于測頻誤差不是隨周期數(shù)（CiR）單調(diào)變化，所以適時調(diào)節(jié)分析點數(shù)運用該新算法測頻就能夠很好地避免這里的奇點效應(yīng)。

（3）圖1、圖4雖然是分別用算例1的3個頻率計算公式計算出的頻率歸一化絕對誤差，但大部分點是重疊的，原因是上文總結(jié)出同一個算例中的3個測頻計算公式線性相關(guān)，計算結(jié)果只有計算機處理數(shù)據(jù)精度上差別，并無原理上大的不同。

（4）圖4在本質(zhì)上是圖1對應(yīng)研究的局部細(xì)化和延伸?？梢园l(fā)現(xiàn)當(dāng)CiR在區(qū)間0.5～2.67時，該文測頻新方法誤差小于10-4數(shù)量級FFT變換分辨率，可以實現(xiàn)頻率的高精度測量（完全滿足測量電網(wǎng)諧波頻率變化的需要[7-9]），在其他CiR大于2.67區(qū)間上雖然誤差比較大，但可通過減少分析點數(shù)將CiR轉(zhuǎn)化至區(qū)間 0.5～2.67。

（5）采樣密度對精度的影響不顯著。采樣率首先要滿足采樣定理的要求，即每個波動周期內(nèi)至少采到兩點。由于固定了采樣點數(shù)N，所以CiR越高，每個周期內(nèi)采集的點數(shù)越少。但是圖1、圖3表明采樣密度對校正精度的影響不具主導(dǎo)性。

（6）為了更加系統(tǒng)地考察新方法測頻的誤差特性，重新選取被測正弦波信號模型頻率f0從0.1Δf掃描到 10Δf，掃描步長為 0.02Δf，其他信號參數(shù)同4.1節(jié)開始第一段敘述不變。在變頻率情況下考慮初相位角φ0變化時，運用算例1的頻率計算公式（11）求出每個頻率的180個不同的初相位角情況下頻率f0的歸一化絕對誤差，而后找到其中的最大值，考察這個最大值隨CiR變化的趨勢，繪制誤差曲線包絡(luò)，如圖6所示（a圖是CiR在區(qū)間0.5～2.67變化時頻率誤差特性圖，b圖畫出了CiR僅在區(qū)間0.5～2.67變化時情況）。該圖表明信號初相位角無論如何變化，運用該文分析的測頻方法均能得到比較高的精度，再次印證了以上總結(jié)出結(jié)論的正確性，說明了該文新方法測頻的有效性。

圖6 新算法的測頻誤差包絡(luò)

5 結(jié)束語

普通頻譜校正方法由于沒有考慮負(fù)頻譜成分的干涉影響，難以應(yīng)用于短記錄的頻譜校正測頻中。該文給出了一種基于變窗函數(shù)方法的頻譜校正的顯式測頻新方法，采用仿真手段對提出的方法進(jìn)行了考核。仿真結(jié)果表明：

（1）采用窗函數(shù)1，2，3進(jìn)行組合的各類頻率計算公式均可以得到比較良好的測頻結(jié)果，并且不用考慮信號中直流分量的影響。當(dāng)CiR在區(qū)間0.5～2.67時頻率誤差特別小，在10-4FFT變換分辨率數(shù)量級上，并且相位變化和采樣密度變化不對測頻精度產(chǎn)生較大影響。

（2）影響該新方法測頻精度最重要的因素是計算公式中有可能出現(xiàn)奇點。工程應(yīng)用中應(yīng)盡量通過選取對應(yīng)同一頻率點頻譜差距較大的一組窗函數(shù)和調(diào)節(jié)分析點數(shù)，使得采樣樣本內(nèi)包含的信號周期個數(shù)CiR在對應(yīng)于測頻誤差最小的范圍內(nèi)，來很好地回避這一問題。

該文的測頻方法只研究了變窗函數(shù)對象為矩形窗、Hanning窗、Blanckman窗、Hamming窗情形，其他窗函數(shù)情形的測頻公式有待進(jìn)一步發(fā)展。此外，在噪聲干擾情形下的測頻誤差的統(tǒng)計特性也有待進(jìn)一步研究。

[1] 丁 康，謝 明，楊志堅.離散頻譜校正理論與技術(shù)[M].北京：科學(xué)出版社，2007.

[2] 斯托伊卡.現(xiàn)代信號譜分析[M].吳任彪，韓 萍，馮 清，等，譯.北京：電子工業(yè)出版社，2007.

[3] 陳奎孚，張森文，郭幸福.消除負(fù)頻率影響的頻譜校正[J].機械強度，2004，26（1）：25-28.

[4] 陳奎孚，王建立，張森文.低頻成分的頻譜校正[J].振動工程學(xué)報，2008，21（1）：289-292.

[5] 陳奎孚，王建立，張森文.短記錄加漢寧窗的頻譜校正[J].振動與沖擊，2008，27（4）：49-51.

[6] 段虎明，秦樹人，李 寧.離散頻譜的校正方法綜述[J].振動與沖擊，2007，26（11）：138-145.

[7] Radil T，Ramos P M，Serra A C.Frequency estimation of power system signals using a new spectrum leakage correction algorithm[J].IEEE MTC.Victoria，BC，Canada，2008，25（3）：2161-2166.

[8] Radil T，Ramos P M.New spectrum leakage correction algorithm for frequency estimation ofpower system signals [J].IEEE Trans.Instrum.Meas，2008，57（5）：1670-1679.

[9] Agrei D.Power system frequency estimation in the shortened measurement time[C]∥IEEE Instrumentation and Measurement Technology Conference：Budapest，Hungary，2007，21（23）：1094-1098.

[10]王正林，劉 明.精通Matlab7[M].北京：電子工業(yè)出版社，2006.