賴(lài)增強(qiáng)，王貴忠，肖智宏，張國(guó)慶，李洪波，劉穎，關(guān)鐵成

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

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

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

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

1.1 Hanning窗與Nuttall窗

Hanning窗為一種升余弦窗，其離散表達(dá)式為

(1)

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

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

(2)

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

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

(3)

式中al為系數(shù)。

式(3)進(jìn)行DFT后的頻域表達(dá)式為

(4)

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

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

本文將采用4項(xiàng)3階Nuttall窗，即

(5)

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

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

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

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

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

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

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

(6)

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

(7)

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

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

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

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

對(duì)時(shí)域信號(hào)x(t)進(jìn)行離散采樣，則

(8)

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

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

(9)

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

設(shè)k1、k2分別為峰值點(diǎn)k0左右兩側(cè)幅值為最大和次大的譜線，可分別表示為y1=|X(k1Δf)|和y2=|X(k2Δf)|，且滿足k1≤k0≤k2=k1+1。同時(shí)引入?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，同時(shí)定義兩譜線加權(quán)平均值為幅值修正式，即

(11)

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

當(dāng)數(shù)據(jù)階段長(zhǎng)度N很大時(shí)，式(11)還可做進(jìn)一步簡(jiǎn)化，即

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

(12)

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

(13)

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

(14)

同時(shí)初相角的修正式也可化簡(jiǎn)為

(15)

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

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

(16)

b)Hanning窗及修正式

(17)

c)Nuttall窗及修正式

(18)

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

3 仿真驗(yàn)證

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

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

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

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

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

通過(guò)程序仿真，得到的具體數(shù)據(jù)見(jiàn)表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類(lèi)窗諧波分析相位誤差對(duì)比Fig.5 Phase errors of harmonic analysis of three kinds of window functions

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

圖6 3類(lèi)窗諧波分析幅值誤差對(duì)比Fig.6 Amplitude errors of harmonic analysis of three kinds of window functions

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

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

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

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

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

在實(shí)際Verilog編程中，流水線乘法器均被混合自卷積窗與采集數(shù)據(jù)所采用，因此可采用乘法單元的知識(shí)產(chǎn)權(quán)核來(lái)設(shè)計(jì)，便能將此算法在FPGA上實(shí)現(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é)束語(yǔ)

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