傅里叶变换FFT算法的介绍及其在微机继电保护中的应用 陆志强

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 为