1、傅里葉變換FFT算法的介紹及其在微機繼電保護中的應用陸志強(溫州大學物理與電子信息工程學院，11電信)摘要：傳統的微機繼電保護算法中 ,一般使用梯形算法來計算周期信號的直流分量和各次諧波的系數 ,此方法計算比較復雜 。本文提出了一種基于 FFT 的算法 。該算法利用 FFT 可以由輸入序列直接計算出輸入信號的直流分量和各次諧波的幅值和相角的特點 ,大大簡化了諧波分析的計算 。與梯形算法相比 ,該算法具有精度高 、計算量小 、更易在數字信號處理器上實現等優點 。因而可以取代梯形算法來計算諧波系數 。針對 FFT計算 ,還介紹了正弦信號采樣頻率的選擇方法 。關鍵字： 傅里葉算法 

2、; FFT; 諧波分析；微機繼電保護。The Introduction of Fourier algorithm based on FFT in Modif ied model of power metering Lu Zhiqiang(School of Physics and Electronic and Information Engineering, Wenzhou University,)Abstract: In microcomputer relay protection of traditional algorithm, coefficient of DC c

3、omponent generally use the trapezoidal algorithm to calculate the periodic signal and harmonic,and this method is very complex. This paper presents an algorithm based on FFT. The algorithm makes use of the FFT and it can be calculated directly from the input sequence characteristics of amplitude and

4、 phase of the DC component of the input signal and harmonic, greatly simplifies the calculation of harmonic analysis. Compared with the trapezoidal algorithm, this algorithm has high precision, small computation, easily realized in digital signal processor. So that you can replace trapezoidal algori

5、thm to calculate the harmonic coefficient. For the FFT calculation, the selection method of sine signal sampling frequency is also presented.Keywords: Fourier algorithm；FFT；harmonic analysis；Modif ied model of power metering.一、傅立葉變換FFT算法簡介：計算離散傅里葉變換的一種快速算法，簡稱FFT。快速傅里葉變換是1965年由J.W.庫利和T.W.圖基提出的。采用這種算法

6、能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數N越多，FFT算法計算量的節省就越顯著。有限長序列可以通過離散傅里葉變換(DFT）將其頻域也離散化快速傅里葉變換成有限長序列。但其計算量太大，很難實時地處理問題，因此引出了快速傅里葉變換(FFT). 1965年，Cooley和Tukey提出了計算離散傅里葉變換（DFT）的快速算法，將DFT的運算量減少了幾個數量級。從此，對快速傅里葉變換（FFT）算法的研究便不斷深入，數字信號處理這門新興學科也隨FFT的出現和發展而迅速發展。根據對序列分解與選取方法的不同而產生了FFT的多種算法，基本算法是基2DIT和基2DI

7、F。FFT在離散傅里葉反變換、線性卷積和線性相關等方面也有重要應用。 快速傅氏變換（FFT），是離散傅氏變換的快速算法，它是根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。它對傅氏變換的理論并沒有新的發現，但是對于在計算機系統或者說數字系統中應用離散傅立葉變換，可以說是進了一大步。 設快速傅里葉變換x(n）為N項的復數序列，由DFT變換，任一X（m）的計算都需要N次復數乘法和N-1次復數加法，而一次復數乘法等于四次實數乘法和兩次實數加法，一次復數加法等于兩次實  快速傅里葉變換數加法，即使把一次復數乘法和一次復數加法定義成一次“運算”（四次實數

8、乘法和四次實數加法），那么求出N項復數序列的X（m），即N點DFT變換大約就需要N2次運算。當N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列（設N=2k,k為正整數），分為兩個N/2項的子序列，每個N/2點DFT變換需要（N/2）2次運算，再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以后，總的運算次數就變成N+2*(N/2）2=N+N2/2。繼續上面的例子，N=1024時，總的運算次數就變成了525312次，節省了大約50%的運算量。而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組

9、的DFT運算單元，那么N點的DFT變換就只需要Nlog2N次的運算，N在1024點時，運算量僅有10240次，是先前的直接算法的1%，點數越多，運算量的節約就越大，這就是FFT的優越性。 二、微機繼電保護簡介：微機繼電保護的輸入信號是電力系統的模擬量，而計算機只能對數字量進行計算和判斷，因此由電力系統經電壓互感器和/或電流互感器輸入的模擬量必先經過預處理繼電保護在大部分情況下取用輸入信號中的基波模擬量。根據采樣定理，如被測信號頻率（或要求保留的最高次諧波頻率）為,則采樣頻率必須大于2,kg1否則由采樣值不可能擬合還原成原來的曲線。對于那些大于0.5頻率的諧波分量,必須在進入采樣器之前，利用模擬

10、式低通濾波器（前置模擬濾波）將其濾掉。 由于輸入信號常常有多個，故設置多路轉換器將輸入模擬信號逐個交與A/D變換器轉化成數字量這些數字量應在存儲器中按先后順序排列，以便后續功能處理程序取用。2.1 濾波：為了保證計算機計算和判斷的正確，實現以某種頻率的正弦電量為基礎的繼電保護原理,必須將經A/D變換后的數字量再經一次濾波。由于數字濾波器精度高、可靠而且調整靈活，通過時分復用可使裝置簡化，因此微機保護中普遍采用數字濾波器。數字濾波器本身可理解為一個計算程序或算法，它將代表輸入信號的數字時間序列轉換為代表輸出信號的數字時間序列，使信號按照預定的形式變化。微機繼電保護中應根據電力系統信號的特點和保護

11、原理的要求設計、選擇相應的數字濾波器。數學濾波器的主要性能指標是頻域特性、時延和計算量。2.2 算法問題：對離散和量化的數字式采樣序列，用數學運算方法實現故障量的測量，這就是微機保護的算法問題。要求運算精度滿足保護的實際需要，同時計算時間又盡可能短。微機繼電保護的研究初期，一些算法是基于被采樣的電壓、電流均系純正弦波的，為此應將輸入信號進行預處理。稍后，相繼提出傅里葉算法和沃爾什函數算法。它們假定輸入信號中含有非周期分量、基波和高次諧波。這些算法本身具有很強的濾去高次諧波的功能，因此無需另設數字濾波器，但對非周期分量必須采取其他措施。由于電力系統中大量應用鐵磁非線性元件，輸電線路分布電容和串聯

12、、并聯電容，以及電壓互感器、電流互感器的暫態特性等因素的影響，使微機繼電保護輸入信號中還含有許多隨機高頻分量，它們起著干擾或噪聲的作用。對此，可采用最小二乘曲線擬合算法或對計算結果采取平滑措施。上述種種算法都是先算出電壓、電流的大小和相位，然后根據保護的動作判據作進一步的運算，最終實現其保護功能。也有一些算法將電量運算與保護動作判據運算直接結合在一起，例如用離散值直接實現的方向阻抗繼電器的算法。2.3 優點與應用前景：由于計算機的優越存儲能力，可以方便地得到保護需要的故障分量并準確地予以保持，這是模擬式保護裝置難以達到的。由于計算機的強大運算能力，可以實現一些以往模擬式保護裝置無法實現的復雜保

13、護動作特性、自適應性的定值或特性改變以及良好的自檢功能。同常規繼電保護相比，微機繼電保護的抗電磁干擾能力較弱，因此，它的廣泛應用受到一定的限制。應用微機繼電保護時，應特別注意解決好電磁兼容性問題。2.4 系統采集電量變換：微機保護中通常要求輸入信號為±5V或±10V的電壓信號， 這是由所采用的模數轉換器所決定的。而從被保護的電力線路或電氣設備的電流互感器、電壓互感器或其它變換器上取得的二次數值對微機電路是不適用的，所以需要進行電量變換。電量變換一般采用中間變換器來實現。2.5 采樣定理和模擬低通濾波：由于輸入信號是模擬量，因此信號在進入微型計算機之前首先進行采樣并

14、保持。采樣就是把一個是時間連續函數信號變換為對時間。2.6 微機繼電保護的特點：1）、 微機繼電保護集測量、控制、監視、保護、通信等多種功能于一體的電力自動化高新技術產品,是構成智能化開關柜的理想電器單元。 2）、多種功能的高度集成，靈活的配置，友好的人機界面，使得該通用型微機綜合保護裝置可作為35KV及以下電壓等級的不接地系統、小電阻接地系統、消弧線圈接地系統、直接接地系統的各類各類電器設備和線路的保護及測控，也可作為部分66KV、110KV電壓等級中系統的電壓電流的保護及測控。 3）、采用32位數字處理器（DPS）具有先進的內核結構，高速運算能力和實時信號處理

15、等優點。 4）、支持常規的RS485總線和及CAN（DEVICENET）現場總線通訊，CAN總線具有也錯帖自動重發和故障節點自動脫離等糾錯機制，保護信息的實施性和可靠性。 5）、完善的自檢能力，發現裝置異常自動報警；具有自保護能力，有效防止接線錯誤和非正常運行引起的裝置永久性損壞；免維護設計，無需在現場調整采樣精度，測量精度不會因為環境改變和長期運行引起誤差增大。三、FFT算法的基本應用：3.1 利用FFT計算連續時間信號的傅里葉變換：設是連續時間信號，并假設時，則其傅里葉變換由下式給出令是一個固定的正實數，是一個固定的正整數。當時，利用FFT算法可計算。已知一個固定的時間

16、間隔，選擇足夠小，使得每一個秒的間隔內，的變化很小，則式中積分可近似為 （1）假設足夠大，對于所有的整數，幅值很小，則式（1）變為 （2）當時，式（28）兩邊的值為 （3）其中代表抽樣信號的點。最后令，則上式變為 （4）首先用FFT算法求出，然后可用上式求出時的。應該強調的是，式（2）只是一個近似表示，計算得到的只是一個近似值。通過取更小的抽樣間隔，或者增加點數，可以得到更精確的值。如果時，幅度譜很小，對應于奈奎斯特抽樣頻率，抽樣間隔選擇比較合適。如果已知信號只在時間區間內存在，可以通過對時的抽樣信號補零，使足夠大。例1 利用FFT計算傅里葉變換如圖1所示的信號其傅里葉變換為：利用下面的命令，

17、可得到的近似值和準確值。 圖1 連續時間信號x(t) N=input('Input N:');T=input('Input T:');%計算X(w)近似值t=0:T:2;x=t-1 zeros(1,N-length(t);X=fft(x);gamma=2*pi/(N*T);k=0:10/gamma;Xapp=(1-exp(-i*k*gamma*T)/(i*k*gamma)*X;%計算真實值X(w)w=0.05:0.05:10;Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w)./(w.*w);plot(k*gamma,abs(Xapp(1

18、:length(k),'o',w,abs(Xact);legend('近似值','真實值');xlabel('頻率(rad/s)');ylabel('|X|')運行程序后輸入N=128，T=0.1，此時，得到實際的和近似的傅里葉變換的幅度譜如圖2所示，此時近似值已經相當準確。通過增加NT可以增加更多的細節，減少T使得到的值更精確。再次運行程序后輸入N=512，T=0.05，此時，得到實際的和近似的傅里葉變換的幅度譜如圖3所示。圖2 N=128，T=0.1時的幅度譜圖3 N=512，T=0.05時的幅度譜3.2 利用

19、FFT計算離散信號的線性卷積：已知兩個離散時間信號與，取，對和右端補零，使得 （5）利用FFT算法可以求得和的L點DFT，分別是和，利用DTFT卷積性質，卷積等于乘積的L點DFT反變換，這也可以通過FFT 算法得到。例2 利用FFT計算線性卷積：已知，其中為單位階躍序列，信號如圖4所示。由于當時，很小，故可以取為17；N取10，。利用下面的Matlab命令，可得到、的卷積圖形如圖4所示。subplot(3,1,1);n=0:16;x=0.8.n;stem(n,x);xlabel('n');ylabel('xn');subplot(3,1,2);n=0:

20、15;y=ones(1,10) zeros(1,6);stem(n,y);xlabel('n');ylabel('yn')subplot(3,1,3);L=26;n=0:L-1;X=fft(x,L);Y=fft(y,L);Z=X.*Y;z=ifft(Z,L);stem(n,z);xlabel('n');ylabel('zn')圖4 信號xn、yn及其卷積zn=xn*yn利用下面的Matlab命令，可得到信號xn、yn的幅度譜與相位譜如圖5所示。subplot(2,2,1);L=26;k=0:L-1;n=0:16;x=0.8.n;X

21、=fft(x,L);stem(k,abs(X);axis(0 25 0 5);xlabel('k');ylabel('|Xk|')subplot(2,2,2);stem(k,angle(X);axis(0 25 -1 1);xlabel('k');ylabel('Angle(Xk)(弧度)')subplot(2,2,3);y=ones(1,10);Y=fft(y,L);stem(k,abs(Y);axis(0 25 0 10);xlabel('k');ylabel('|Yk|')subplot(2,

22、2,4);stem(k,angle(Y);axis(0 25 -3 3);xlabel('k');ylabel('Angle(Yk)(弧度)')圖5 信號xn、yn的幅度譜與相位譜3.3 利用FFT進行離散信號壓縮：利用FFT算法對離散信號進行壓縮的步驟如下：1）通過采樣將信號離散化；2）對離散化信號進行傅里葉變換；3）對變換后的系數進行處理，將絕對值小于某一閾值的系數置為0，保留剩余的系數；4）利用IFFT算法對處理后的信號進行逆傅里葉變換。例3 對單位區間上的下列連續信號以采樣頻率進行采樣，將其離散化為個采樣值.用FFT分解信號，對信號進行小波壓縮，然后重構

23、信號。令絕對值最小的80%系數為0，得到重構信號圖形如圖6 a)所示，均方差為0.0429，相對誤差為0.0449；令絕對值最小的90%系數為0，得到重構信號圖形如圖6 b)所示，均方差為0.0610，相對誤差為0.0638。 a) 絕對值最小的80%系數為0的重構信號（FFT） b) 絕對值最小的90%系數為0的重構信號（FFT）圖6 用FFT壓縮后的重構信號相關Matlab程序如下：function wc=compress(w,r)%壓縮函數compress.m%輸入信號數據w,壓縮率r%輸出壓縮后的信號數據if(r<0)|(r>1) error('r 應該介于0和1之

24、間!');end;N=length(w);Nr=floor(N*r);ww=sort(abs(w);tol=abs(ww(Nr+1);wc=(abs(w)>=tol).*w;function unbiased_variance,error=fftcomp(t,y,r)%利用FFT做離散信號壓縮%輸入時間t,原信號y,以及壓縮率r%輸出原信號和壓縮后重構信號的圖像,以及重構均方差和相對l2誤差if(r<0)|(r>1) error('r 應該介于0和1之間!');end;fy=fft(y);fyc=compress(fy,r); %調用壓縮函數compr

25、ess.myc=ifft(fyc);plot(t,y,'r',t,yc,'b');legend('原信號','重構信號');unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);輸入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=fftcomp(t,f,0.8)unbiased_variance = 0.0429error =0.044

26、9如果用Harr尺度函數和Harr小波分解信號，對信號進行小波壓縮，然后重構信號。令絕對值最小的80%系數為0，得到重構信號圖形如圖7 a)所示，均方差為0.0584，相對誤差為0.0611；令絕對值最小的90%系數為0，得到重構信號圖形如圖7 b)所示，均方差為0.1136，相對誤差為0.1190。 a) 絕對值最小的80%系數為0的重構信號（Harr） b) 絕對值最小的90%系數為0的重構信號（Harr）圖7 用Harr小波壓縮后的重構信號相關Matlab程序如下：function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechie

27、s系列小波做離散信號壓縮%輸入時間t,原信號y,分解層數n,以及壓縮率r%輸出原信號和壓縮后重構信號的圖像,以及重構均方差和相對l2誤差if(r<0)|(r>1) error('r應該介于0和1之間!');end;c,l=wavedec(y,n,'db1');cc=compress(c,r); %調用壓縮函數compress.myc=waverec(cc,l,'db1');plot(t,y,'r',t,yc,'b');legend('原信號','重構信號');unbias

28、ed_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);輸入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=daubcomp(t,f,8,0.8)unbiased_variance = 0.0584error = 0.0611結論：在信號沒有突變、快變化或者大致上具有周期性的信號，用FFT可以處理得很好（甚至比小波還要好）。四、FFT算法在微機繼電保護中的應用：在微機繼電保護中 ,有兩種形式的濾波器可供選擇：一種

29、是模擬濾波器 ,另一種是數字濾波器 。同模擬濾波器相比 ,由于數字濾波器具有高精確性、高靈活性和高穩定性以及便于時分復用等優點 ,因此目前所研制的電力監控產品中 ,絕大多數都用到數字濾波算法 。其中傅里葉算法因能夠有效地去除直流分量和諧波干擾 ,并且可以有選擇地單獨計算諧波分量 ,所以被廣泛地應用于諧波分析中 。FFT 由于具有原位性 ,計算量較小并且易于流水線操作等特點 ,所以非常適合用數字信號處理器(DSPs) 進行處理 。我們可以通過一定的轉換和計算 ,用 FFT 來實現傅里葉算法 ,可以大大減小運算量 ,而且使其更易于通過 DSPs 實現 。傅里葉算法的基本思想源于傅里葉級數 。該算法

30、假設輸入信號為一周期性信號 ,即輸入信號中除基頻分量外 ,只包含恒定的直流分量和各種整次諧波分量 。此時電壓 ( 電流) 輸入信號可表示為：也可以合并為：其中：式中：T1 為周期信號的周期, c0為直流分量, ck 為 k次諧波的幅值, ck2為k次諧波的有效值。對于周期連續信號x ( t ) , 式 ( 2)和式 ( 3)的積分可用梯形法則 1 求得 :其中： N 為一周期內采樣的點數 ; x ( n ) 為第 n 次采樣值 , n = 0 ,1 ,2 , , N - 1 。(當輸入為電壓 ( 電流 ) 信號時 , 由式 （4）、（5）、（6）、（7） 得出的 ck 和 分別對應著電壓 (

31、電流) 的 k(k次諧波的幅值 Uk ( Ik ) ,和 k 次諧波的相位uk (ik ) ,由此可計算出電壓 ( 電流) 的 k 次諧波的有效值 。在此基礎上還可以計算出 k 次諧波的有功功率 Pk ,無功功率 Qk ,視在功率 S k 2 。同時也可以計算出 k 次諧波的電壓 ( 電流 ) 含有率 HRUk ( HRIk )同理 也 可 以 算 出 電 壓 ( 電 流 ) 諧 波 總 畸 變 率THD u ( THD i )：4.1 基于 FFT 的傅里葉算法的實現： 在傅里葉算法中 ,每計算 1 次 ak 或 bk 就要計算 1 次式 ( 6) 或 ( 7) ,很不方便 ;而且當需要計算

32、的諧波次數很高時 , 就會造成很大的計算量 。為了克服這些缺點 ,可以利用傅里葉級數和離散傅里葉變換(的關系 ,通過 FFT 代替梯形法則 ( 式 ( 6) 、7) ) 來計算a k 和 bk 。離散傅里葉變換 (DFT)和快速傅里葉變換(FFT) 實質上是同種變換,FFT只不過是利用 DFT 系數e N 的對稱性 、周期性和可約性等性質將長序列的 DFT 分解若干個短序列的 DFT 計算 ,然后再按一定規則將其合并 ,從而得到整個的 DFT。因此對 FFT 的研究 ,實際上就是對DFT的研究 。根據離散傅里葉變換有：其中 : N 表示時域中一周期的采樣點數 , n = 0 , 1 , 2 ,

33、 N - 1 ; k = 0 ,1 ,2 , , N - 1 。將 x ( t ) 表示成傅里葉級數的指數形式：其中：根據傅里葉指數形式和三角形式的關系有：根據傅里葉級數性質不難得到：要將連續的周期信號的傅里葉級數和 DFT 聯系起來 ,就需要在時域內對 x ( t ) 進行抽樣 , 抽樣間隔為 T 。一周期內的抽樣點數為 N , 則：。根據信號的時域和頻域的對稱關系 , 當信號在時域被抽樣后 ,其頻域內的頻譜以抽樣頻率 f s =1/T做周期性延拓 。一周期內的角頻率間隔為 1 =2/T1.頻率f 1 表示為 f 1 =1/T1. 則頻域內一周期的抽樣點數為 N 。帶入式 ( 15) 得：其

34、中 : n = 0 ,1 ,2 , N - 1 ; k = 0 ,1 ,2 , N - 1 。比較式 ( 8) 不難看出 :式 ( 17) 表明了連續的周期信號被抽樣后其離散傅里葉變換序列和傅里葉級數系數序列的關系 。比較式 ( 12) 和 ( 17) ,可得 :可以看出 ak 和 bk 分別與 X ( k ) 的實部和虛部相對應 ( 不是相等) 。(將式 ( 17) 帶入式 ( 13) 、14) 得 :得到了 ck 和k ,就可以按照 1 中所介紹的公式進行功率計算和諧波分析了。在現場測量中 ,要得到精確的計算結果,采樣頻率的選擇很重要。如果采樣頻率過高,雖提高了計算精度,但增加了計算量,會

35、影響到實時性;如果采樣頻率過低,會造成其頻域的混疊,而無法如實地反映出原來的信號。 對于一般的連續信號,根據時域抽樣定理,應有:fs m2f，其中: f s為采樣頻率, f m為奈奎斯特頻率,但是對于正弦信號,由于其頻譜是譜線(在±f 0處的函數) ,既不能簡單地視為帶限信號,也不能簡單地視為窄帶信號。當其初相 不確定時, 若選取f s =2f m ,有可能導致波形的嚴重失真。對于正弦信號若選擇抽樣頻率 f s = 2f 1 則會出現以下三種情況 3 :當 < =/ 2 時 ,可以由 x ( n) 重建 x ( t ) ;當 = 0 時 ,無法由 x ( n) 重建 x ( t

36、 ) ; 0 < < / 2 時 , 由 x ( n ) 重 建 出 的 不 是當x ( t ) ,而是幅值為 x ( t ) = A sin () 、初相為零的同頻余弦信號 。若 確定 , 可以得到原信號 x ( t ) ; 若 不確定 ,則無法得到原信號 x ( t ) 。只有當 f s m 時 , 才可以保證任何初相位3f情況下 ,由 x ( n) 重建 x ( t ) 。結論 顯而易見 : 若 = 0 , 則一個周期內抽得的兩點全是零 ,自然無法重建 x ( t ) ; 結論 可以通過圖 1 說明 。注 : a 是 x ( t ) = sin ( 2 f 1 + ,其中 不確定 ; b 是 f s = 2f m 時的采樣結果 ; c 是由 b 重建的信號 。圖不同初相的正弦信號的重建Fig. Rebuilt of sine signal with different phases由圖可以看到 , 由抽樣的信號 b 重建的信號c 即 x ( t ) = sin ( 2 f 1 ) 不是原信號 a ,而是幅值變cos 為