利用加窗FFT和IFFT實(shí)現(xiàn)Hilbert變換*

魏 國(guó),張 蓓,撒文彬,孫金瑋

(哈爾濱工業(yè)大學(xué)電氣工程學(xué)院,哈爾濱 150001)

電力系統(tǒng)無(wú)功功率是電力系統(tǒng)中一個(gè)非常重要的技術(shù)指標(biāo)。控制無(wú)功的基礎(chǔ)是無(wú)功功率測(cè)量的實(shí)時(shí)性、有效性和準(zhǔn)確性。基于C.Budeanu對(duì)無(wú)功功率的定義[1],目前主要的無(wú)功功率測(cè)量方法有移相法[2]、平均功率瞬時(shí)檢測(cè)法[3]、基于快速傅里葉變換FFT(fast Fourier transform)的方法[4],還有基于H ilbert變換[5,6]和基于H ilbert數(shù)字濾波器[7～10]的方法等。其中Hilbert變換法較其他方法精度更高,能更精確地滿足無(wú)功功率的定義,是應(yīng)用較廣泛的一種方法。

傳統(tǒng)離散時(shí)間Hilbert變換是基于離散傅里葉變換DFT(discrete Fourier transform)和離散傅里葉反變換IDFT(inverse discrete Fourier transform)實(shí)現(xiàn)的。根據(jù)這種方法,在實(shí)際應(yīng)用中只需離線求出一個(gè)Hilbert變換矩陣H,然后用采樣序列與H相乘即可得到移相90°后的信號(hào)。經(jīng)實(shí)際驗(yàn)證,該算法具有較高的精度,在信號(hào)的采樣點(diǎn)數(shù)少的情況下,測(cè)量較為方便,計(jì)算時(shí)間也能基本滿足實(shí)際工程需求。但進(jìn)一步測(cè)試發(fā)現(xiàn),若采樣點(diǎn)數(shù)增多,H矩陣的維數(shù)迅速增大,信號(hào)的Hilbert變換需要消耗大量時(shí)間,無(wú)功測(cè)量的實(shí)時(shí)性和有效性已經(jīng)不能滿足實(shí)際要求。經(jīng)分析,若采用加窗FFT和IFFT方法實(shí)現(xiàn)信號(hào)Hilbert變換,可有效縮短計(jì)算時(shí)間,提高無(wú)功測(cè)量的實(shí)時(shí)性,其變換效果也優(yōu)于傳統(tǒng)方法。

1 Hilber t變換原理

理想的H ilbert變換是幅頻特性為1、負(fù)頻率成分進(jìn)行90°相移、正頻成分進(jìn)行-90°相移的線性變換,其頻域傳遞函數(shù)為

設(shè)有實(shí)信號(hào)f(t),其傅氏變換為F(ω)。令R1(ω)、I1(ω)分別表示ω＞0時(shí)的頻譜實(shí)部與虛部,有:

令f′(t)=f(t)+j(t),則f′(t)的傅里葉變換為

由式(4)可知,f′(t)的正頻分量為f(t)的兩倍,負(fù)頻分量為0。實(shí)際應(yīng)用中,可先求出信號(hào)f(t)的傅氏變換F(ω),然后根據(jù)式(4)求得F′(ω),進(jìn)而得到F′(ω)對(duì)應(yīng)的時(shí)域信號(hào)f′(t)。f′(t)的虛部即為經(jīng)過(guò)H ilbert變換后的信號(hào)(t)。

2 傳統(tǒng)的離散時(shí)間Hilbert變換方法

以采樣頻率f s對(duì)電壓信號(hào)采樣m個(gè)周期,每周期采樣M點(diǎn),得N點(diǎn)序列u(n)。令u(n)基波頻率真值為f1,所含最高諧波次數(shù)為h,傳統(tǒng)的離散時(shí)間H ilbert變換方法如下。

1)對(duì)u(n)進(jìn)行DFT運(yùn)算,得到信號(hào)各次諧波分量。

式中U(0),U(m),…,U(mh)等效為序列u(n)的離散時(shí)間傅里葉變換DTFT(discrete time Fourier transform)在[0,f s]區(qū)間上以m×為間隔的h+1點(diǎn)采樣,分別表示信號(hào)的直流分量、基波、…、h次諧波分量。

2)對(duì)各次諧波分量加權(quán),然后進(jìn)行IDFT運(yùn)算。

3)綜合前兩步,序列u(n)的Hilbert變換過(guò)程等效于H矩陣與u(n)乘積的過(guò)程,即u^(n)=H u(n),其中H為N×N階實(shí)矩陣。

傳統(tǒng)方法簡(jiǎn)單直觀,易于軟件實(shí)現(xiàn),在采樣點(diǎn)數(shù)較少時(shí)不失為一種好方法。但信號(hào)采樣點(diǎn)數(shù)增多后,H矩陣維數(shù)變大導(dǎo)致計(jì)算量快速增大,降低了無(wú)功測(cè)量的實(shí)時(shí)性。下文對(duì)傳統(tǒng)方法進(jìn)行了改進(jìn)。

3 加窗FFT和IFFT實(shí)現(xiàn)Hilbert變換

3.1 用FFT替代傳統(tǒng)方法中DFT運(yùn)算

對(duì)于長(zhǎng)度為N的序列x(n),其N點(diǎn)DFT為

令序列長(zhǎng)度為N=2M,對(duì)x(n)定義另一種形式的DFT運(yùn)算:

則有

考慮到

且N=2M,式(10)變?yōu)?

若x(n)長(zhǎng)度為N=m×M,參照上述推導(dǎo)有

就等同于

容易看出,Xs(k)序列等效為x(n)的DTFT在區(qū)間[0上以m×為間隔的M點(diǎn)采樣。

參照上述推導(dǎo),對(duì)于式(5)中的電壓序列u(n),若將其m個(gè)周期上的采樣點(diǎn)對(duì)應(yīng)累加到一個(gè)周期并做M點(diǎn)FFT,所得M點(diǎn)序列Us(k)即等效為u(n)的DTFT在[0,f s]上以f為間隔的M點(diǎn)采樣。其中,顯然,U s(k)與式(5)中U(mk)的意義完全相同。由此,式(5)可變?yōu)?

求取U(0)～U(mh)的過(guò)程簡(jiǎn)化為:(1)將N點(diǎn)序列u(n)累加到一個(gè)周期,得M點(diǎn)序列us(n);(2)對(duì)u s(n)做M點(diǎn)FFT,得U s(k),k=0,1,…,M-1。Us(0),Us(1),…,Us(h)即為U(0),U(m),…,U(mh),分別表示信號(hào)的直流分量、基波、…、h次諧波分量。

3.2 用IFFT替代傳統(tǒng)方法中IDFT運(yùn)算

在式(6)中,令:

對(duì)B做如下處理。

1)在B的每?jī)蓚€(gè)元素間補(bǔ)m-1個(gè)零,在最后一個(gè)元素后補(bǔ)N-1-mh個(gè)零,得:

2)在B后補(bǔ)上M-h-1個(gè)零,得:

寫(xiě)出N點(diǎn)IFFT的標(biāo)準(zhǔn)系數(shù)矩陣A1。

易得:

可見(jiàn),式(6)的運(yùn)算過(guò)程等效于對(duì)B1的N點(diǎn)IFFT所得時(shí)域序列b1(n)取虛部的過(guò)程,即

由于B1中部分元素為0,b1(n)可寫(xiě)為

易得:

參照式(22)、(24)和(26),式(6)的運(yùn)算最終等效為:1)對(duì)B2做M點(diǎn)IFFT運(yùn)算并將所得M點(diǎn)時(shí)域序列除以m,得序列b1(n)(n=0,1,…,M-1);2)根據(jù)采樣周期數(shù)m延拓b1(n)(n=0,1,…,M-1),得b1(n)(n=0,1,…,N-1);3)對(duì)N點(diǎn)長(zhǎng)度的b1(n)序列取其虛部即得到電壓移相信號(hào)(n)。

3.3 窗函數(shù)應(yīng)用于Hilbert變換

傳統(tǒng)的H ilbert變換方法需先求出信號(hào)各次諧波參數(shù)。由于電網(wǎng)信號(hào)的頻率通常會(huì)在額定頻率附近波動(dòng),同步采樣不可能實(shí)現(xiàn)。對(duì)應(yīng)式(14)有:

采樣的非同步性導(dǎo)致信號(hào)FFT分析時(shí)各諧波間出現(xiàn)頻譜泄漏。同時(shí)FFT算法本身又會(huì)引入柵欄效應(yīng),這兩因素都將降低諧波參數(shù)的測(cè)量精度。對(duì)應(yīng)式(15)中,U(0)～U(mh)的值不再準(zhǔn)確。

若選擇合適的窗函數(shù)則可有效地抑制頻譜泄漏,同時(shí),柵欄效應(yīng)可以通過(guò)基于窗函數(shù)的插值算法加以克服。但插值算法要求精確求解信號(hào)采樣同步誤差,其計(jì)算量通常較大,消耗時(shí)間多,影響無(wú)功功率測(cè)量的實(shí)時(shí)性。實(shí)際應(yīng)用中,只需根據(jù)信號(hào)的采樣周期數(shù),選擇合適的窗函數(shù)進(jìn)行運(yùn)算,所得諧波參數(shù)測(cè)量結(jié)果就足以滿足精度要求。目前比較常用的窗函數(shù)為組合余弦窗函數(shù),如Hanning窗、Blackman窗、Blackman-H arris窗、矩形自卷積窗[6]等。

設(shè)長(zhǎng)度為N的窗序列為:

在H ilbert變換中,對(duì)信號(hào)u(n)加窗得序列:

綜上所述,只需對(duì)長(zhǎng)度為m×M的周期離散序列加窗,然后進(jìn)行1次M點(diǎn)的FFT和1次M點(diǎn)的IFFT運(yùn)算即可準(zhǔn)確實(shí)現(xiàn)其H ilbert變換。與傳統(tǒng)的離散Hilbert變換方法相比,改進(jìn)方法有效簡(jiǎn)化了計(jì)算,并且對(duì)信號(hào)進(jìn)行了預(yù)加窗處理,改善了變換的效果。稱上述算法為基于加窗FFT和IFFT的H ilbert變換方法。

4 計(jì)算量比較及驗(yàn)證

令1次實(shí)數(shù)乘法所需運(yùn)算時(shí)間為a,1次實(shí)數(shù)加法運(yùn)算時(shí)間為b,1次復(fù)數(shù)乘法運(yùn)算時(shí)間為a1,1次復(fù)數(shù)加法運(yùn)算時(shí)間為b1。忽略加窗所需要的計(jì)算時(shí)間。

但進(jìn)一步分析發(fā)現(xiàn),式(7)中的H矩陣具有一定的對(duì)稱特性,具體來(lái)說(shuō),H可表示為

式中,H′是H的子矩陣,為M階方陣。

若將電壓序列u(n)寫(xiě)成如下形式:

其中,uMi表示u(n)中的第i個(gè)周期的采樣點(diǎn)構(gòu)成的序列,為M長(zhǎng)度的列向量。

可見(jiàn),使用傳統(tǒng)方法對(duì)u(n)移相的過(guò)程最終可簡(jiǎn)化成使用H′矩陣與u s(n)做乘積然后將所得時(shí)域序列周期延拓的過(guò)程。實(shí)際應(yīng)用中,只需離線求出H′矩陣,即可實(shí)現(xiàn)信號(hào)的Hilbert變換,這樣可以提高計(jì)算速度、節(jié)省大量的存儲(chǔ)空間。

考慮到應(yīng)用傳統(tǒng)方法時(shí)可以做如上調(diào)整,其計(jì)算時(shí)間縮短為

下文提到的傳統(tǒng)算法都是指調(diào)整后的傳統(tǒng)算法。

改進(jìn)算法不需要離線計(jì)算,但需對(duì)信號(hào)進(jìn)行1次M點(diǎn)FFT和1次M點(diǎn)IFFT,其計(jì)算時(shí)間約為

1次復(fù)乘相當(dāng)于4次實(shí)乘和2次實(shí)加,1次復(fù)加相當(dāng)于2次實(shí)加,則

式(35)變?yōu)?/p>

參考(34)、(37),傳統(tǒng)算法與改進(jìn)算法計(jì)算時(shí)間之差為

若只考慮乘法運(yùn)算:

可見(jiàn),ΔT只與單周期采樣點(diǎn)數(shù)M有關(guān)。

用ΔT/T′o表示改進(jìn)算法相對(duì)于傳統(tǒng)方法節(jié)省計(jì)算量的相對(duì)值,表1給出了其與M之間的關(guān)系(單位:1)。

表1 ΔT/-M關(guān)系表Tab.1 Relationship betweenΔT/ and M

表1 ΔT/-M關(guān)系表Tab.1 Relationship betweenΔT/ and M

642560 a 4096 a 62.512812800 a 16384 a 78.125657344 a 65536 a 87.5512243712 a 262144 a 93.010241007616 a 1048576 a 96.120484104192 a 4194304 a 97.9

由表1知,隨著M的增加,ΔT/T′o的值快速增大,改進(jìn)算法的速度優(yōu)勢(shì)更加明顯。

改進(jìn)算法已應(yīng)用于“三相數(shù)字式多功能測(cè)控電表”的無(wú)功功率測(cè)量中,該電表使用自帶12位A/D轉(zhuǎn)換器的C8051F系列單片機(jī)。實(shí)驗(yàn)表明,改進(jìn)算法相比傳統(tǒng)算法,無(wú)功功率計(jì)算速度快出很多且測(cè)量精度亦有所提高。單周期采樣128點(diǎn),采樣2周波時(shí),傳統(tǒng)算法實(shí)現(xiàn)Hilbert變換所需時(shí)間為89 m s,而改進(jìn)算法只需要42 ms,計(jì)算時(shí)間減少了約52.8%。雖然根據(jù)理論推導(dǎo),在M=128時(shí),改進(jìn)算法節(jié)省時(shí)間量應(yīng)為78.1%,但上述推導(dǎo)只考慮了FFT蝶形運(yùn)算所需時(shí)間,而程序具體實(shí)現(xiàn)中需對(duì)采樣信號(hào)加窗,并進(jìn)行碼位倒置和寄存器賦值等一系列額外運(yùn)算,所以說(shuō)該實(shí)驗(yàn)結(jié)果與理論推導(dǎo)保持一致。

5 結(jié)語(yǔ)

本文對(duì)基于DFT和IDFT實(shí)現(xiàn)H ilbert變換的傳統(tǒng)算法進(jìn)行改進(jìn),提出了一種基于加窗FFT和IFFT的H ilbert變換快速算法。該算法只需對(duì)加窗信號(hào)做1次M點(diǎn)FFT運(yùn)算和1次M點(diǎn)IFFT運(yùn)算(M為單周波采樣點(diǎn)數(shù))即可實(shí)現(xiàn)采樣序列的準(zhǔn)確90°相移。理論分析表明,采用此算法能有效縮短信號(hào)移相計(jì)算時(shí)間,提高無(wú)功測(cè)量的精度。改進(jìn)算法已經(jīng)應(yīng)用于工程實(shí)踐,在應(yīng)用中提高了程序運(yùn)行的效率和無(wú)功測(cè)量的實(shí)時(shí)性、準(zhǔn)確性。

[1] 鄭小平(Zheng Xiaoping).關(guān)于無(wú)功功率的定義及其計(jì)算方法(The definition and arithmetic for the reactive pow er)[J].電測(cè)與儀表(ElectricalM easurement&Instrumentation),2006,43(6):1-4,16.

[2] 陳國(guó)通,吳杰康,張宏亮(Chen Guotong,Wu Jiekang,Zhang Hongliang).無(wú)功功率和電能的移相算法(Phase-shift algorithms for reactive electric power and energy in pow er system s)[J].電力學(xué)報(bào)(Journal of Electric Pow er),2007,22(4):428-431.

[3] 薛蕙,楊仁剛(Xue Hui,Yang Rengang).改進(jìn)的瞬時(shí)無(wú)功和諧波電流檢測(cè)理論(A novel detection theory of instaneous reactive and harmonic current)[J].電力系統(tǒng)及其自動(dòng)化學(xué)報(bào)(Proceedings of the CSUEPSA),2002,14(2):8-11.

[4] 邱海鋒,周浩(Q iu H aifeng,Zhou H ao).電力系統(tǒng)無(wú)功測(cè)量方法綜述(Summary on reactive pow er measurement in pow er system)[J].電測(cè)與儀表(Electrical M easurement&Instrumentation),2007,44(1):5-9.

[5] 鄭常寶,王群京,方斌,等(Zheng Changbao,W ang Qun jing,Fang Bin,et a l).用小波變換和希爾伯特變換測(cè)量無(wú)功功率(Reactive power measurement by wavelet transform and H ilbert transform)[J].系統(tǒng)仿真學(xué)報(bào)(Journalof System Simulation),2005,17(4):822-824.

[6] 黃純,江亞群(Huang Chun,Jiang Yaqun).電功率微機(jī)測(cè)量新算法(Novel algorithm for measuring active and reactive pow er)[J].中國(guó)電機(jī)工程學(xué)報(bào)(Proceedings o f the CSEE),2008,28(4):54-58.

[7] 俎云霄,龐浩,李東霞,等(Zu Yunxiao,Pang H ao,Li Dongxia,et a l).一種基于H ilbert數(shù)字濾波的無(wú)功功率測(cè)量方法(A method of reactive power measurement based on H ilbert digital filtering)[J].電力系統(tǒng)自動(dòng)化(Automation of Electric Pow er System s),2003,27(16):50-52,70.

[8] 龐浩,王贊基,陳建業(yè),等(Pang Hao,Wang Zan ji,Chen Jianye,etal).基于2對(duì)H ilbert移相濾波器的無(wú)功功率測(cè)量方法(Method of reactive pow er measurementbased on tw o pairs of H ilbert digital filters)[J].電力系統(tǒng)自動(dòng)化(Automation of Electric Pow er Systems),2006,30(18):45-48.

[9] AnsariRashid.IIR discrete-time H ilbert transformers[J].IEEE Trans on Acoustics,Speech,and Signal Processing,1987,35(8):1116-1119.

[10]Ansari Rashid.Ellip tic filter design for a class of generalized halfband filters[J].IEEE T rans on A-coustics,Speech,and Signal Processing,1985,33(5):1146-1150.