有限字長對(duì)滑動(dòng)窗DFT穩(wěn)定性的影響研究

劉華吾　胡海兵　邢　巖

(南京航空航天大學(xué)自動(dòng)化學(xué)院 江蘇省新能源發(fā)電與電能變換重點(diǎn)實(shí)驗(yàn)室　南京　210016)

有限字長對(duì)滑動(dòng)窗DFT穩(wěn)定性的影響研究

劉華吾胡海兵邢巖

(南京航空航天大學(xué)自動(dòng)化學(xué)院 江蘇省新能源發(fā)電與電能變換重點(diǎn)實(shí)驗(yàn)室南京210016)

摘要滑動(dòng)窗離散傅里葉變換(DFT)具體數(shù)字實(shí)現(xiàn)過程中，有限字長效應(yīng)引入的誤差會(huì)被累積和放大，易造成系統(tǒng)發(fā)散和不穩(wěn)定。針對(duì)該問題和現(xiàn)象，應(yīng)用隨機(jī)信號(hào)的統(tǒng)計(jì)特性，詳細(xì)分析了信號(hào)量化噪聲、旋轉(zhuǎn)因子量化誤差及數(shù)學(xué)運(yùn)算舍入誤差對(duì)滑動(dòng)窗DFT穩(wěn)定性的影響，并指出舍入誤差所引入的直流分量是導(dǎo)致系統(tǒng)不穩(wěn)定的根本原因。仿真和實(shí)驗(yàn)均驗(yàn)證了理論分析的正確性。為了解決此不穩(wěn)定問題，提出兩種通過交換運(yùn)算次序來消除誤差累積和放大的解決方案，實(shí)驗(yàn)和仿真結(jié)果驗(yàn)證了這兩種方法的有效性。

關(guān)鍵詞：有限字長效應(yīng)滑動(dòng)窗DFT穩(wěn)定性舍入誤差

0引言

離散傅里葉變換(Discrete Fourier Transform，DFT)是數(shù)字信號(hào)處理中的重要工具[1-4]。傳統(tǒng)的傅里葉變換或快速傅里葉變換(FFT)運(yùn)算量大，很難滿足實(shí)時(shí)性要求較高的應(yīng)用場合[5，6]。近年來提出的迭代形式DFT即滑動(dòng)窗DFT[7，8]，由于其運(yùn)算量小及實(shí)時(shí)性好等優(yōu)點(diǎn)，得到了廣泛關(guān)注，被應(yīng)用于實(shí)時(shí)頻譜分析、繼電保護(hù)和電力系統(tǒng)諧波提取等領(lǐng)域[9-14]。

實(shí)際數(shù)字系統(tǒng)中，有限字長效應(yīng)引入的信號(hào)抽樣誤差、濾波器系數(shù)量化誤差、計(jì)算舍入誤差等會(huì)給系統(tǒng)的精度和穩(wěn)定性帶來不良影響?；瑒?dòng)窗DFT的極點(diǎn)分布在離散域的單位圓上，是一個(gè)臨界穩(wěn)定系統(tǒng)，文獻(xiàn)[7，14]指出，旋轉(zhuǎn)因子系數(shù)的量化誤差可能會(huì)將系統(tǒng)極點(diǎn)移出單位圓，致使系統(tǒng)不穩(wěn)定，但文中沒有考慮其他誤差對(duì)系統(tǒng)穩(wěn)定性的影響。文獻(xiàn)[15]指出，滑動(dòng)窗DFT大量乘法的舍入誤差會(huì)在系統(tǒng)中無衰減傳播和累積，文中應(yīng)用隨機(jī)信號(hào)的方差特征分析誤差的放大效應(yīng)，沒有考慮誤差正負(fù)的相互抵消因素。文獻(xiàn)[16]詳細(xì)描述了浮點(diǎn)數(shù)實(shí)現(xiàn)滑動(dòng)窗DFT的不穩(wěn)定現(xiàn)象，但未給出清晰的理論解釋。

針對(duì)上述穩(wěn)定性問題，現(xiàn)有文獻(xiàn)提出了許多解決措施。引入衰減因子法[7，17]將滑動(dòng)窗DFT的極點(diǎn)強(qiáng)行移至單位圓內(nèi)，對(duì)誤差的傳播形成衰減和阻尼，這種方法簡單可靠，但會(huì)引入額外運(yùn)算量且易導(dǎo)致輸出失真?；贕oertzel算法的改進(jìn)DFT[15]和并行DFT[16，18，19]迭代形式的DFT并行運(yùn)算，在每個(gè)變換周期結(jié)束后，用非迭代DFT的輸出更新迭代DFT的數(shù)據(jù)，避免誤差的累積。該方法輸出結(jié)果準(zhǔn)確可靠，但實(shí)現(xiàn)結(jié)構(gòu)復(fù)雜、運(yùn)算量大。優(yōu)化DFT法[16，20]通過優(yōu)化算法運(yùn)算順序避免誤差的累積和放大，具有結(jié)構(gòu)簡單和運(yùn)算量少的優(yōu)點(diǎn)，但需要占用2倍的存儲(chǔ)空間[21，22]，不適合需要同時(shí)提取多個(gè)頻譜的應(yīng)用場合。

本文首先給出了滑動(dòng)窗DFT的數(shù)學(xué)表達(dá)式，指出具體數(shù)字實(shí)現(xiàn)過程中有限字長效應(yīng)引入的誤差，提出了一種穩(wěn)定性分析模型，全面詳細(xì)分析了各種誤差對(duì)系統(tǒng)穩(wěn)定性的影響，得出舍入誤差引入的直流分量是導(dǎo)致不穩(wěn)定現(xiàn)象的原因的結(jié)論，并進(jìn)行了詳細(xì)的仿真和實(shí)驗(yàn)對(duì)比，驗(yàn)證了理論分析的正確性。最后，針對(duì)此穩(wěn)定性問題，基于優(yōu)化算法運(yùn)算順序的思想，提出了兩種簡潔有效的解決措施，并進(jìn)行了仿真和實(shí)驗(yàn)驗(yàn)證。

1滑動(dòng)窗DFT

用采樣頻率N/T0采樣輸入信號(hào)x(t)，生成離散序列{x(n)}，從中選取N點(diǎn)作為樣本。設(shè)第n-1次樣本序列為{x(n-N)，x(n-N+1)，…，x(n-1)}，其m次諧波設(shè)為Xm(n-1)。一個(gè)采樣周期之后，第n次的樣本序列在上次的樣本中移除x(n-N)，加入x(n)，生成{x(n-N+1)，x(n-N+2)，…，x(n)}，該新樣本的m次諧波為Xm(n)。則Xm(n-1)和Xm(n)有以下關(guān)系[23]

Xm(n)=Xm(n-1)+[x(n)-x(n-N)]ejω0(n-1)m

(1)

(2)

θm(n)=∠Xm(n)

(3)

式(1)～式(3)即是滑動(dòng)窗DFT的計(jì)算公式，以提取基波頻譜為例，實(shí)現(xiàn)框圖如圖1所示。

圖1　滑動(dòng)窗DFT實(shí)現(xiàn)框圖Fig.1　Block diagram of recursive DFT

理論上，若輸入是周期信號(hào)且頻率等于1/T0(本文T0=1/50 s)，則x(n)-x(n-N)=0；X1(n)保持上個(gè)采樣周期的值不變。當(dāng)輸入信號(hào)頻率波動(dòng)Δf而DFT采樣頻率保持不變時(shí)，x(n)-x(n-N)≠0，DFT輸出X1(n)也會(huì)有相應(yīng)的脈動(dòng)，但系統(tǒng)穩(wěn)定，不會(huì)隨時(shí)間推進(jìn)而發(fā)散[24]。

然而具體數(shù)字實(shí)現(xiàn)過程中，有限字長的效應(yīng)會(huì)引入各種誤差，包括模擬信號(hào)抽樣產(chǎn)生的量化噪聲、旋轉(zhuǎn)因子的量化誤差、數(shù)學(xué)運(yùn)算的舍入誤差[25]。這些誤差在離散系統(tǒng)中傳遞，經(jīng)滑動(dòng)窗DFT無限次迭代后，易被累積和放大，最終導(dǎo)致系統(tǒng)輸出發(fā)散和不穩(wěn)定。

2有限字長對(duì)滑動(dòng)窗DFT的穩(wěn)定性影響

為分析滑動(dòng)窗DFT的穩(wěn)定性問題，可將圖1的實(shí)現(xiàn)框圖簡化為圖2所示的分析模型。

圖2　滑動(dòng)窗DFT的簡化模型Fig.2　Simplified analysis model of the sliding DFT

圖2中M表示運(yùn)算

vq(n)=[x(n)-x(n-N)]ejω0(n-1)

(4)

式(4)對(duì)應(yīng)的離散域傳遞函數(shù)是一個(gè)有限脈沖響應(yīng)(FIR)系統(tǒng)，無輸出到輸入的反饋，不會(huì)累積和放大誤差，也不會(huì)導(dǎo)致系統(tǒng)輸出不穩(wěn)定。而圖2中的數(shù)據(jù)累加部分包含輸出到輸入的反饋，本質(zhì)是一個(gè)積分器，對(duì)直流分量有無窮大增益。這意味著，如果輸入vq(n)含有直流分量，則會(huì)被無限累積和放大而導(dǎo)致輸出發(fā)散?？梢姡瑒?dòng)窗DFT的迭代形式是其實(shí)際數(shù)字過程中存在穩(wěn)定性問題的前提。相反的，如果系統(tǒng)輸出不穩(wěn)定，必然是由于vq(n)含有直流分量。因而下文將圍繞是否會(huì)在vq(n)引入直流分量，對(duì)有限字長效應(yīng)引入的各種誤差展開詳細(xì)分析。

滑動(dòng)窗DFT可采用定點(diǎn)數(shù)或浮點(diǎn)數(shù)實(shí)現(xiàn)，穩(wěn)定性分析方法類似，本文針對(duì)定點(diǎn)數(shù)場合。由圖1可知，滑動(dòng)窗DFT提取頻譜的實(shí)部和虛部對(duì)稱，選取實(shí)部的運(yùn)算過程作為研究對(duì)象，考慮數(shù)字實(shí)現(xiàn)過程中所有可能的誤差，可得到如圖3所示的穩(wěn)定性分析模型。模型被劃分為信號(hào)量化、梳狀濾波器、旋轉(zhuǎn)因子相乘、數(shù)據(jù)舍入和累加等五部分，其中es(n)表示信號(hào)抽樣的量化噪聲，eR(n)表示旋轉(zhuǎn)因子的量化誤差，ev(n)表示乘法運(yùn)算的舍入誤差。

圖3　滑動(dòng)窗DFT的穩(wěn)定性分析模型Fig.3　Analysis model of the practically implemented sliding DFT

2.1信號(hào)量化

滑動(dòng)窗DFT的輸入vs(t)經(jīng)采樣量化成離散序列vt(n)，實(shí)際的采樣量化過程可等效成理想的抽樣疊加量化噪聲[25]

vt(n)=vs(n)+es(n)

(5)

式中，vs(n)為理想抽樣結(jié)果；es(n)為量化噪聲，是一個(gè)平穩(wěn)隨機(jī)序列[25，26]，其統(tǒng)計(jì)均值μes=μ(和時(shí)間無關(guān)的常數(shù))。

量化噪聲es(n)經(jīng)過梳狀濾波器、旋轉(zhuǎn)因子相乘和數(shù)據(jù)舍入等模塊最終傳遞到vq(n)，它對(duì)系統(tǒng)穩(wěn)定性的影響在2.5節(jié)論述。

2.2梳狀濾波器

梳狀濾波器離散域傳遞函數(shù)為

H1(z)=1-z-N

(6)

其頻率響應(yīng)為

H1(ejω)=1-e-jNω

(7)

易知，當(dāng)ω=kω0=2kπ/N(k=0，±1，…，±N-1)時(shí)，H1(ejω)=0，這表明梳狀濾波器完全濾除直流分量和頻率為ω0及其整數(shù)倍諧波。

信號(hào)量化模塊的輸出vt(n)，經(jīng)梳狀濾波器后為vr(n)

vr(n)=vt(n)h1(n)=[vs(n)+es(n)]h1(n)

=vs(n)h1(n)+es(n)h1(n)

=vrs(n)+vrv(n)

(8)式中，h1(n)=δ(n)-δ(n-N)，表示H1(z)的單位序列響應(yīng)；vrs(n)是輸入確定信號(hào)vs(n)經(jīng)梳狀濾波器的響應(yīng)輸出，由梳狀濾波器的特性知，vrs(n)不含有直流分量和頻率為ω0及其整數(shù)倍諧波分量；vrv(n)是一個(gè)隨機(jī)信號(hào)，由量化噪聲es(n)經(jīng)梳狀濾波器的輸出，其統(tǒng)計(jì)平均值[25]

(9)

式(9)表明隨機(jī)信號(hào)vrv(n)的統(tǒng)計(jì)均值為零，和量化噪聲es(n)的統(tǒng)計(jì)數(shù)字特征μ無關(guān)。

2.3旋轉(zhuǎn)因子相乘

定點(diǎn)數(shù)表示旋轉(zhuǎn)因子時(shí)，需要先定標(biāo)(乘以2Q，然后取整)，定標(biāo)取整引入的量化誤差記為eR(n)，如圖3所示。eR(n)是一個(gè)平穩(wěn)的隨機(jī)序列[25，26]，統(tǒng)計(jì)均值E[eR(k)]=ε(常數(shù))。

vr(n)經(jīng)旋轉(zhuǎn)因子相乘后的vp(n)為

vp(n)=vr(n){2Qcos[ω0(n-1)]+2QeR(n)]}

(10)

2.4數(shù)據(jù)舍入和累加

旋轉(zhuǎn)因子相乘后的舍入誤差記為ev(n)。由圖3可知，最終累加部分的輸入vq(n)為

(11)

結(jié)合式(8)和式(10)，將式(11)展開得

vq(n)=y1(n)+y2(n)+y3(n)+y4(n)+y5(n)

(12)

式中

y1(n)=vrs(n)cos[ω0(n-1)]

y2(n)=vrs(n)eR(n)

y3(n)=vrv(n)cos[ω0(n-1)]

y4(n)=vrv(n)eR(n)

y5(n)=ev(n)

(13)

式中，y1(n)是確定信號(hào)；y2(n)～y5(n)是隨機(jī)信號(hào)，表示旋轉(zhuǎn)因子量化誤差、抽樣量化噪聲和數(shù)字運(yùn)算舍入誤差等對(duì)滑動(dòng)窗DFT的穩(wěn)定性影響。

2.5穩(wěn)定性分析

y1(n)是確定信號(hào)，2.2節(jié)指出vrs(n)不含頻率為ω0分量，則乘積vrs(n)cos[ω0(n-1)]必然沒有直流分量(詳細(xì)數(shù)學(xué)證明見附錄A)，不影響系統(tǒng)穩(wěn)定性。

隨機(jī)信號(hào)y2(n)～y5(n)的直流分量通常用統(tǒng)計(jì)均值表示[25]。

y2(n)的統(tǒng)計(jì)均值為

E[y2(n)]=E[vrs(n)eR(n)]=vrs(n)E[eR(n)]

=εvrs(n)

(14)

由于vrs(n)是無直流分量的確定信號(hào)，因而y2(n)不會(huì)影響系統(tǒng)的穩(wěn)定性(詳細(xì)分析見附錄B)。

y3(n)的統(tǒng)計(jì)均值為

E[y3(n)]=E{vrv(n)cos[ω0(n-1)]}

=cos[ω0(n-1)]E[vrv(n)]

=μrvcos[ω0(n-1)]

(15)

由式(9)知μrv=0，因而E[y3(n)]=0，從而y3(n)無直流分量，不影響系統(tǒng)的穩(wěn)定性。

y4(n)的統(tǒng)計(jì)均值為

E[y4(n)]=E[vrv(n)eR(n)]

(16)

考慮到隨機(jī)信號(hào)vrv(n)和eR(n)不相關(guān)，有

E[y4(n)]=E[vrv(n)]E[eR(n)]=μrvE[eR(n)]=0

(17)

可見，y4(n)同樣無直流分量，不影響系統(tǒng)的穩(wěn)定性。

y5(n)是乘法運(yùn)算的舍入誤差，其統(tǒng)計(jì)均值和舍入方式有關(guān)。直接截尾和四舍五入是兩種最常用的舍入方式，直接截尾法將尾數(shù)直接舍去，實(shí)現(xiàn)最簡單。四舍五入法根據(jù)尾數(shù)的最高位數(shù)字來決定是否進(jìn)位，實(shí)現(xiàn)相對(duì)復(fù)雜。

y5(n)的實(shí)際統(tǒng)計(jì)均值與抽樣樣本大小、隨機(jī)實(shí)驗(yàn)次數(shù)以及數(shù)值的波動(dòng)范圍等有關(guān)。這和y3(n)、y4(n)不同，由式(5)～式(17)可知，y3(n)、y4(n)的統(tǒng)計(jì)均值直接由vrv(n)的統(tǒng)計(jì)特征決定；梳狀濾波器完全濾除直流分量的特性保證了vrv(n)的統(tǒng)計(jì)均值μrv=0(式(9))而和隨機(jī)過程無關(guān)。

采用直接截尾法時(shí)，雖然無法準(zhǔn)確求取其統(tǒng)計(jì)均值，但每次舍入時(shí)均將尾數(shù)直接舍去，y5(n)<0恒成立，從而

E[y5(n)]<0

(18)

這說明直接截尾法的舍入誤差會(huì)引入負(fù)的直流分量，必然導(dǎo)致系統(tǒng)不穩(wěn)定。

采用四舍五入法時(shí)，每次舍入的尾數(shù)可正可負(fù)，當(dāng)抽樣樣本無限大、隨機(jī)實(shí)驗(yàn)次數(shù)無限多、乘法運(yùn)算的數(shù)值均勻分布時(shí)，y5(n)的統(tǒng)計(jì)均值趨于零。這些條件在實(shí)際系統(tǒng)中顯然很難保證。相反的，很可能在一段時(shí)間內(nèi)E[y5(n)]>0，而在另一段時(shí)間內(nèi)E[y5(n)]<0。y5(n)統(tǒng)計(jì)均值的這種正負(fù)波動(dòng)將致使滑動(dòng)窗DFT輸出頻譜的實(shí)虛部往+∞或-∞方向發(fā)散，引起誤差累積和放大，易導(dǎo)致輸出不可靠。

總結(jié)上述分析，得到以下結(jié)論：

1)輸入信號(hào)波動(dòng)、A-D量化噪聲和旋轉(zhuǎn)因子系數(shù)量化誤差不影響系統(tǒng)穩(wěn)定性。

2)直接截尾法會(huì)引入負(fù)直流偏置，必然導(dǎo)致系統(tǒng)輸出的頻譜實(shí)虛部往-∞方向發(fā)散。

3)四舍五入法引入的誤差，在實(shí)際系統(tǒng)中很難保證其統(tǒng)計(jì)均值為零，易導(dǎo)致系統(tǒng)輸出不可靠。

2.6穩(wěn)定性問題的仿真和實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證上文分析的正確性，搭建了仿真模型和實(shí)驗(yàn)平臺(tái)，對(duì)滑動(dòng)窗DFT進(jìn)行詳細(xì)仿真和實(shí)驗(yàn)。

圖4為Matlab/Simulink中的仿真模型。量化系數(shù)取214/5，輸入信號(hào)疊加白噪聲，模擬實(shí)際采樣電路和A-D量化引入的噪聲，量化模塊的步長取5/2-14。計(jì)算過程采用Q14定標(biāo)，旋轉(zhuǎn)因子定標(biāo)值采用直接尾數(shù)截?cái)喾?。?jì)算模塊尾數(shù)的舍入，為對(duì)比研究，分別采用四舍五入或直接截尾法。

圖4　滑動(dòng)窗DFT Matlab仿真模型Fig.4　Matlab simulation model of the sliding DFT

搭建了實(shí)驗(yàn)平臺(tái)，平臺(tái)采用現(xiàn)場可編程門陣列(FPGA，EP4CE15F17C8)作為控制器，任意函數(shù)發(fā)生器(DG4072)作為信號(hào)輸入源。輸入信號(hào)通過14bit-A-D(ADS7945)輸入至FPGA，A-D的基準(zhǔn)由芯片(AD586)生成。FPGA計(jì)算結(jié)果經(jīng)16bit-D-A(DAC8501)轉(zhuǎn)換成模擬量，輸出至示波器。實(shí)驗(yàn)中滑動(dòng)窗DFT實(shí)現(xiàn)架構(gòu)和圖4仿真模型一致。

仿真和實(shí)驗(yàn)中，用滑動(dòng)窗DFT提取輸入信號(hào)的基波頻譜，DFT變換點(diǎn)數(shù)N為128，采樣頻率為6 400 Hz。

1)實(shí)例一：頻率49.5 Hz正弦波輸入，旋轉(zhuǎn)因子乘法運(yùn)算尾數(shù)舍入采用直接截尾法。

仿真中，輸入信號(hào)為頻率49.5 Hz的正弦波

vs(t)=sin(49.5·2πt)

(19)

實(shí)驗(yàn)中，由信號(hào)發(fā)生器生成幅值1 V、頻率49.5 Hz的正弦波，經(jīng)A-D量化后輸入至滑動(dòng)窗DFT。

圖5為系統(tǒng)運(yùn)行初始的仿真和實(shí)驗(yàn)結(jié)果，其中A和θ分別表示DFT輸出的幅值和相角。

圖5　滑動(dòng)窗DFT的仿真和實(shí)驗(yàn)結(jié)果(實(shí)例一)Fig.5　Simulation and experimental results of the recursive DFT under conditions of the first case

從圖5中可明顯看出，DFT輸出的幅值迅速發(fā)散，相角也隨之明顯畸變。頻譜的實(shí)虛部含有幅值保持不變的交流分量和往-∞發(fā)散的直流分量。交流分量由輸入信號(hào)為0.5 Hz的頻率偏差引起，不影響系統(tǒng)穩(wěn)定性。直流分量由舍入誤差引入，經(jīng)DFT輸出累積部分迭代后，迅速累積和放大，致使系統(tǒng)發(fā)散。

2)實(shí)例二：頻率49.5 Hz正弦波輸入，旋轉(zhuǎn)因子乘法運(yùn)算尾數(shù)舍入采用四舍五入法。

在該實(shí)例中，滑動(dòng)窗DFT的輸入和實(shí)例一保持一致，僅將舍入方式改成四舍五入法。

仿真模型和實(shí)驗(yàn)平臺(tái)運(yùn)行10 s、1 h和20 h結(jié)果分別如圖6a～圖6f所示。可以看出，仿真和實(shí)驗(yàn)剛開始時(shí)，DFT輸出的幅值和相角是準(zhǔn)確可靠的，但隨著時(shí)間的推進(jìn)，幅值A(chǔ)上疊加的低頻脈動(dòng)被逐漸放大，只是發(fā)散的速度非常緩慢。實(shí)驗(yàn)20 h后，低頻脈動(dòng)的幅值約0.1 V，比仿真的結(jié)果大，這是由于實(shí)驗(yàn)和仿真中A-D采樣量化噪聲的不一致造成的。此外，受四舍五入引入的舍去誤差統(tǒng)計(jì)均值的正負(fù)波動(dòng)影響，該低頻脈動(dòng)不是一直被放大，存在被縮小的時(shí)候，因文章篇幅限制，此處略去相應(yīng)的實(shí)驗(yàn)和仿真波形?？偠灾?，即便采用較為復(fù)雜的四舍五入舍入方式，滑動(dòng)窗DFT輸出依然不可靠。

圖6　滑動(dòng)窗DFT的仿真和實(shí)驗(yàn)結(jié)果(實(shí)例二)Fig.6　Simulation and experimental results of the recursive DFT under conditions of the second case

3解決方法

3.1改進(jìn)DFT法

從圖7可看出，改進(jìn)DFT法的實(shí)現(xiàn)簡潔方便，不會(huì)引入額外的存儲(chǔ)空間和運(yùn)算量，僅會(huì)使數(shù)據(jù)累加部分的加法運(yùn)算位寬擴(kuò)大Q位。

圖7　改進(jìn)DFT法Fig.7　Modified DFT technique

3.2改進(jìn)并行DFT法

圖2分析指出，滑動(dòng)窗DFT穩(wěn)定的前提是數(shù)據(jù)累加部分的輸入不含直流分量。3.1節(jié)提出的改進(jìn)DFT法的核心是通過交換運(yùn)算次序?qū)ⅰ爸绷鞣至康漠a(chǎn)生源——數(shù)據(jù)舍入”移至DFT最終輸出，避免其進(jìn)入數(shù)據(jù)累加部分。此方法適合定點(diǎn)數(shù)運(yùn)算，但在浮點(diǎn)數(shù)計(jì)算場合，由于無法直接和方便地控制數(shù)據(jù)舍入，其不再適用。在這種情況下，本文采用另外一種思路：濾除直流分量。

為濾除直流分量，保證系統(tǒng)穩(wěn)定性，最直接的方法是在數(shù)據(jù)累加部分之前插入額外的低阻濾波器，但這會(huì)增加運(yùn)算量且可能影響DFT提取頻譜的精度。而從對(duì)圖3的分析可知，滑動(dòng)窗DFT中已存在可以完全濾除直流分量的梳狀濾波器，因而無須再引入額外濾波器，僅須調(diào)整運(yùn)算次序，將其移至數(shù)據(jù)累加部分之前，如圖8所示。

圖8　優(yōu)化并行DFT法Fig.8　Optimized parallel DFT technique

優(yōu)化并行DFT法中，梳狀濾波器緊挨數(shù)據(jù)累加部分，相比于滑動(dòng)窗DFT的傳統(tǒng)實(shí)現(xiàn)方式(見圖1或圖3)，其不僅可剔除輸入信號(hào)、A-D采樣量化等引入的直流偏置，還可以濾除數(shù)據(jù)舍入引入的直流分量，從而無直流偏置進(jìn)入數(shù)據(jù)累加部分，保證系統(tǒng)穩(wěn)定性。

優(yōu)化并行DFT法可靠有效，其僅將梳狀濾波器移至數(shù)據(jù)累加之前，相比于改進(jìn)DFT法，無須對(duì)數(shù)據(jù)舍入部分作特別處理，因此適用于定點(diǎn)數(shù)和浮點(diǎn)數(shù)運(yùn)算場合。然而，這種方法占用2N存儲(chǔ)空間，且當(dāng)需要同時(shí)提取多次頻譜時(shí)，占用的存儲(chǔ)空間將會(huì)成倍數(shù)增加，因而需要進(jìn)一步優(yōu)化該方法。

圖8所示優(yōu)化并行DFT法對(duì)應(yīng)的數(shù)學(xué)表達(dá)式為

X1(n)=X1(n-1)+vt(n)ejω0(n-1)-vt(n-N)ejω0(n-1-N)

(20)

考慮到

ejω0(n-1-N)=ejω0(n-1)-2π=ejω0(n-1)

(21)

因而，式(20)可簡化為

X1(n)=X1(n-1)+vt(n)ejω0(n-1)-vt(n-N)ejω0(n-1)

(22)

式(22)對(duì)應(yīng)的具體實(shí)現(xiàn)框圖如圖9所示。改進(jìn)并行DFT法和圖8所示的優(yōu)化并行DFT法等價(jià)，只是前者利用旋轉(zhuǎn)因子的周期對(duì)稱性質(zhì)，進(jìn)一步優(yōu)化了具體實(shí)現(xiàn)過程。此外，該算法雖然額外增加了兩次乘法和一次減法，但僅占用N存儲(chǔ)空間，且存儲(chǔ)空間不會(huì)隨提取頻譜次數(shù)的增加而增加。因而該方法特別適用于存儲(chǔ)資源有限或需要同時(shí)提取多個(gè)頻譜的應(yīng)用場合。

圖9　改進(jìn)并行DFT法Fig.9　Modified parallel DFT technique

3.3實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證本文提出的兩種算法的可行性，進(jìn)行了實(shí)驗(yàn)研究?？紤]到圖7所示的改進(jìn)DFT法僅適用于定點(diǎn)數(shù)計(jì)算，而圖9的改進(jìn)并行DFT法適用于定點(diǎn)數(shù)和浮點(diǎn)數(shù)，因此本文驗(yàn)證滑動(dòng)窗DFT的3種具體實(shí)現(xiàn)方式如下：

1)采用定點(diǎn)數(shù)運(yùn)算和直接截尾的舍入方式，實(shí)現(xiàn)改進(jìn)DFT法。

2)采用定點(diǎn)數(shù)運(yùn)算和直接截尾的舍入方式，實(shí)現(xiàn)改進(jìn)并行DFT法。

3)運(yùn)用FPGA中的自帶浮點(diǎn)數(shù)運(yùn)算IP核，采用單精度(single)型計(jì)算方式，實(shí)現(xiàn)改進(jìn)并行DFT法。

實(shí)驗(yàn)研究主要分為兩部分：第一部分以信號(hào)發(fā)生器作為信號(hào)輸入源，第二部分將本文提出的算法應(yīng)用于有源濾波器(APF)中。

1)以信號(hào)發(fā)生器作為信號(hào)輸入源

該部分的實(shí)驗(yàn)平臺(tái)和2.6節(jié)的一致，如圖10所示，以Altera公司的FPGA(EP4CE15F17C8)作為核心控制單元，任意函數(shù)發(fā)生器(DG4072)作為信號(hào)輸入源。

圖10　基于信號(hào)發(fā)生器的實(shí)驗(yàn)平臺(tái)Fig.10　Experimental setup based on signal generator

信號(hào)發(fā)生器輸出1 V/49.5 Hz正弦波，在FPGA中分別采用上述滑動(dòng)窗DFT的3種實(shí)現(xiàn)方式，提取50 Hz基波，驗(yàn)證算法的可靠性。平臺(tái)運(yùn)行約20 h后，3種算法的實(shí)驗(yàn)波形一致，如圖11所示，其中vs為輸入信號(hào)，vf為DFT提取的基波時(shí)域信號(hào)，A和θ分別為DFT輸出的幅值和相角。從圖中可看出，輸入信號(hào)頻率偏差、舍入和量化誤差等因素均未導(dǎo)致DFT輸出發(fā)散，系統(tǒng)穩(wěn)定可靠。

圖11　滑動(dòng)窗DFT運(yùn)行20 h的實(shí)驗(yàn)結(jié)果Fig.11　20 h-post simulation and experimental results of the recursive DFT

2)并聯(lián)型選擇性APF實(shí)驗(yàn)研究

該部分實(shí)驗(yàn)將滑動(dòng)窗DFT應(yīng)用到實(shí)驗(yàn)室20 A單相并聯(lián)型選擇性APF的樣機(jī)平臺(tái)上，如圖12所示。實(shí)驗(yàn)平臺(tái)采用FPGA作為核心控制器，主電路開關(guān)器件為F4-50R06W1E3。FPGA控制器利用滑動(dòng)窗DFT提取非線性負(fù)載中的各次諧波，然后控制APF輸出相應(yīng)的電流以抵消電網(wǎng)中的諧波，提高電能質(zhì)量。APF電流的控制采用比例諧振調(diào)節(jié)器(PR)，以實(shí)現(xiàn)無靜差諧波基準(zhǔn)跟蹤。

圖12　APF實(shí)驗(yàn)平臺(tái)Fig.12　APF experimental setup

將上述滑動(dòng)窗DFT的3種實(shí)現(xiàn)方式(定點(diǎn)數(shù)運(yùn)算的改進(jìn)DFT法，定點(diǎn)數(shù)和浮點(diǎn)數(shù)運(yùn)算的改進(jìn)并行DFT法)分別在實(shí)驗(yàn)平臺(tái)上驗(yàn)證。3種實(shí)現(xiàn)方式下APF輸出波形均一致，如圖13所示。圖中，iL為非線性負(fù)載的諧波波形，iC為APF輸出電流波形，iS為電網(wǎng)電流波形。圖13a中，APF僅補(bǔ)償3次諧波的波形，電網(wǎng)電流中3次諧波成分從補(bǔ)償前的5.32 A降到0.05 A；圖13b為APF補(bǔ)償3～21次諧波的波形，電網(wǎng)電流THD從補(bǔ)償前的55.96%降到6.03%。圖13說明了本文提出的方法能夠準(zhǔn)確提取諧波，驗(yàn)證了算法的可行性。

圖13　APF實(shí)驗(yàn)結(jié)果Fig.13　APF experimental results

為進(jìn)一步驗(yàn)證算法的穩(wěn)定性，以同時(shí)補(bǔ)償3～21次諧波為例，將滑動(dòng)窗DFT的3種實(shí)現(xiàn)方式下的APF分別進(jìn)行了約10 h的拷機(jī)實(shí)驗(yàn)，拷機(jī)過程中APF能夠穩(wěn)定運(yùn)行，系統(tǒng)最終輸出和圖13所示保持一致，說明了算法的可靠性。

以上實(shí)驗(yàn)結(jié)果表明，采用本文提出的改進(jìn)DFT法和改進(jìn)并行DFT法，能夠準(zhǔn)確提取諧波頻譜，有效避免有限字長效應(yīng)引起的穩(wěn)定性問題，保證系統(tǒng)可靠輸出，從而驗(yàn)證了算法的可行性和有效性。

4結(jié)論

深入分析了有限字長效應(yīng)對(duì)滑動(dòng)窗DFT穩(wěn)定性的影響。理論分析和仿真實(shí)驗(yàn)表明，乘法運(yùn)算的舍入誤差會(huì)在系統(tǒng)數(shù)據(jù)累加部分引入直流分量，導(dǎo)致輸出發(fā)散，而輸入信號(hào)量化噪聲和旋轉(zhuǎn)因子系數(shù)量化誤差等均不會(huì)影響系統(tǒng)穩(wěn)定性。

針對(duì)滑動(dòng)窗DFT的不穩(wěn)定問題，基于優(yōu)化運(yùn)算順序的思想，提出了改進(jìn)DFT法和改進(jìn)并行DFT法，前者占用資源最少，適用于定點(diǎn)數(shù)運(yùn)算，后者適用于浮點(diǎn)或定點(diǎn)數(shù)運(yùn)算。此外，這兩種方法均不會(huì)占用額外存儲(chǔ)空間，適用于需要同時(shí)提取多個(gè)頻譜的應(yīng)用場合。

附錄

附錄 A：證明y1(n)=vrs(n)cos[ω0(n-1)]不含直流分量。

由三角函數(shù)性質(zhì)可知

y1(n)=vrs(n)cos[ω0(n-1)]

(A1)

設(shè)vrs(n)和y1(n)對(duì)應(yīng)的傅里葉變換分別為Vrs(jω)和Y1(jω)，則由式(A1)知

(A2)

從而y1(n)的直流分量

(A3)

本文中2.2節(jié)已經(jīng)指出，vrs(n)不含頻率為ω0及其整數(shù)倍諧波分量，即Vrs(j0)=Vrs(-j0)=0，從而Y1(j0)=0，表明y1(n)無直流分量。

附錄 B：證明隨機(jī)信號(hào)y2(n)=vrs(n)eR(n)不影響滑動(dòng)窗DFT的穩(wěn)定性。

設(shè)y2(n)經(jīng)滑動(dòng)窗DFT的累加部分后輸出為R2(n)，則

(B1)

本文中第2節(jié)指出，輸出累加部分僅對(duì)直流分量增益無窮大，因而可通過分析R2(n)的直流分量來評(píng)估y2(n)對(duì)系統(tǒng)穩(wěn)定性的影響。

隨機(jī)信號(hào)y2(n)的統(tǒng)計(jì)均值

(B2)

式中，ε是一個(gè)常數(shù)，表示eR(n)的統(tǒng)計(jì)均值?？紤]到vrs(n)是一個(gè)無直流分量的確定信號(hào)，∑vrs(k)和E[R2(n)]必然是一個(gè)有界值，不會(huì)隨時(shí)間推進(jìn)(n變大)而發(fā)散，故而y2(n)不影響系統(tǒng)穩(wěn)定性。

參考文獻(xiàn)

[1]王堯，李奎，任伯飛，等.基于全相位傅里葉變換的磁調(diào)制交直流漏電電流檢測方法[J].電工技術(shù)學(xué)報(bào)，2015，30(18)：254-260.

Wang Yao，Li Kui，Ren Bofei，et al.Study of fluxgate current detecting method for AC-DC earth leakage current based on apFFT[J].Transactions of China Electrotechnical Society，2015，30(18)：254-260.

[2]孫曉云，同向前，高鑫.柔性直流輸電系統(tǒng)中IGBT閥的故障診斷方法[J].電工技術(shù)學(xué)報(bào)，2014，29(8)：235-241.

Sun Xiaoyun，Tong Xiangqian，Gao Xin，et al.Research on the fault diagnosis of IGBT valve in VSC-HVDC[J].Transactions of China Electrotechnical Society，2014，29(8)：235-241.

[3]曾博，唐求，卿柏元，等.基于Nuttall自卷積窗的改進(jìn)FFT譜分析方法[J].電工技術(shù)學(xué)報(bào)，2014，29(7)：59-65.

Zeng Bo，Tang Qiu，Qing Baiyuan，et al.Spectral analysis method based on improved FFT by Nuttall selfconvolution window[J].Transactions of China Electrotechnical Society，2014，29(9)：59-65.

[4]王興貴，劉正英.基于載波變幅移相調(diào)制方法的串聯(lián)型微網(wǎng)功率平衡控制[J].電力系統(tǒng)保護(hù)與控制，2015，43(13)：38-44.

Wang Xinggui，Liu Zhengying.Series micro-grid power balance control based on carrier amplitude variation and[J].Power System Protection and Control，2015，43(13)：38-44.

[5]徐建，張語勍，李彥斌，等.短時(shí)傅里葉變換和S變換用于檢測電壓暫降的對(duì)比研究[J].電力系統(tǒng)保護(hù)與控制，2014，12(16)：44-48.Xu Jian，Zhang Yuqing，Li Yanbin，et al.Comparative study of STFT and S transform on detecting voltage sag[J].Power System Protection and Control，2014，12(16)：44-48.

[6]郁祎琳，徐永海，劉曉博.滑窗迭代DFT的諧波電流檢測方法[J].電力系統(tǒng)保護(hù)與控制，2011，39(13)：78-90.

Yu Yilin，Xu Yonghai，Liu Xiaobo.Study of harmonic current detection based on sliding-window iterative algorithm of DFT[J].Power System Protection and Control，2011，39(13)：78-90.

[7]Jacobsen E，Lyons R.The sliding DFT[J].IEEE Signal Processing Magazine，2003，20(3)：74-80.

[8]Jacobsen E，Lyons R.An update to the sliding DFT[J].IEEE Signal Processing Magazine，2004，21(1)：110-111.

[9]王宏偉.基于傅里葉變換的數(shù)字信道化及相關(guān)技術(shù)[D].西安：西安電子科技大學(xué)，2010.

[10]Carlos M O，Ignacio C，Sebastian M，et al.Harmonics measurement with a modulated sliding discrete Fourier transform algorithm[J].IEEE Transactions on Instrumen-tation and Measurement，2014，63(4)：781-793.

[11]周柯，羅安，湯賜，等.一種大功率混合注入式有源電力濾波器的工程應(yīng)用[J].中國電機(jī)工程學(xué)報(bào)，2007，27(22)：80-86.

Zhou Ke，Luo An，Tang Ci，et al.High-power hybrid injection active power filter’s engineering application[J].Proceedings of the CSEE，2007，27(22)：80-86.

[12]Ni Ruoshui，Li Yunwei，Zhang Ye，et al.Virtual impedance based selective harmonic compensation(VI-SHC)PWM for current source rectifiers[J].IEEE Transactions on Power Electronics，2014，29(7)：3346-3356.

[13]Reza M S，Ciobotaru M，Agelidis V G.Accurate estimation of single-phase grid voltage parameters under distorted conditions[J].IEEE Transactions on Power Delivery，2014，29(3)：138-1146.

[14]Kim J H，Chang T G.Analytic derivation of the finite word-length effect of the twiddle factors in recursive implementation of the sliding DFT[J].IEEE Transactions on Signal Processing，2000，48(5)：1485-1488.

[15]王宏偉.滑動(dòng)離散傅立葉算法輸出穩(wěn)定性研究[J].電波科學(xué)學(xué)報(bào)，2012，27(4)：773-779.

Wang Hongwei.Output stabilization of sliding discrete Fourier transform algorithm[J].Chinese Journal of Radio Science，2012，27(4)：773-779.

[16]Darwish H A，F(xiàn)ikri M.Practical considerations for recursive DFT implementation in numerical relays[J].IEEE Transactions on Power Delivery，2007，22(1)：42-49.

[17]Duda K.Accurate，guaranteed stable，sliding discrete fourier transform[J].IEEE Signal Processing Magazine，2010，27(6)：124-127.

[18]Neves F A S，Souza H E P，Cavalcanti M C，et al.Digital filters for fast harmonic sequence component separation of unbalanced and distorted three-phase signals[J].IEEE Transactions on Industrial Electronics，2012，59(10)：3847-3859.

[19]Neves S，Arcanjo C，Azevedo S，et al.The SVFT-Based Control[J].IEEE Transactions on Industrial Electronics，2014，61(8)：4152-4160.

[20]Borisov K，Ginn H，Chen G.A computationally efficient RDFT-based reference signal generator for active compensators[J].IEEE Transactions on Power Delivery，2009，24(4)：2396-2404.

[21]Ginn H L，Chen Guangda.Digital control method for grid-connected converters supplied with nonideal voltage[J].IEEE Transactions on Industrial Informatics，2014，10(1)：127-136.

[22]Chen Guangda，Jiang Yingying，Zhou Haiguo.Practical issues of recursive DFT in active power filter based on CPC power theory[C]//Proceedings of Power and Energy Engineering Conference，Wuhan，2009：1-5.

[23]McGrath B P，Holmes D G，Galloway J J H.Power converter line synchronization using a discrete Fourier transform(DFT)based on a variable sample rate[J].IEEE Transactions on Power Electronics，2005，20(4)：877-884.

[24]Yang Junzhe，Liu Chihwen.A precise calculation of power system frequency[J].IEEE Transactions on Power Delivery，2001，16(3)：361-366.

[25]胡廣書.數(shù)字信號(hào)處理——理論、算法與實(shí)現(xiàn)[M].3版.北京：清華大學(xué)出版社，2012：406-425.

[26]常建平，李海林.隨機(jī)信號(hào)分析[M].北京：科學(xué)出版社，2006：82-94.

Research of Finite-Word-Length Effects on the Stability of the Sliding DFT

Liu HuawuHu HaibingXing Yan

(Jiangsu Key Laboratory of New Energy Generation and Power Conversion Nanjing University of Aeronautics and AstronauticsNanjing210016China)

AbstractThe sliding discrete Fourier transform(DFT)can experience instability issues due to the accumulation and amplification of the errors introduced by the finite-word-length effects in real practical discrete systems.In this paper,the quantization noise of A-D,the quantization error of the twiddle factors,and the rounding error of arithmetic operations are taken into consideration to analyze their influences on the stability of the sliding DFT based on the statistical properites of random signals.Based on the detailed analysis,it is revealed that the rounding error of arithmetic operations is the key reason to cause the unstable phenomenon.Both simulation and experimental results validate the theoretical analysis of instability issue of the sliding DFT.To deal with this instability problem,two methods are proposed to eliminate the error accumulation and amplification through swapping the calculation sequences of the sliding DFT.Extensive simulations and experiments are provided to verify the effectiveness of the proposed solutions.

Keywords：Finite-word-length effects，sliding DFT，stability，rounding errors

收稿日期2015-04-02改稿日期2015-07-20

作者簡介E-mail：liuhuawu@nuaa.edu.cn(通信作者) E-mail：huhaibing@nuaa.edu.cn

中圖分類號(hào)：TM935

江蘇省產(chǎn)學(xué)研聯(lián)合創(chuàng)新資金前瞻性研究項(xiàng)目資助(BY2014003-12)。

劉華吾男，1989年生，博士研究生，研究方向?yàn)殡娏﹄娮优c電力傳動(dòng)。

胡海兵男，1973年生，教授，研究方向?yàn)殡娏﹄娮訑?shù)字控制裝置及電力電子系統(tǒng)集成。