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

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

（第二炮兵工程學院，陜西 西安 710025）

1 引 言

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

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

2 考慮負頻響應的頻譜模型

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

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

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

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

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

3 算法推導

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

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

構造運算：

并且記：

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

算例1：取窗函數1，2，3進行研究，利用各自窗函數幅值恢復系數對自身進行補償。

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

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

整理可得：

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

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

當W1（f）=Wr（f），W2（f）=ahnWhn（f），

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

當W1（f）=Wr（f），W2（f）=abWb（f），

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

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

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

4 仿真研究與結果討論

4.1 仿真參數與計算條件

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

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

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

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

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

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

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

圖5 零頻率頻譜值曲線

4.2 結果與討論

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

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

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

（4）圖4在本質上是圖1對應研究的局部細化和延伸。可以發現當CiR在區間0.5～2.67時，該文測頻新方法誤差小于10-4數量級FFT變換分辨率，可以實現頻率的高精度測量（完全滿足測量電網諧波頻率變化的需要[7-9]），在其他CiR大于2.67區間上雖然誤差比較大，但可通過減少分析點數將CiR轉化至區間 0.5～2.67。

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

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

圖6 新算法的測頻誤差包絡

5 結束語

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

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

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

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

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

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

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

[4] 陳奎孚，王建立，張森文.低頻成分的頻譜校正[J].振動工程學報，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].北京：電子工業出版社，2006.