基于ZoomFFT 和多譜線插值算法的諧波分析方法*

朱 坤 陳樂柱 許 勝 史琛筠

(1.安徽工業大學電氣與信息工程學院 馬鞍山 243000；2.泰州市電能變換與控制工程技術研究中心 泰州 225300)

1 引言

隨著電力電子行業的發展和非線性負荷設備的增加，電網諧波污染進一步惡化，從而影響電力系統安全平穩運行。為了能夠有效治理電網中存在的電能質量問題，首先需要對電網的電力參數進行檢測，為保證對電力系統的有效濾波或補償，需要保證檢測參數的精準性，同時為保證對電力系統的及時治理，還需要保證參數的實時性。

目前，基于快速傅里葉變換(Fast Fourier transform，FFT)的算法因具備計算速度快等優點被廣泛運用，但仍存在頻譜泄漏、柵欄效應和頻譜混疊等問題[1-4]。

理想狀態下時域下的信號序列通常認為是無限長，而實際上在分析時域信號時都是截取某一段信號時序，需要加窗函數，相當于頻域卷積，并且通常情況下在進行時域截取信號時，為非整數信號周期，在分析的頻譜圖上就會存在旁瓣。為了減少頻譜泄漏，需要乘以加權的窗函數，例如Hanning 窗、Hamming 窗、Kaiser 窗、Blackman 窗、Nuttall 窗等。實際在對信號進行處理時，時域和頻域都是離散的，連續的信號不會造成頻率分量丟失，但離散信號時序經過FFT 從時域到頻域會造成頻譜中頻率分量丟失，導致頻譜不完整。

結合工程應用，采集的電信號中會存在與整數次諧波相近頻率的間諧波含量，采用傳統的快速傅里葉變換很難將其區分，此時需要提高頻率分辨率，可以采用提高采樣點數N或者減小采樣頻率的方法。其中也有學者提出頻譜細化算法和加窗插值算法對頻譜進行分析，例如線性調頻 Z 變換(或CZT)[5]、小波變換[6]、復調制頻譜細化(或ZFFT)[7-8]、Kaiser 窗雙譜線插值算法[9]、Kaiser 窗三譜線插值算法[10]和Nuttall 窗三譜線插值[11]等。這些算法都能夠提高頻率分辨率，但仍有其局限性，例如CZT 算法不適用于多頻電信號且相對密集的有效分析，且實時性相對較差；Kaiser 窗和Nuttall 窗的多譜線插值算法測量精度不能滿足所有條件。

針對電網信號測量分析中需要數據精度的同時也需要滿足其實時性的要求，本文結合復調制頻譜細化(ZFFT)算法和基于Blackman-Harris 窗四譜線插值算法對測量信號進行分析，并結合主瓣干擾判定條件對算法進行處理，最后通過仿真驗證改進型方法的精準度和實時性。

2 基于Blackman-Harris 窗四譜線插值校正算法

針對基于Blackman-Harris 窗雙譜線插值算法[12]和基于Blackman-Harris 窗三譜線插值校正算法[13]存在的誤差，本文提出基于Blackman-Harris 窗四譜線插值算法來減少校正誤差，提高測量精確性。

2.1 Blackman-Harris 窗

加窗函數能夠在發生頻譜泄漏時有效改善并降低頻譜泄漏帶來的影響，同時提高頻譜分辨率。通常窗函數具有如下幾種特征參數：主瓣寬度、旁瓣衰減速率、最高旁瓣高度、幅值失真度等。因此，在選擇窗函數時要結合實際測量電信號和分析窗函數的旁瓣特性。

表1 所示為部分窗函數的特性[14]，根據窗函數的特性和電信號的性質來選擇所需的窗函數。

表1 部分窗函數特性

由表1 可知，Blackman-Harris 窗函數的旁瓣峰值電平相比較于其他窗函數較低(-92 dB)，旁瓣衰減速率為6 dB/oct。衰減速率越大，旁瓣峰值電平越小，對頻譜泄漏的治理能力越好。綜合考慮實際情況，本文選取Blackman-Harris 窗函數。

Blackman-Harris 窗的本質是一個4 項系數的加權余弦窗，其函數的時域表達式為

式中，K為余弦窗的項數；N為窗函數的數據長度，n=1,2,3,…,N-1。取K=3，得到一般表達式為

式中，a1=0.358 75，a2=0.488 29，a3=0.141 28，a4=0.011 68。

2.2 四譜線插值算法

為了方便理解和計算，對電信號x(t)采取只含單一次諧波分量的信號進行分析(式(3))，其中以采樣頻率Fs對信號x(t)進行采樣得到的離散時間信號如式(4)所示

式中，f0表示單一次諧波的頻率；A表示單一次諧波的幅值；φ0表示單一次諧波的相位[15]。

現對信號x(n)進行加窗函數w(n)處理，即得到式(5)

代入式(4)，對式(5)進行離散傅里葉變換，得到表達式為

忽略負頻點的旁瓣影響，最終得到簡化表達式為

式中，W[·]為窗函數的離散傅里葉變換形式；Δf為頻率分辨率，Δf=Fs/N，kr=f0/Δf，其中kr為峰值頻點，由于是離散采樣，所以通常kr都是非整數。

對式(2)中Blackman-Harris 窗的時域表達式進行離散傅里葉變換，得到表達式為

式中，b1=0.358 75，b2=0.224 145，b3=0.070 64，b4=0.005 8。WR(·)為矩形窗的頻譜幅度函數表達式。式(8)中W[·]的參數根據現有條件可以轉換為[16]

WR(·)的表達式為

由式(8)可知ω的表達式(10)，將其代入式(12)進行離散采樣

因式(9)中含有WR(·)形式，為利于后面計算所以寫成WR(ω±2πm/N)的形式。

通常情況下對信號采樣很難做到同步采樣，所以f0表示的峰值頻率通常不在離散的頻譜點上，如圖1 所示，在峰值頻點kr兩側分別有兩條最靠近kr頻點的譜線，分別是kr1、kr2、kr3、kr4。分析單譜線、雙譜線和三譜線算法可以知道，靠近峰值頻點kr的譜線理論上擁有更多的頻譜信息，所以要盡可能多地分析譜線信息。

圖1 信號頻譜圖

其中它們之間的關系為：kr1+1=kr2，kr2+1=kr3，kr3+1=kr4。每條譜線幅值分別表示為y1=|XD(kr1Δf)|，y2=|XD(kr2Δf)|，y3=|XD(kr3Δf)|，y4=|XD(kr4Δf)|。引入參數γ=kr-kr2-0.5，根據條件可知γ的取值范圍為[-0.5，0.5]。

引入變量ρ，ρ代表位于峰值頻點kr兩邊的幅值之比，表達式為

將式(8)、(10)一并代入式(13)得到表達式為

從式(14)可以看出，ρ可以表示成γ的函數表達式，ρ=y(γ)，其中反函數可以表示為γ=y-1(ρ)，可以根據ρ求出γ，利用Matlab 中的polyfit 多項式曲線擬合函數擬合出表達式，在調用polyfit 函數時需要選擇擬合多項式的系數n，n取值過大會增加運算成本，通過仿真遍歷可知，當n取值大于7 后，多項式中階數大于7 的項對計算結果影響已經很小，且會增加運算時間，所以通過polyfit 多項式曲線擬合函數擬合的通用表達式為

式中，c1、c3、c5、c2n+1分別代表每個奇數次項(2n+1)前面的系數。

其中信號頻率的校正公式為

由條件可知測量信號的初始相位修正表達式為

測量信號的幅值是根據圖1 中的四條譜線通過加權推算得來，且因為kr2、kr3是最靠近峰值頻點的兩條譜線，所以需要乘以更大的權值，這四條譜線的權值分別取1、3、3、1。由此得到信號幅值修正表達式為

其中，當頻譜分析中N取值較大時，此時可以表達為

g(γ)經過擬合后的表達式可表示為

式中，b0、b2、b2n分別代表偶數次項(2n)前面的系數。

將式(9)和式(12)代入式(14)，利用Matlab 中的曲線擬合函數polyfit()進行求解，為保證數據的精度，在選擇擬合次數時不能過小，此次擬合選取的次數為7，在γ取值范圍內取值代入，通過函數polyfit(ρ，γ，7)擬合求得Blackman-Harris 窗四譜線插值算法中γ為

如圖2 所示，圖2a 為通過擬合得到的曲線圖，圖2b 為真實數值與擬合曲線數值的差值。

圖2 擬合曲線和差值圖

將式(9)、式(12)和式(18)結合，利用曲線擬合函數polyfit(γ，g(γ)，6)求得Blackman-Harris 窗四譜線插值算法中g(γ)為

將求得的式(21)、(22)代入式(16)、(17)、(19)，即可得到基于Blackman-Harris 窗函數四譜線插值算法的頻率、相位和幅值。

3 主瓣干擾判定

3.1 頻率相近分量頻譜分析

一般測量得到的電信號中不僅含有整數次諧波，還含有間諧波，且當間諧波與整數次諧波相近時就會發生頻譜混疊，導致的結果就是不能區分相近的頻率分量。

設待測信號的表達式為

式中，A1、A2表示兩個頻率分量的幅值；f1、f2表示信號中含有的兩個頻率分量；φ1、φ2表示兩個頻率分量的初相位。相近頻率頻譜圖如圖3 所示。

圖3 相近頻率頻譜圖

由圖3 可以看出，當兩個頻率相差比較近時，在頻譜圖中只能看到一個波峰，也就是發生了頻譜混疊，只存在一個主瓣，且發生了主瓣干擾，使得被干擾頻率的譜線幅值也發生了失真現象。如果采用加窗插值算法對測量信號進行分析處理，得到的數據準確度會產生偏差，且不能反映出真實的頻譜情況。

3.2 主瓣干擾分析

現就需要對測量的電信號進行判定，判定的依據為是否發生了主瓣間的干擾。一般可以認為，當沒有發生主瓣干擾時，且不受信號噪聲影響時，主瓣內譜線間的相位和幅值有關聯，主瓣內相鄰譜線間的相位差絕對值等于π，所以在判定主瓣干擾時可以作為依據。

余弦組合窗是目前運用最多的一類窗函數，以此為例，其時域表達式為式(1)，可知其頻譜表達式為

式中，Wo(·)為矩形窗的離散傅里葉變換形式。

將式(25)代入式(24)化簡可得

將式(26)寫成

設加窗信號經過離散傅里葉變換為X(k)，其中在測量第i項諧波時，若不考慮負頻率和周邊諧波的泄漏影響可知

設在主瓣內，k1和k2是相鄰的兩根譜線，k2=k1+1，可知

可知兩條譜線的相位差為

由此可知在不受信號噪聲及主瓣干擾的情況下，主瓣內相鄰譜線間的相位差為π。

引入參數Δθ，根據表達式

式中，δ1、δ2分別是靠近峰值頻點最近的兩個譜線相位。按照實際工程應用來說，Δθ通常都不等于0，所以需要選定一個恰當的值￡來權衡是否發生了主瓣干擾，￡值的大小直接決定了測量時長和精準度。如果￡的取值過大，就會導致發生了主瓣干擾但是沒有判別出來的情況，從而導致漏判；若￡的取值過小，將會導致過于靈敏，使得計算時長增加[17]。本文結合實際工程應用，對不同頻率差信號分量進行仿真遍歷，在引入較少運算量的同時對諧波分量進行有效分離檢測，最終選取￡=8×10-6。

4 復調制頻譜細化算法

當測量電信號的某一頻段上判定了主瓣干擾，說明在這一頻段上存在不止一個頻率分量，且它們的頻率值相差很小，已經發生了頻譜混疊的現象，這時就需要利用頻率分辨率更高的處理算法來區分在這一頻段上存在的頻率分量[18]。本文選取的是復調制頻譜細化算法ZoomFFT(或ZFFT)。

ZFFT 算法能夠實現在某一頻段上對頻譜進行放大處理，使其在這一頻段上頻率分辨率提高，放大D倍[19]。對于頻率分辨率Δf=Fs/N，根據公式可知想要提高頻率分辨率，可以改變Fs和N這兩個值。

(1) 保持N值不變，減小采樣頻率Fs的值，使得Δf減小，也就是頻率分辨率提高，但是根據奈奎斯特采樣定理可知，如果采樣頻率Fs小于兩倍信號中的最高頻率，測量的信號頻譜將會發生頻譜混疊現象，即要保證Fs≥2fmax。

(2) 保持采樣頻率Fs不變，增大采樣點數N，也能提高頻率分辨率，但是隨著N值增大，帶來的是計算量的增大，從而導致計算時長的增加，這對一些實時性要求高的場合就不太適用，且會占用大量內存。

ZoomFFT 算法實現的流程圖如圖4 所示。

圖4 ZFFT 算法流程圖

ZFFT 算法實現的過程總共分為7 個部分，其具體步驟如下所示。

(1) 濾波。根據奈奎斯特采樣定理可知，FFT的頻率分析范圍有限，2/Fs為最高的分析頻率，所以需要濾除其頻率以上的頻率，防止在分析頻譜時出現混疊。

(2) 離散采樣。對時域信號x(t)進行離散化處理，得到x(n)，采樣長度為放大倍數D與采樣點數的乘積。

(3) 頻移。經過處理后的離散信號x(n)乘以旋轉因子exp(-j2πfn)，目的是為了讓所需分析的頻段中心頻率f移至頻率零點，得到新的信號。

(4) 低通濾波。目的是為了濾除所需頻譜細化頻段以外的信號分量，只保留所需細化的頻段。

(5) 重采樣。對經過低通濾波器之后保留的頻段進行重新采樣，采樣時每個點的間隔為D，采樣頻率由之前的Fs變為Fs/D，由此可見，采樣頻率降低了D倍。

(6) FFT 分析。對重采樣后的離散信號進行FFT分析，分析的點數仍然為N，但是采樣頻率由之前的Fs變為Fs/D，所以頻率分辨率也會提升D倍。

(7) 頻率調整。由步驟(3)可知信號經過頻移，所以需要對最終得到的頻譜信號進行調整。

綜上，通過ZFFT 算法處理過后的頻譜在以頻率f為中心的頻段上，將頻率分辨率由之前的Δf=Fs/N降為Δf=Fs/ND，也就是將頻率分辨率提升了D倍。放大的頻段范圍是ΔF=f±Fs/2D。

觀察ZFFT 算法流程步驟可知，ZFFT 算法只能針對想要觀察的目標頻段進行放大處理，不能夠對全頻域內進行放大。結合如今的電力系統測量分析可知，大多分析為50 次諧波以內，所以一般分析0～2.5 kHz 諧波含量，若只對存在主瓣干擾的頻段采用頻譜細化分析，對只含有整數次諧波的頻段采用加窗插值算法分析，這將會大大減少計算時間，不僅提高了諧波檢測精度，同時也提高了實時性。放大倍數的選擇也決定了ZFFT 算法的有效性，如果放大倍數過小，則不能夠對存在主瓣干擾的頻率實現有效區分，若放大倍數過大，則會導致采樣的長度過長。

綜上所述，對測量電能質量提出了改進方法，流程如圖5 所示。

圖5 改進型測量方法流程圖

根據圖5 可知，結合Blackman-Harris 窗插值算法和ZFFT 算法，利用條件主瓣干擾判定選用算法，彌補了ZFFT 算法不能全頻域對頻譜進行細化，基于Blackman-Harris 窗四譜線插值不能夠精確區分頻率相近分量的缺點。

5 仿真驗證與分析

為了驗證所提改進型方法的準確度和實時性，這里對整數次諧波和間諧波進行了頻譜分析。主要的電力參數設置如表2 所示。

表2 電力參數

結合實際電力系統工程分析，此次試驗分析的諧波范圍為0～50 次，相應頻譜范圍為0～2.5 kHz，最高的頻率為2.5 kHz。因此，根據奈奎斯特采樣定理選取采樣頻率為5 kHz，采樣點數N=1 024，ZoomFFT 選取的放大倍數D=10。

此外，在表2 中，參數中除了含有基波和其他整數次諧波含量，還設置了間諧波且頻率相近。

圖6 為改進型方法仿真局部頻譜圖。根據仿真結果實測數據得到參數結果如表3 所示。

圖6 改進型方法仿真頻譜圖

表3 改進型方法仿真結果數據

同時分別對傳統的 FFT 諧波檢測方法和Blackman-Harris 窗插值算法進行仿真驗證對比，圖7為基于Blackman-Harris 窗四譜線插值算法仿真頻譜圖，圖8 為FFT 算法仿真頻譜圖。信號仿真設置參數如表2 所示。

圖7 基于Blackman-Harris 窗四譜線插值算法仿真頻譜圖

圖8 傳統FFT 算法仿真頻譜圖

對比圖6、圖7、圖8 和表3、表4 可知，信號分量X2在基于Blackman-Harris 窗四譜線插值算法中無法被檢測出，不能將所有的間諧波區分開來，同樣條件下，采用本文所提方法能夠將所有間諧波全部區分并取得精確的參數值。

表4 基于Blackman-Harris 窗四譜線插值算法仿真結果數據

通過對比不同諧波檢測分析的方法可知，對電信號加不同的窗函數和采用不同的插值算法得到的時間成本都不相同[20]，其中ZoomFFT 算法在提高頻率分辨時是以分析更多的采樣點數得到的。如果只使用ZFFT 算法對信號進行分析，將會耗費大量的時間和內存單元。通過仿真對比，不同加窗插值算法的平均計算時長如表5 所示。由表5 可知，ZFFT算法精確度高，但計算時間過長，不利于測量要求的實時性，且會占據大量內存。因此，本文采用的方法兼容了ZFFT 和加Blackman-Harris 窗插值算法的優點。

表5 仿真平均計算時長

6 結論

本文結合實際的工程應用，對電網電能質量參數進行快速及準確的采集和處理。針對電網諧波分量發生的頻譜泄漏和混疊現象，造成難以精確區分相近諧波分離檢測的問題，本文采用基于ZoomFFT 和多譜線插值算法的諧波分析方法，得出如下結論。

(1) 基于Blackman-Harris 窗四譜線插值算法，能夠滿足諧波高精度測量且計算時長少的要求。

(2) 在發生主瓣干擾的情況下，采用頻率分辨率更高的頻譜細化ZoomFFT 算法精準分析。

(3) 通過仿真數據對比分析驗證了本文方法的精確性和實時性，結合兩種算法各自優點，保證了測量精度，同時也滿足了更少的運算時長，從而確保了現代電力系統發展所需的要求。