混合自卷積窗諧波算法研究與FPGA方案

賴增強，王貴忠，肖智宏，張國慶，李洪波，劉穎，關(guān)鐵成

(1.哈爾濱工業(yè)大學(xué) 電氣工程及自動化學(xué)院，黑龍江 哈爾濱 150090；2.哈工大(張家口) 工業(yè)技術(shù)研究院，河北 張家口 075421；3.國網(wǎng)經(jīng)濟技術(shù)研究院有限公司，北京 102209；4.國網(wǎng)營口供電公司，遼寧 營口 115002)

隨著電力系統(tǒng)電力電子化進程的推進，大量電力電子裝置的并網(wǎng)使用使得電網(wǎng)中的諧波檢測與治理問題得到人們的廣泛關(guān)注，諧波分析在電力系統(tǒng)中顯得尤為重要。電力諧波分析算法一般分為如下4種[1-2]：①時域仿真分析法；②頻域仿真分析法；③基于數(shù)學(xué)變換理論算法；④基于人工智能算法。迄今為止文獻數(shù)量最多、研究最深入、運用最廣泛當屬基于數(shù)學(xué)變換理論算法，它主要包括快速傅里葉變換(fast Fourier transform, FFT)、小波變換分析法、矢量分析法、Prony法、廣義S變換法以及相關(guān)改進算法。其中，F(xiàn)FT是把時域與頻域相關(guān)聯(lián)的算法，具有完備、正交等眾多優(yōu)勢。盡管FFT不如小波變換那樣特別適合處理突變和不平穩(wěn)信號的分析，但它的改進算法因其計算方便特別適合現(xiàn)場可編程門陣列(field-programmable gate array, FPGA)等硬件編程實現(xiàn)，是電能質(zhì)量裝置中諧波分析的首選算法。然而FFT不足之處是非同步采樣和非整數(shù)周期截斷時易出現(xiàn)“頻譜泄漏”與“柵欄效應(yīng)”，致使不能精確地測量出電網(wǎng)的諧波數(shù)據(jù)[3]；同時對于電網(wǎng)的頻率波動，即使是鎖相環(huán)也難以實現(xiàn)電網(wǎng)的同步采樣[4]，因此需要采用一種性能優(yōu)良的窗函數(shù)對FFT算法進行加窗處理來抑制頻譜效應(yīng)。常用的傳統(tǒng)窗函數(shù)有矩形窗、Hanning窗、Hamming窗和Blackman窗等[5-9]。

電力諧波測量標準要求幅值誤差不應(yīng)大于5%，相角誤差不大于5°[10]。為了提高諧波分析精確度，文獻[11-12]提出了Nuttall窗及其自卷積窗算法，文獻[13]提出了Hanning自卷積窗算法，自卷積窗有利于提高窗的本身優(yōu)勢，但是同時也放大了窗的一些劣勢，如Nuttall窗自卷積后主瓣寬度增加。文獻[14]提出了不同窗混合卷積方法，即將兩種窗的優(yōu)勢互補，進而提高計算精度。文獻[15-16]提出了五點加權(quán)法，該方法雖精度高但因?qū)FT輸出序列進行加權(quán)變換導(dǎo)致計算量大，難以硬件實現(xiàn)。為了得到頻譜特性更為優(yōu)越的窗函數(shù)，本文提出一種新型Nuttall混合自卷積窗函數(shù)，即先將Nuttall窗與Hanning窗混合卷積后再L階自卷積，并對其進行雙譜線插值與數(shù)值擬合修正。該算法利用Nuttall窗旁瓣峰值衰減速度快和Hanning窗主瓣較窄的特點，經(jīng)自卷積進一步提升優(yōu)勢，從而得到的混合自卷積窗十分適用于諧波分析。最后通過仿真對比，驗證混合自卷積窗諧波算法的準確性。同時對算法進行了硬件結(jié)構(gòu)方案設(shè)計與FPGA關(guān)鍵技術(shù)研究，為諧波分析儀和電能質(zhì)量裝置的研制提供參考。

1 混合自卷積窗設(shè)計

1.1 Hanning窗與Nuttall窗

Hanning窗為一種升余弦窗，其離散表達式為

(1)

式中：n為采樣樣點序列號;N為采樣樣點總數(shù);RN(n)為矩形窗函數(shù)。

式(1)進行離散傅里葉變換(discrete Fourier transform, DFT)后的頻域表達式為

(2)

式中WR(·)為矩形窗的DFT表達式。

Nuttall窗也為一種余弦窗，函數(shù)表達式為

(3)

式中al為系數(shù)。

式(3)進行DFT后的頻域表達式為

(4)

由于Nuttall具有多項多階，即不同的L值具有不同的al系數(shù)。Nuttall窗的多項多階與系數(shù)對應(yīng)關(guān)系見表1。

表1 Nuttall窗的多項多階與系數(shù)對應(yīng)關(guān)系Tab.1 Corresponding relation between multiple orders and coefficients of Nuttall window

本文將采用4項3階Nuttall窗，即

(5)

給出Hanning窗與Nuttall窗的幅頻響應(yīng)曲線，且幅頻圖中以歸一化后的頻率做橫坐標，以最大值為0分貝值做縱坐標，如圖1所示。

圖1中，不僅表明了窗函數(shù)的幅頻關(guān)系，更包含了一些十分重要的指標，如主瓣寬度、第一旁瓣衰減以及旁瓣峰值衰減等。根據(jù)圖1，不難發(fā)現(xiàn)Nuttall窗旁瓣峰值衰減速度快及Hanning窗主瓣比較窄。實際上，窗函數(shù)第一旁瓣衰減大和旁瓣峰值衰減速度快最能有效抑制頻譜效應(yīng)，但是主瓣寬度較大往往又給信號分析帶來不良影響，所以這是一個綜合評估并折中選擇的過程[17-18]。

-32 dB—Hanning窗第一旁瓣衰減；-83 dB—Nuttall窗第一旁瓣衰減；-18 dB/oct—Hanning窗旁瓣峰值衰減速度；-30 dB/oct—Nuttall窗旁瓣峰值衰減速度，dB/oct是衰減斜率的單位分頻斜率，用來反映分頻點以下頻響曲線的下降斜率,用分貝/倍頻程表示，即dB/oct。

圖1 Hanning窗與 Nuttall窗的幅頻響應(yīng)對比
Fig.1 Contrast of amplitude frequency response ofHanning and Nuttall window

1.2 混合自卷積窗函數(shù)設(shè)計

Nuttall混合自卷積窗函數(shù)，是先將Nuttall窗與Hanning窗混合卷積，得到混合窗函數(shù)

wNHC(n)=wN(n)*wH(n).

(6)

再將得到的wNHC(n)進行L個時域卷積計算，即L階自卷積，推導(dǎo)出混合自卷積窗函數(shù)

(7)

圖2給出了不同L階混合自卷積窗函數(shù)的幅頻響應(yīng)曲線。從圖2可以看出，隨著L的增加，混合自卷積窗函數(shù)的主瓣寬度緩慢增寬，但旁瓣特性卻迅速提升，所以對于電力諧波信號，選擇合適的L階混合自卷積窗函數(shù)對信號加窗，可有效提升諧波的FFT分析精度。結(jié)合窗函數(shù)選擇原則，若無法同時滿足主瓣寬度與旁瓣特性，常常以犧牲主瓣的寬度來重視旁瓣特性的要求[19-20]，故以下將選擇L=4作為混合自卷積窗函數(shù)的階數(shù)來做進一步的研究。

2 混合自卷積窗雙譜線插值算法研究

電力信號實際上并非同步采樣，而且很可能是非整周波采樣，因此在離散頻率點上不一定存在峰值頻率為基波整數(shù)倍的諧波,但雙譜線插值算法能很好地解決這一問題[11,21]。

圖2 不同L階混合自卷積窗函數(shù)的幅頻響應(yīng)對比Fig.2 Contrast of amplitude frequency ofNHSC window function with different L orders

對時域信號x(t)進行離散采樣，則

(8)

式中：f0為頻域離散信號x(n)的初始頻率；A為x(t)幅值；φ0為x(t)初相角;fs為采樣頻率;t為時間變量。

對信號x(n)加窗及快速傅里葉變換后經(jīng)離散時間傅里葉變換(discrete-time Fourier transform，DTFT)離散抽樣，若抽樣時間間隔為Δf=fs/N,則

(9)

式中：X(kΔf) 為第k條譜線的DTFT的表達式;W(·)為窗函數(shù)的DTFT表達式。

設(shè)k1、k2分別為峰值點k0左右兩側(cè)幅值為最大和次大的譜線，可分別表示為y1=|X(k1Δf)|和y2=|X(k2Δf)|，且滿足k1≤k0≤k2=k1+1。同時引入?yún)?shù)α，并令α=k0-k1-0.5，則α的取值范圍為[-0.5,0.5]。定義β=(y2-y1)/(y2+y1)，則

(10)

由式(10)定義α和β的函數(shù)關(guān)系，記做β=g(α)或α=g-1(β)，并得到了頻率修正式f0=k0Δf=(α+k1+0.5)Δf，同時定義兩譜線加權(quán)平均值為幅值修正式，即

(11)

式中A1、A2分別為第k1和k2條譜線所對應(yīng)的幅值。

當數(shù)據(jù)階段長度N很大時，式(11)還可做進一步簡化，即

A=N-1(y1+y2)v(α).

(12)

式中v(α)為下文式(14)的偶函數(shù)形式。同時，由式(9)，可得初相位修正式為

(13)

綜上所述，雙譜線插值算法計算繁瑣，特別是應(yīng)用在嵌入式實時監(jiān)測系統(tǒng)中，加之混合自卷積窗對系統(tǒng)的運算量提出了很高的要求。為此，采用一種多項式逼近的數(shù)值分析算法求取α=g-1(β)。設(shè)α系列α1,α3,…,α2m+1取值范圍為[-0.5,0.5]，代入式(10)，得到與之對應(yīng)的β值，通過數(shù)據(jù)擬合從而得到

(14)

同時初相角的修正式也可化簡為

(15)

針對混合自卷積窗，為了保證精度的同時提升其計算速度，本文將采用雙譜線插值五次多項式逼近的修正式，同時給出修正式α、β的擬合曲線，如圖3所示；不僅如此，為了對比分析，也分別給出了性能較優(yōu)的Hanning窗和Nuttall窗的修正式，即：

a)混合自卷積窗函數(shù)及修正式

(16)

b)Hanning窗及修正式

(17)

c)Nuttall窗及修正式

(18)

圖3 混合自卷積窗的修正式擬合曲線Fig.3 Correction formula fitting curve of NHSC window

3 仿真驗證

為了驗證理論分析，采用MATLAB編程仿真，對表2中一組諧波數(shù)據(jù)信號分別用混合自卷積窗、Nuttall窗和Hanning窗進行處理。通過觀察各類窗雙譜線插值后的數(shù)據(jù)結(jié)果，對比分析各類窗的優(yōu)劣。

設(shè)信號基波頻率為49.5 Hz，采樣頻率為3 200 Hz。

MATLAB仿真流程如圖4所示，其中輸入信號為x(n)，分別為混合自卷積窗、Nuttall窗、Hanning窗,同時令基波頻率、幅值及相位角分別為f0,A0,φ0，其他各次諧波頻率、幅值及相位角分別為fi,Ai,φi。

表2 諧波數(shù)據(jù)表Tab.2 Harmonic data table

圖4 仿真程序流程Fig.4 Simulation program flow chart

通過程序仿真，得到的具體數(shù)據(jù)見表3—表5，其中相位及幅值誤差如圖5和圖6所示。

表3 混合自卷積窗諧波數(shù)據(jù)Tab.3 Harmonic data of NHSC window

表4 Hanning窗諧波數(shù)據(jù)Tab.4 Harmonic data of Hanning window

表5 Nuttall窗諧波數(shù)據(jù)Tab.5 Harmonic data of Nuttall window

圖5 3類窗諧波分析相位誤差對比Fig.5 Phase errors of harmonic analysis of three kinds of window functions

從表3—表5數(shù)據(jù)及圖5和圖6誤差中可以分析得出：無論在相位還是幅值檢測上，較之性能較好的Nuttall窗與Hanning窗，混合自卷積窗更能體現(xiàn)出誤差小、精度高等強大優(yōu)勢。其中，相位誤差穩(wěn)定在0.000 01°，幅值誤差幾乎為0，完全滿足諧波測量標準要求，該算法對弱幅值諧波也有較強的檢測精度。

圖6 3類窗諧波分析幅值誤差對比Fig.6 Amplitude errors of harmonic analysis of three kinds of window functions

4 硬件結(jié)構(gòu)與FPGA設(shè)計方案

FPGA作為集微處理器、存儲設(shè)備和輸入輸出接口等模塊于一身的系統(tǒng)級集成電路，因其功能強大、保密性能好以及良好的可重復(fù)編程等特性，深受工業(yè)界的青睞。因此本算法的硬件設(shè)計方案采用FPGA作為主控芯片。本算法的硬件全局結(jié)構(gòu)如圖7所示。

圖7中，標準100 V/5 V電壓互感器和5 A/2.5 mA電流互感器分別把三相電壓信號和三相電流信號轉(zhuǎn)換為6路小電壓及小電流信號，并通過信號調(diào)制電路得到適合數(shù)模轉(zhuǎn)換的量程；其次FPGA控制數(shù)模轉(zhuǎn)換后的數(shù)據(jù)信號緩存到輸入雙口數(shù)據(jù)存儲器中，然后運用數(shù)據(jù)處理模塊對數(shù)據(jù)進行類型轉(zhuǎn)換、混合自卷積窗處理、FFT運算以及有效值處理等，并且整個處理過程是由有限狀態(tài)機控制，并緩存到輸出雙口數(shù)據(jù)存儲器中，作進一步處理運算；最后進行存儲，并利用串口將處理好的數(shù)據(jù)輸出到上位機液晶屏中顯示，從而實現(xiàn)電能質(zhì)量諧波分析與監(jiān)測的功能。

基于混合自卷積窗函數(shù)雙譜線插值算法的硬件移植，包含混合自卷積窗和FFT運算兩大部分。主要思路是運用摩爾狀態(tài)機的控制將256個點的采樣數(shù)據(jù)緩存到雙口數(shù)據(jù)存儲器中進行加窗、雙峰譜線插值修正和FFT運算，其具體FFT程序控制流程如圖8所示。

對于混合自卷積窗及插值修正，可以事先利用MATLAB獲得修正后的混合自卷積窗的.mif數(shù)據(jù)文件，并將其存儲到硬件電路的芯片中，可利用Quartus II 13.0設(shè)置為數(shù)據(jù)存儲，具體程序架構(gòu)如圖9所示。

在實際Verilog編程中，流水線乘法器均被混合自卷積窗與采集數(shù)據(jù)所采用，因此可采用乘法單元的知識產(chǎn)權(quán)核來設(shè)計，便能將此算法在FPGA上實現(xiàn)電力諧波分析。

圖7 硬件裝置結(jié)構(gòu)Fig.7 Hardware device structure

圖8 FFT程序流程Fig.8 FFT program flow chart

圖9 混合自卷積加窗插值程序架構(gòu)Fig.9 NHSC Window interpolation program architecture

5 結(jié)束語

本文提出了一種新型的混合自卷積窗，分析了Hanning窗與Nuttall窗各自的優(yōu)勢，且利用二者混合卷積后再L階自卷積的方式，使得新型窗—混合自卷積窗旁瓣特性大幅度提升，且性能不斷優(yōu)化。在此基礎(chǔ)上，因考慮實際檢測諧波是非周期且非同步采樣，故對混合自卷積窗進行了雙譜線插值，進一步提高了FFT算法的準確度，且抑制了一些頻譜效應(yīng)。對所采用的算法進行對比仿真驗證，結(jié)果表明本文提出的基于混合自卷積窗的FFT插值修正算法誤差小、精度高，適合電力系統(tǒng)諧波分析。最后給出了FPGA的硬件設(shè)計方案，運用MATLAB獲取一系列窗的數(shù)值表，與采集后離散化數(shù)據(jù)做乘法加窗處理，在保證算法計算速度的同時，比起傳統(tǒng)FFT算法，提高了測量精度，為相關(guān)諧波分析儀與電能質(zhì)量裝置的研制提供了一定的理論與技術(shù)支撐。