適用于非整數次冪的高精度混合基FFT諧波測量算法

王保帥 肖 勇 胡珊珊 趙 云

（南方電網科學研究院有限責任公司 廣州 510663）

0 引言

隨著環境和能源問題的日益突出，采用太陽能、風能等新能源分布式發電取得了很大的發展[1]。但其中的電力電子裝置等非線性設備會向電網注入諧波[2]，導致電網中的諧波污染日益嚴重，給電力系統的安全穩定運行造成了不可忽視的影響。對于電網中電壓電流等電參數的準確量測是諧波分析治理的關鍵所在。快速傅里葉變換（Fast Fourier Transform, FFT）因易于嵌入實現的特點，是目前諧波測量分析的最常用算法[3-5]。

國際標準IEC61850-9-2LE中規定，數字化電能計量系統的電參數采樣率為256fr或80fr（fr為工頻50Hz）[6]。在國際標準IEC61000-4-7中規定在工頻50Hz的電網信號采集分析中，進行一次快速傅里葉變換運算的時間窗必須是10個信號周期。這就產生了兩個方面的影響：一方面，實際的電力系統頻率總是圍繞在fr附近波動，同時采樣率固定，電壓電流的同步采樣則得不到保障，在非同步采樣下就會產生頻譜泄露；另一方面，按照標準規定，采樣數據序列長度必定不是2的整數次冪，在進行諧波分析時，無法采用傳統的基2快速傅里葉分析算法。

目前非整周期截斷和非同步采樣引起的頻譜泄露問題多采用加窗插值FFT算法解決，國內外學者先后運用余弦窗[7]、三角窗[8]、卷積窗[9]、乘法窗[10]等窗函數對信號序列進行截斷，抑制長范圍泄露。在插值方面，學者們相繼提出雙譜線插值[11-12]、三譜線插值[13-14]、四譜線插值[15-16]、六譜線插值[17-18]等算法修正分析結果，抑制短范圍泄露。窗函數越來越復雜，插值譜線數目越來越多，雖然提升了諧波分析計算的準確度，但計算的復雜度也隨之提升，并且都沒有解決FFT對數據序列長度有所限制的問題。如應用較多的基2FFT算法，要求數據序列長度必須是2的整數次冪，加窗插值算法無法解決這個問題。當采樣數據序列不是2的整次冪時，可以采用離散傅里葉變換（Discrete Fourier Transform, DFT）運算分析，但是DFT的運算量過于龐大，不利于工程應用。已有文獻研究中通常采用補零的方法使數據序列長度滿足基2的要求，目前多數補零研究的文獻大多直接采用FFT進行計算，較少分析相應的影響。文獻[19]分析指出，補零前后的離散傅里葉變換結果不同，補零需要滿足一定的規則。文獻[20]分析了補零對三譜線插值算法修正公式的影響，以加漢寧窗為例分析研究得出補零后歸一化頻率、幅值和相位的修正公式。文獻[21-24]等探討研究了混合基FFT在電能質量分析和諧波分析中的應用，得出混合基FFT在非2整次冪數據序列分析中具有一定的應用價值。

針對諧波分析計算中采樣數據長度為非2整次冪的情況，本文重點分析了混合基FFT的頻譜分布情況，提出一種適用于非整數次冪的高精度混合基FFT諧波測量算法，闡述了其原理和校正過程。該算法的優越性在于對采樣數據序列長度沒有2整次冪的要求，同時改善了非同步采樣下混合基FFT的頻譜泄露問題。最后結合仿真及試驗分析，驗證了其可行性和有效性。

1 常規混合基FFT算法

當對采樣數據序列進行FFT分析計算時，序列數據的長度需要滿足一定的要求，如基2FFT算法要求數據長度必須是2的整數次冪，基3FFT算法要求數據長度必須是3的整數次冪等。當數據長度不滿足要求時，只能采用DFT算法、補零或混合基FFT算法對采樣數據序列進行分析。若對采樣數據序列直接進行DFT運算，則會導致運算量過于龐大。對于補零運算，在最壞的情況下，需要補加N/2-1個零，也浪費了大量計算資源[25-26]。

對于非2的整數次冪數據序列，只要序列數據長度N是復合數，就可以應用混合基FFT算法[21,27]對其進行分析。N點序列長度x(n)的DFT為

式中，WN為旋轉因子，。當序列數據長度N是2的整數次冪時，可以直接采用基2 FFT算法進行計算；當N不等于2的整數次冪時，N可以表示為兩個因子的乘積，即N=N1N2，此時

式中，n為N2進制數，nβ為其進位，nα為其末位；k為N1進制數，kβ為其進位，kα為其末位。將n和k代入式（1）得

利用旋轉因子的特性，可將式（3）化簡為

式中，X1(kα,nα)為以nβ為輸入、kα為輸出的N1點 DFT，；X2(kα,nβ)為以nα為輸入，kβ為輸出的N2點DFT，。從式（4）可以看出，在計算N點DFT時，可以將N拆分為N1與N2兩個因子的乘積，先計算N1點DFT，再計算N2點DFT即可得到結果。若N2仍可表示為多個因子的乘積，可繼續向下層分解。當N1或N2存在2的整數次冪時，就可以調用基2FFT算法進行快速計算。

2 高精度混合基FFT諧波測量算法

2.1 混合基FFT頻譜譜線分析

一般的諧波信號在時域中的表達式為

式中，m為諧波次數；H為諧波總次數；f0為基波頻率；Am、φm分別為第m次諧波的幅值和初相位。

按照國際標準規定，選取采樣10個周波共2 560點數據對式（5）所示信號進行混合基FFT運算。根據常規混合基FFT算法，首先將2 560點數據劃分為5組，每組512點數據，需要完成5組512點的基2FFT運算，然后對其結果進行512組5點DFT運算。常規混合基FFT運算的頻譜譜線分布如圖1所示。

由圖1可知，當基波頻率偏離50Hz發生變化時，常規混合基FFT算法產生了較大的頻譜泄露，使得諧波分析的準確度降低。為說明常規混合基FFT算法目標頻點譜線群的分布情況，對圖1所示的譜線群加窗，混合基FFT運算加窗后的第m次諧波頻譜為

圖1 常規混合基FFT運算后的諧波頻譜譜線分布Fig.1 Harmonic spectrum line distribution after conventional mixed-base FFT operation

式中，Δf為頻譜分辨率，；fs為采樣速率；km為目標頻點即第m次諧波的譜線號；W(ω)為余弦窗的頻譜函數。在加窗后的頻譜中，目標頻點處一般沒有離散譜線，依據目標頻點譜線、幅值最大譜線和幅值次最大譜線的分布情況，有兩種典型的譜線群分布情況，如圖2所示。如果采用常規混合基FFT算法對非2整次冪采樣數據序列進行分析，則只是用目標頻點附近幅值最大譜線的信息估算目標頻點的信息，忽略了目標頻點附近譜線群中包含的豐富諧波信息，勢必會產生較大的分析誤差。

圖2 諧波頻譜譜線群分布情況Fig.2 Distribution of harmonic spectral line groups

2.2 高精度混合基FFT測量算法

在圖2中，km是目標頻點的譜線號，k1、k2、k3、k4依次為目標頻點附近幅值按從大到小排列的四根譜線的譜線號。由圖2可以明顯看出，在對非2整次冪數據序列進行混合基加窗FFT運算后的頻譜中，目標頻點附近的4根譜線幅值都比較高，含有豐富的諧波信息。設這四條譜線按幅值從大到小的頻率為k1Δf、k2Δf、k3Δf、k4Δf，對應的幅值為y1、y2、y3、y4。在圖2a所示的情形下，目標頻點處幅值最大譜線在目標頻點的左側，有k3＜k1＜km＜k2＜k4，并且 0＜km-k1＜1。令α=km-k1- 0.5，則α∈(-0.5,+0.5)，引入變量

在窗函數的頻譜函數W(ω)中，令ω=得

假設

則由式（7）可得

利用多項式曲線擬合得到α=g(β)，從而得出

按照式（6）修正幅值為

其中，幅值修正式為。

由式（5）和式（6）可得到相位修正式為

式中，arg為求復數輻角函數。

同理，在圖2b所示的情況下，目標頻點處幅值最大譜線在目標頻點的右側，有k4＜k2＜km＜k1＜k3，并且顯然有0＜k1-km＜1，同理令α= (k1-km- 0.5)∈ (-0.5,+0.5)。類似以上可推導得出頻率、幅值和相位的修正式分別為

利用多項式曲線擬合函數polyfit()，導出本文采用的四項五階Nuttall窗的頻率、幅值修正式為

3 仿真試驗分析

本文所提算法的原理如圖3所示，為驗證本文提出算法的有效性，選取典型的21次諧波信號進行仿真試驗分析，仿真信號模型為

圖3 算法原理Fig.3 The principle diagram of the proposed algorithm

信號參數見表1。依據國際標準IEC61850-9-2LE的規定，本文采用每周波256點的采樣速率，依據國際標準IEC61000-4-7規約的規定，本文采樣10個周波進行分析。整體算法仿真試驗的流程如圖4所示。

圖4 仿真試驗流程Fig.4 Flow chart for simulation experiment

表1 信號參數Tab.1 Signal parameters

3.1 測量準確度分析

取基波頻率f0=50.3Hz，依次采用本文所提算法、補零至4 096點采用基2FFT算法（以下簡稱補零FFT算法）、常規混合基FFT算法對諧波信號進行分析，得到幅值、頻率和相位的測量分析結果如圖5所示。其中補零FFT算法采用4 096點四譜線插值公式，測量結果都以相對誤差的形式給出。圖中，m為諧波次數，E(A)為幅值相對誤差，E(f)為頻率相對誤差，E(φ)為相位相對誤差。

圖5 仿真試驗結果 Fig.5 Emulation results

從圖5的仿真結果可知，常規混合基FFT算法與補零FFT算法的誤差較大，不滿足諧波分析的精度要求。說明補零操作對插值算法產生了影響，補零后不能直接采用補零前的插值算法，同時也說明了常規混合基FFT算法在實際諧波測試環境下的實用性不大。本文算法的準確度較高，仿真結果顯示，本文算法在21次諧波幅值測量的相對誤差不超過7.28×10-4%，相位測量的相對誤差不超過1.13×10-3%，頻率測量的相對誤差不超過1.05×10-5%。故本文提出的混合基FFT多譜線插值算法能有效克服復雜環境的影響，滿足高準確度諧波測量的要求。

3.2 基頻波動影響分析

在實際電網環境中，基波頻率難以避免地存在波動。為驗證本文提出的算法在基頻波動影響下的準確度，調整基波頻率在49.5～50.5Hz區間內逐漸變化，分別采用本文算法、補零FFT算法、常規混合基FFT算法進行仿真試驗分析，得到基波頻率波動影響下基波頻率的測量相對誤差見表2。

表2 基波頻率波動下的基波頻率相對誤差分布Tab.2 Relative errors of fundamental frequency with frequency changing of fundamental

由表2可見，在基波頻率波動影響下，常規混合基FFT算法和補零FFT算法已經失準，本文算法的測量相對誤差整體穩定在10-5%以下。

在不同基波頻率下，采用本文算法對各次諧波幅值進行測量分析，得到各次諧波幅值的相對誤差如圖6所示。

由表2和圖6可以看出，在基波頻率波動影響下，采用本文算法對諧波信號測量分析時，基波頻率相對誤差不超過1.25×10-5%，21次諧波幅值相對誤差保持在10-3%以下，與基波頻率穩定時相比，諧波幅值的測量精度有所下降，但依然遠高于相關標準的測量分析要求。

圖6 基波波動下的諧波幅值計算相對誤差Fig.6 Relative errors of harmonic amplitude with frequency changing of fundamental

3.3 白噪聲影響分析

在實際測量環境下，不可避免地會受到外界干擾噪聲的影響。令白噪聲的信噪比為60dB，分別采用各算法進行測試分析，得到基波、2次諧波和3次諧波的幅值、頻率和相位測量結果見表3。

表3 白噪聲影響下的諧波測量分析結果Tab.3 Analysis results of harmonic measurement under the influence of white noise

由表3可見，在噪聲影響下本文算法對諧波幅值、頻率和相位的測量精度有所降低，但依然高于其余算法，并且滿足相關標準的測量分析要求。

為進一步驗證本文算法在白噪聲影響下的測量準確度，對諧波信號在不同信噪比的白噪聲下進行了仿真分析，仍采用式（18）的諧波信號模型，信噪比SNR范圍30～100dB，信噪比遞進步長為10dB，白噪聲影響下諧波幅值的測試分析結果如圖7所示。

圖7 白噪聲影響下的諧波幅值相對誤差分布Fig.7 Relative error distribution of amplitude under white noise

由圖7可知，當信噪比較低（SNR＜50dB）時，算法的誤差明顯增加，特別是高次諧波（如19次、20次諧波等）的幅值相對誤差較大。當信噪比較高時，21次諧波幅值相對誤差不超過0.02%，各次諧波參數的分析精度幾乎不受影響，遠高于相關標準的要求，故本文算法能有效降低白噪聲對諧波檢測分析的影響。

3.4 運算量分析

由式（1）可見，進行N點DFT需要N2次復乘運算和N(N-1)次復加運算，進行基2 FFT則需要次復乘運算和 2logNN次復加運算。按照國際標準規定，對于采樣10個周波共2 560點數據的情況，采用本文算法的復乘和復加運算量分別為

算法運算量的對比見表4。從表4中可以看出，與DFT算法相比，本文算法的運算量大大減少，與補零至基2FFT算法的運算量相差不多，但精度卻大大增加。故對于非2整次冪的數據序列分析而言，本文算法具有一定的應用價值。

表4 算法運算量對比Tab.4 Comparison of algorithm computation

4 試驗分析

采用德州儀器公司的MSP430F6459作為主控MCU搭建了試驗平臺，其結構如圖8所示。其中，高精密程控源發出的諧波信號經過預處理電路后，被模數轉換器AD7606采樣得到數字信號采樣序列。主控MCU采用本文算法對信號采樣序列分析處理，得到的分析結果將在上位機上顯示。

圖8 試驗平臺框圖Fig.8 Diagram of the test platform

試驗中諧波信號的基波頻率為50.3Hz，在IAR Embedded Workbench編譯軟件中編寫算法后，在此試驗平臺上對多種算法進行了試驗驗證，上位機得到的試驗結果見表5。

表5 試驗結果Tab.5 Experiment results

由表5可以看出，在非2整數次冪數據序列的諧波分析中，常規混合基FFT算法與補零FFT算法的誤差較大，不能滿足諧波分析的精度要求。在此試驗平臺上，本文算法對基波和諧波幅值的分析準確度達到10-3%、10-4%的數量級。若在0.2級的數字化智能電表中進行諧波電能計量，本文算法誤差的影響將會很小。

此外，本文算法在試驗分析中的準確度與仿真分析結果相比有所降低，并非完全一致，其差異主要由信號源自身的誤差、模數轉換器的量化誤差以及電路板上的電磁干擾所引起的。如模數轉換器AD7606的基準溫度系數為±10×10-6℃，若試驗過程中AD7606芯片的溫度變化范圍為±5℃，則引入的誤差為±0.005%，再加上AD7606的積分非線性和微分非線性誤差，模數轉換器引入的誤差約為±0.006 5%。此外，信號預處理電路中的儀表運算放大器在固定增益下也存在精度誤差和溫漂，這些誤差的引入最終都會導致上位機接收到的數據精度有所降低。因此對試驗平臺的進一步完善改進也是后續的研究工作之一。

5 結論

針對諧波分析中采樣數據序列是非2整數次冪的情形，在分析常規混合基FFT算法的不足，并充分考慮目標頻點附近譜線信息的基礎上，提出一種適用于非整數次冪的高精度混合基FFT的諧波測量算法。該算法以2 560點采樣序列為例，將其分解為5點DFT運算和基2FFT運算，依據混合基FFT頻譜分布情況采用四條譜線插值校正計算，并且根據不同情況分別推導出了不同的修正公式。仿真及試驗結果表明，對于非2整數次冪數據序列，相比于其余分析算法，本文算法具有更高的諧波分析準確度，并能有效地改善基波偏移和白噪聲產生的影響，具有很高的應用研究前景。