FFT(离散傅氏变换的快速算法)

1、FFT（离散傅氏变换的快速算法）目录1算法简介2DFT算法3源码表示4MATLAB中FFT的使用方法1算法简介编辑FFT（Fast Fourier Transformation），即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的FFT算法图(Bufferfly算法)发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一

2、次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出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次

3、，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。2DFT算法编辑For length N input vector x, the DFT is a length N vector X,with elementsNX(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 = k = N.n=1The inverse DFT (compu

4、ted by IFFT) is given byNx(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 = n 0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数/ pr: l=0时，存放N点采样数据的实部/ l=1时, 存放傅立叶变换的N个实部/ pi: l=0时，存放N点采样数据的虚部/ l=1时, 存放傅立叶变换的N个虚部/ 出口参数：/ fr: l=0, 返回傅立叶变换的实部/ l=1, 返回逆傅立叶变换的实部/ fi: l=0, 返回傅立叶变换的虚部/ l=1, 返回逆傅立叶变换的虚部/ pr: il=1,i=0 时，返回傅立

5、叶变换的模/ il=1,i=1 时，返回逆傅立叶变换的模/ pi: il=1,i=0 时，返回傅立叶变换的辐角/ il=1,i=1 时，返回逆傅立叶变换的辐角void fft(double pr, double pi, int n, int k, double fr, double fi, int l, int il)int it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;for(it=0;it=n-1;m=it+)is=0;for(i=0;i=k-1;i+)j=m/2;is=2*is+(m-2*j);m=j;frit=pris;fiit=pi

6、is;/-pr0=1.0;pi0=0.0;p=6.283185306/n;pr1=cos(p);pi1=-sin(p);if (l)pi1=-pi1;for(i=2;i=n-1;i+)p=pri-1*pr1;q=pii-1*pi1;s=(pri-1+pii-1)*(pr1+pi1);pr=p-q;pi=s-p-q;for(it=0;it=0;l0-)m/=2;nv=1;for(it=0;it=(m-1)*nv;it+=nv)for(j=0;j=(nv/2)-1;j+)p=prm*j*frit+j+nv/2;q=pim*j*fiit+j+nv/2;s=prm*j+pim*j;s*=(frit+j

7、+nv/2+fiit+j+nv/2);poddr=p-q;poddi=s-p-q;frit+j+nv/2=frit+j-poddr;fiit+j+nv/2=fiit+j-poddi;frit+j+=poddr;fiit+j+=poddi;if(l)for(i=0;i=n-1;fr/=n,fii+/=n);if(il)for(i=0;i=n-1;i+)pr=sqrt(fr*fr+fi*fi);if(fabs(fr)0?90.0-90.0;elsepi=atan(fi/fr)*360.0/6.283185306;return;源码（2）ps:可以运行的/ The following line mu

8、st be defined before including math.h to correctly define M_PI#define _USE_MATH_DEFINES#include #include #include #define PI M_PI /* pi to machine precision, defined in math.h */#define TWOPI (2.0*PI)/*FFT/IFFT routine. (see pages 507-508 of Numerical Recipes in C)Inputs:data : array of complex* dat

9、a points of size 2*NFFT+1.data0 is unused,* the nth complex number x(n), for 0 = n = length(x)-1, is stored as:data2*n+1 = real(x(n)data2*n+2 = imag(x(n)if length(Nx) = length(x).isign: if set to 1,computes the forward FFTif set to -1,computes Inverse FFT - in this case the output values haveto be m

10、anually normalized by multiplying with 1/NFFT.Outputs:data : The FFT or IFFT results are stored in data, overwriting the input.*/void four1(double data, int nn, int isign)int n, mmax, m, j, istep, i;double wtemp, wr, wpr, wpi, wi, theta;double tempr, tempi;n = nn 1;j = 1;for (i = 1; i i) tempr = dat

11、aj; dataj = datai; datai = tempr;tempr = dataj+1; dataj+1 = datai+1; datai+1 = tempr;m = n 1;while (m = 2 & j m) j -= m;m = 1;j += m;mmax = 2;while (n mmax) istep = 2*mmax;theta = TWOPI/(isign*mmax);wtemp = sin(0.5*theta);wpr = -2.0*wtemp*wtemp;wpi = sin(theta);wr = 1.0;wi = 0.0;for (m = 1; m mmax;

12、m += 2) for (i = m; i = n; i += istep) j =i + mmax;tempr = wr*dataj - wi*dataj+1;tempi = wr*dataj+1 + wi*dataj;dataj = datai - tempr;dataj+1 = datai+1 - tempi;datai += tempr;datai+1 += tempi;wr = (wtemp = wr)*wpr - wi*wpi + wr;wi = wi*wpr + wtemp*wpi + wi;mmax = istep;在C+环境下的源码bool FFT(complex * TD,

13、 complex * FD, int r)/一维快速Fourier变换。/complex * TD 指向时域数组的指针; complex * FD 指向频域数组的指针; r 2的幂数，即迭代次数LONG count; / Fourier变换点数int i,j,k; / 循环变量int bfsize,p; / 中间变量double angle; / 角度complex *W,*X1,*X2,*X;count = 1 r; / 计算Fourier变换点数为1左移r位W = new complexcount / 2;X1 = new complexcount;X2 = new complexcoun

14、t; / 分配运算所需存储器/ 计算加权系数（旋转因子w的i次幂表）for(i = 0; i count / 2; i+)angle = -i * PI * 2 / count;W i = complex (cos(angle), sin(angle);/ 将时域点写入X1memcpy(X1, TD, sizeof(complex) * count);/ 采用蝶形算法进行快速Fourier变换for(k = 0; k r; k+)for(j = 0; j 1 k; j+)bfsize = 1 (r-k);for(i = 0; i bfsize / 2; i+)p = j * bfsize;X2

15、i + p = X1i + p + X1i + p + bfsize / 2 * Wi * (1k);X2i + p + bfsize / 2 = X1i + p - X1i + p + bfsize / 2 * Wi * (1k);X = X1;X1 = X2;X2 = X;/ 重新排序for(j = 0; j count; j+)p = 0;for(i = 0; i r; i+)if (j&(1i)p+=1(r-i-1);FDj=X1p;/ 释放内存delete W;delete X1;delete X2;return true;在Matlab环境下的源码function X=myfft(

16、x)%myfft函数 用递归实现N=length(x);t=log2(N);t1=floor(t);t2=ceil(t);if t1=t2; %若x的长度N不为2的整数次幂，则补0至最接近的2的整数次幂x=x zeros(1,2t2-N);N=2t2;endw0=exp(-j*2*pi/N);X=zeros(1,N);if N=2X(1)=x(1)+x(2);X(2)=x(1)-x(2);elsen=1:N/2;xe(n)=x(2*n-1);xo(n)=x(2*n);XE=myfft(xe); %递归调用XO=myfft(xo);for n=1:N/2X(n)=XE(n)+XO(n)*(w0(

17、n-1);X(n+N/2)=XE(n)-XO(n)*(w0(n-1);endend4MATLAB中FFT的使用方法编辑2一.调用方法X=FFT(x)；X=FFT(x，N)；x=IFFT(X);x=IFFT(X,N)用MATLAB进行谱分析时注意：（1）函数FFT返回值的数据结构具有对称性。例：N=8;n=0:N-1;xn=4 3 2 6 7 8 9 0;Xk=fft(xn)Xk =39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.292

18、9iXk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。（2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。二.FFT应用举例例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。clf;fs=100;N=128; %采样频率和数据点数n=0:N-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=f

19、ft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求得Fourier变换后的振幅f=n*fs/N; %频率序列subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=128);grid on;subplot(2,2,2),plot(f(1:N/2),mag(1:N/2); %绘出Nyquist频率之前随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=128);grid on;%对信号采样数据为1024点的处理fs=100;N=1024;n=0:N-

20、1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求取Fourier变换的振幅f=n*fs/N;subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=1024);grid on;subplot(2,2,4)plot(f(1:N/2),mag(1:N/2); %绘出Nyquist频率之前随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=102

21、4);grid on;运行结果：fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率01进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。

22、例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制：（1）数据个数N=32，FFT所用的采样点数NFFT=32；（2）N=32，NFFT=128；（3）N=136，NFFT=128；（4）N=136，NFFT=512。clf;fs=100; %采样频率Ndata=32; %数据长度N=32; %FFT的数据长度n=0:Ndata-1;t=n/fs; %数据对应的时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %时间域信号y=fft(x,N); %信号的Fourier变换mag=abs(y); %求取振幅

23、f=(0:N-1)*fs/N; %真实频率subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel(频率/Hz);ylabel(振幅);title(Ndata=32 Nfft=32);grid on;Ndata=32; %数据个数N=128; %FFT采用的数据长度n=0:Ndata-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:N-1)*fs/N; %真实频率subplot(2,2,2),plot(

24、f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel(频率/Hz);ylabel(振幅);title(Ndata=32 Nfft=128);grid on;Ndata=136; %数据个数N=128; %FFT采用的数据个数n=0:Ndata-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:N-1)*fs/N; %真实频率subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel(频率/Hz);ylabel(振幅);title(Ndata=136 Nfft=128