利用MATLAB编写FFT快速傅里叶变换..doc

实验目的 利用MATLAB编写FFT快速傅里叶变换。 比较编写的myfft程序运算结果与MATLAB中的FFT的有无误差。 二、实验条件 PC机，MATLAB7.0 实验原理 FFT（快速傅里叶变换）原理： 将一个N点的计算分解为两个N/2点的计算，每个N/2点的计算再进一步分解为N/4点的计算，以此类推。根据DFT的定义式，将信号x[n]根据采样号n分解为偶采样点和奇采样点。设偶采样序列为y[n]=x[2n]，奇采样序列为z[n]=x[2n+1]。 上式中的为旋转因子。下式则为y[n]与z[n]的表达式： 蝶形变换的原理： 下图给出了蝶形变换的运算流图，可由两个N/2点的FFT（Y[k]和Z[k]得出N点FFT X[k]）。同理，每个N/2点的FFT可以由两个N/4点的FFT求得。按这种方法，该过程可延迟后推到2点的FFT。 下图为N=8的分解过程。图中最右边的为8个时域采样点的8点FFTX[k]，由偶编号采样点的4点FFT和奇编号采样点的4点得到。这4点偶编号又由偶编号的偶采样点的2点FFT和奇编号的偶采样点的2点FFT产生。相同的4点奇编号也是如此。依次往左都可以用相同的方法算出，最后由偶编号的奇采样点和奇编号的偶采样点的2点FFT算出。图中没2点FFT成为蝶形，第一级需要每组一个蝶形的4组，第二级有每组两个蝶形的两组，最后一级需要一组4个蝶形。 实验内容 1.定义函数disbutterfly，程序根据FFT的定义：、，将序列x分解为偶采样点y和奇采样点z。 function [y,z]=disbutterfly(x) N=length(x); n=0:N/2-1; w=exp(-2*1i*pi/N).^n; x1=x(n+1); x2=x(n+1+N/2); y=x1+x2; z=(x1-x2).*w; 2.定义函数rader，纠正输出序列的输出顺序。 function y=rader(x,N) n=[0:N-1]; bn=dec2bin(n); rbn=fliplr(bn); rn=bin2dec(rbn); y=x(rn+1); 3.定义函数myfft，程序中套了两个循环。 function X=myfft(x) N=length(x); h=log2(N); %h=3 for i=1:h %第一次i=1;第二次i=2 s=[]; for j=1:2^(i-1);%i=1时，j=1;i=2时,j=1:2 M=2^(h-i+1);%M:M=8；M=4 xj=x([1:M]+(j-1)*M);%xj=x([1:8]+(1-1)*8)=x(1)+x(2)...+x(8); %j=1:xj=x([1:4]);j=2:xj=x([1:4]+4) [y,z]=disbutterfly(xj); s=[s,y,z]; end x=s; end X=rader(x,N); 4.主程序，将myfft与fft相减，比较之间的误差。 a=[1,2,3,4,5,6,7,8]; X=fft(a); X1=myfft(a); X0=fft(a)-myfft(a); subplot(4,1,1); stem(a); title(a序列); subplot(4,1,2); stem(X); title(a序列的fft); subplot(4,1,3); stem(X1); title(a序列的myfft); subplot(4,1,4); stem(X0); title(fft(a)-myfft(a)); 图中可看出fft与myfft的图几乎一模一样，且fft-myfft所得到的值几乎为零（虽然在4时有不等于零的结果，但是其值为数量级，误差极小）。 实验结论和讨论 通过自行编写FFT的运行程序，更加深入的了解了FFT的运算方法。如果设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”，那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多时，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数