基于MATLAB的全波傅氏算法仿真

1、文章編號 :1005 - 0523 (2006) 01 - 0113 - 04基于 MATLAB 的全波傅氏算法仿真何志勤 ,樊江濤( 華東交通大學 電氣與電子工程學院 ,江西 南昌 330013)摘要 :運用 MATLAB 對傳統(tǒng)全波傅氏算法和 2 種改進算法進行仿真 ,通過對 3 種算法頻譜圖的比較分析 ,說明 2 種改進傅氏算法能夠有效濾掉故障電流中衰減的直流分量 ,從而獲得更為精確的基波和諧波的幅值 .關(guān) 鍵 詞 :傅氏算法 ;衰減直流分量 ;MATLAB中圖分類號 : TM774文獻標識碼 :A運用全波傅氏算法 , 有 : 21引言ak = 0 i ( t) cos ( kt) dt

2、TT 2Tbk = T 0 i ( t) sin ( kt) dt任何一種保護功能的實現(xiàn)都必須要有相應的算法來支持 . 我們對算法性能優(yōu)劣的評價也取決于該算法能否在較 短數(shù)據(jù)窗中從交流信號的若干采樣值中獲得基波分量或者 某次諧波分量較為精確的估計值 . 傅氏算法是當今較常用 的算法 ,該算法具有很強的濾除諧波分量的能力 ,但缺點是 其本身不能濾除衰減的直流分量 . 為了克服信號中衰減的非周期分量的影響 , 國內(nèi)外很多學者作了大量的分析研究工作 ,并提出了一些相應的保護算法 . 本文利用 MATLAB 對 傳統(tǒng)全波傅氏算法和 2 種改進后的傅氏算法進行仿真 ,并對 仿真結(jié)果進行比較和探討.經(jīng)過

3、A/ D ( 模擬量/ 數(shù)字量) 變換采樣量化后 , 連續(xù)量變?yōu)殡x散量 , 連續(xù)量求積變?yōu)殡x散量求和 , 從而有 : 2N2ak = in cos ( nk ) 式 (1)N n = 1N 2N2bk = N in sin ( nk ) 式 (2)Nn = 1上式中 , N 為一個周期 T 內(nèi)的采樣點數(shù) ; k 是諧波次數(shù) ;n 為離散的采樣點 .2. 2傳統(tǒng)傅氏算法的仿真在用 MATLAB 進行仿真模擬時 ,要注意 N 數(shù)目的選擇 ,這將對仿真效果產(chǎn)生影響. 仿真程序的流程圖如圖 1.2 傳統(tǒng)傅氏算法及其仿真2. 1傳統(tǒng)傅氏算法的推導我們以電流為例 , 設(shè)定故障電流波形為以下形式 :ni (

4、 t) = I0 e - at + Ik sin ( kt + k )k = 1其中 k 為諧波次數(shù) ,t 為基波的角頻率 . 設(shè) IRk = Ik sin(k ) , Iik = Ik cos (k )則原式轉(zhuǎn)化為 :ni ( t) = I0 e - at + IRk cos ( kt) + I1 k sin ( kt) k = 1仿真程序的部分代碼如下 :%FFTglobal step %運算步數(shù)global bitbutterfly %每個蝴蝶結(jié)中所包含的點數(shù)收稿日期 :2005 - 05 - 25作者簡介 :何志勤 (1982 - ) ,男 ,江西省九江市人 ,華東交通大學 03 級碩

5、士研究生 ,研究方向為電力系統(tǒng)及其自動化 .global frequency %每步運算中蝴蝶結(jié)的個數(shù)for step = 1 :5for bitbutterfly = 1 :2( step - 1)i = (2 (5 - step) ) 3 ( bitbutterfly - 1) ; RW = cos (2 3 pi 3i/ 32) ; % W 的實部IW = ( - 1) 3 sin (2 3 pi 3 i/ 32) ; % W 的虛部for frequency = 1 :2(5 - step)temp = (frequency - 1) 3 (2step) + bitbutterfly

6、; TR = dataR(temp) ; TI = data I(temp) ;dataR (temp) = dataR (temp) + RW 3 dataR (temp + 2 ( step -1) ) - IW 3 data I(temp + 2( step - 1) ) ;data I(temp) = data I ( temp ) + RW 3 data I ( temp + 2 ( step -1) ) + IW 3 dataR (temp + 2( step - 1) ) ;temp1 = dataR(temp + 2( step - 1) ) ;dataR(temp + 2(

7、step - 1) ) = TR - ( RW 3 dataR ( temp + 2 ( step - 1) ) - IW 3 data I(temp + 2( step - 1) ) ) ;data I(temp + 2 ( step - 1) ) = TI - ( RW 3 data I ( temp + 2 ( step - 1) ) + IW 3 temp1) ;end end end筆者選擇 N 為 32 點 , 并選取兩組故障電流進行仿真 ,從而比較非周期的衰減直流分量對運行結(jié)果的影響 . 最后 觀察的是每次諧波的幅值和真實值的差距 , 根據(jù)式 ( 1) 和式仿真程序的部分代碼如下

8、 :p = 3 ; sum = 0 ;for j = 1 :32 %計算 I (0)sum = x (j) + sum ;endy = sum/ 32 ; d = x(1 + p) - x (33 + p) ;z = y/ d ; %計算 Kfor k = 1 :32dataR ( k) = 0 ; data I( k) = 0 ;for n = 1 + p :32 + pdatar ( n) = x ( n) 3 cos (2 3 pi 3 n 3 k/ 32) ; datai ( n) = x ( n)3 ( - 1) 3 sin (2 3 pi 3 n 3 k/ 32) ;datarr

9、( n) = 2 3 y/ (1 + 4 3 pi 3 pi 3 k 3 k 3 (z2) ) ; %計算 實部誤差dataii ( n) = k 3 2 3 pi 3 z 3 datarr ( k) ; %計算虛部誤差dataR( k) = dataR ( k) + datar ( n) - datarr ( n) ; data I ( k) =data I( k) + datai ( n) - dataii ( n) ;end end程序中由 p 來控制起始 采 樣 點 的 位 置 , 通 過 對 p 的 調(diào)整 ,來測試其對運行結(jié)果的影響 . 結(jié)果說明 p = 3 時的仿真效 果比較好 .

10、3. 3濾除衰減直流分量的改進全波算法 2我們對輸入信號進行等間隔采樣 . 選取三個數(shù)據(jù)窗 ,時 間間隔為一個采樣周期 . 可以得到 :a2 + b2 .(2) , 幅值即為 :| X ( k) | =k k3 兩種改進后的傅氏算法及其仿真 .3. 1傳統(tǒng)傅氏算法的誤差來源傳統(tǒng)傅氏算法的誤差主要來源于兩方面 :1) 用離散值累 加代替連續(xù)積分 ,其結(jié)果受到采樣頻率的影響 . 此外計算要用到全部 N 個采樣值 ,因此 ,計算須在故障后第 N 個采樣值 出現(xiàn)時 ,才正確 . 在此之前 ,N 個采樣值中有一部分吻故障前的數(shù)值 ,從而使結(jié)果不精確 . 2) 傳統(tǒng)傅氏算法無法濾掉衰減的直流分量 .3.

11、2濾除衰減直流分量的改進全波算法 1所謂改進就是在全波傅氏變換提取出基波或各次諧波 分量的基礎(chǔ)上 , 減去直流衰減分量帶來的誤差 . 我們將實部a= I cos ( + k t) + ak = Ikcosk + abk = Ik sink + bkkka;ak = IRk + a 式 (3)bk = Iik + b 式 (4)b= I sin ( + kt) + k kkb和虛部的誤差設(shè)為 a 和b , 則有 :; 設(shè) ak = Ik cos (k + 2 k t) + a;令 : A = ab -k a + k b ; B = i (0) - i ( N) , 也就是說 , 在傳統(tǒng)傅氏算法的

12、數(shù)據(jù)窗 32 個點的基礎(chǔ)上再加一個點 , 取第一個點減去第 N + 1 個點 . 通過 計算 , 可得校正量 :a = 2/ ( 1 + 42 k2 K2 ) I 3 ( 0) ;b = 2 ka ;c k s kbk = Ik sin (k + 2 kt) + bbk -kc bk + ks ak ; C = ak -kc ak + ks bk ; D = bk -kc bk + ks ak ;A ( k T - kc ) - B ks| C| + | D|其中 K = I 3 (0) / . 在交流采樣算法中 , I 3 (0) = 1 i ( n) ,N - 1+ | B | ;a =;

13、b進一步得出 : k T = | A|1 + k2T - 2 kc kTN n = 0即采樣點數(shù)值之和與總采樣點數(shù)之商 . 將計算出來的 a 和b 回帶入式 (3) 和式 (4) , 計算 ak 和 bk . 仿真程序的流程圖如圖 2 :B ( k T - kc ) - Aks= ;21 + k T - 2 kc kT 22其中 : ks = sin ( N k) ; kc = cos ( N k) ; 將計算出來的 a 和 b回帶入式 (3) 和式 (4) ,計算 ak 和 bk . 仿真程序的流程圖如圖 3.仿真程序的部分代碼如下 :for n = (2 + p) : (33 + p) %

14、 取數(shù)據(jù)窗 3datar3 ( n) = x ( n) 3 cos (2 3 pi 3 n 3 k/ 32) ; datai3 ( n) = x( n) 3 ( - 1) 3 sin (2 3 pi 3 n 3 k/ 32) ;data R3 ( k) = data R3 ( k ) + datar3 ( n) ; data I3 ( k ) = data I3( k) + data i3 ( n) ;endbb = data I1 ( k) - data I2 ( k) ; % 計算 Bcc = data R2 ( k) - data R3 ( k) ; % 計算 C dd = data I

15、2 ( k) - data I3 ( k) ; % 計算 D kt = abs ( cc) / ( abs ( kc 3 aa + ks 3 bb) ) ;kc1 = ( bb 3 ( kc 3 kt - 1) + aa 3 ks 3 kt ) / ( 1 + kt2 - 2 3kt 3 kc) ; %虛部誤差ks1 = ( aa 3 ( kc 3 kt - 1) -kt 3 kc) ; %實部誤差bb 3 ks 3 kt ) / ( 1 + kt2 - 2 3data R ( k) = data R1 ( k ) + ks1 ; data I ( k ) = data I1 ( k ) +a

16、a = data R1 ( k) - data R2 ( k) ;% 計算 Akc1 ;end筆者在編程時 , 對誤差計算做了修改 , 結(jié)果更佳 . 要注意 誤差和初值的加減問題 , 在算法 2 的程序中誤差是和初值相 加 , 在算法 1 中是相減 , 否則結(jié)果誤差會反而增大 , 此處不再 做深入分析 .4傳統(tǒng)算法和改進算法的性能比較 ( N 統(tǒng)一取 32)4. 1含有較大的衰減直流分量根據(jù)程序 , 我們?nèi)?:i ( t) = 30 e - 40 t + 30sin (100t +/ 3) + 6sin (200t +/ 4) +15sin (300t +/ 6) + 5sin (400t +

17、/ 3) + 10sin (500t +/ 4) ;仿真結(jié)果如表 1.表 1含較大衰減值流分量仿真結(jié)果表基波二次諧波三次諧波算法類型幅值誤差幅值誤差幅值誤差傳統(tǒng)傅氏算法 ( FFT)改進算法 1改進算法 233 . 8 21530 . 4 49530 . 0 00012. 74 %1 . 50 %0 . 00 %8. 4 2006. 1 9096. 0 00040 . 33 %3 . 18 %0 . 00 %16 . 7 70814 . 4 79615 . 0 00011 . 81 %- 3 . 47 %0 . 00 %計算后的仿真數(shù)值波形如圖 4 :改進算法 2 仿真波形4. 2含有較小的衰

18、減直流分量根據(jù)程序 , 我們?nèi)?:i ( t) = 8 e - 20 t + 30sin (100t +/ 3) + 6sin (200t +/ 4) +15sin (300t +/ 6) + 5sin (400t +/ 3) + 10sin (500t +/ 4) ;仿真結(jié)果如下表 2 :表 2 含較小衰減值流分量仿真結(jié)果表基波二次諧波三次諧波算法類型幅值誤差幅值誤差幅值誤差傳統(tǒng)傅氏算法 ( FFT)改進算法 1改進算法 230 . 5 41330 . 0 17030 . 0 0001. 8 %0 . 06 %0 . 00 %6. 3 6455. 9 3876. 0 0006 . 08 %-

19、 1 . 02 %0 . 00 %15 . 2 79614 . 8 94715 . 0 0001 . 86 %- 0 . 70 %0 . 00 %計算后的仿真數(shù)值波形如圖 5.5 全波傅氏算法小結(jié)全波富氏算法能濾除所有的整次諧波 分 量 , 且 穩(wěn) 定 性好 , 但其數(shù)據(jù)窗較長 , 所以其響應速度較慢 . 通過上面的仿真 結(jié)果 , 我們可以得出以下結(jié)論 :1) 對照同一故障波形進行分析 , 可以看出 :2 種改進算 法所算結(jié)果都比較傳統(tǒng)傅氏算法要精確 .2) 在故障波形諧波含量都相同的情況下 , 衰減直流分 量的初始值越小 , 衰減時間常數(shù)越大時 , 衰減非周期分量對基波以及各次諧波的影響越小

20、 .3) 筆者將采樣點數(shù) N 從 16 提高到 32 后 , 3 種算法誤差 均有一定程度的減小 . 如果進一步提高采樣率到 64 點/ 周或128 點/ 周時 , 對提高計算精度作用不大 , 而計算時間相應增 加 . 因此 , 將采樣點定為 32 點/ 周對減少誤差有益 .雖然改進算法總體要比傳統(tǒng)算法精確 , 但都還存在一個 無法避免的計算誤差 , 那就是所取數(shù)據(jù)窗內(nèi)的正常數(shù)據(jù)和故 障數(shù)據(jù)的分辨問題 . 為了全部使用故障后的采樣值 , 往往取i N/ 2 , i 表示從故障起始時刻開始第 i 個采樣點 . 本文中 是通過對改進算法中 p 的調(diào)整 , 模擬故障發(fā)生后 , 起始采樣 點的位置 . 理論上如果有充足的判據(jù)是可以的 , 但無疑增加了保護計算的時延和難度 , 而且跳過的半個周期采樣值不予計算 , 這肯定會給故障時刻電流的計算帶來一定的誤差 . 因 此 , 這個問題還有待進一步的研究 .參考文獻 :1劉建剛 ,孫同景. 基于 FFT 的傅立葉算法在微機繼電保護中的應用J . 繼電器. 2004 ,10 :24 - 26 .李孟秋 ,王耀南 ,王輝. 基于全周波富氏算法濾除衰減直流分量 新方法J . 湖南大學學報. 2001 ,1 :59 - 63 .周大敏. 一種消除非周期分量對非遞推富氏算法影響的精確方法J . 繼電器. 1998 ,4