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

1、 实验目的1. 利用MATLAB编写FFT快速傅里叶变换。2. 比较编写的myfft程序运算结果与MATLAB中的FFT的有无误差。二、实验条件 PC机，MATLAB7.03、 实验原理1. FFT（快速傅里叶变换）原理： 将一个N点的计算分解为两个N/2点的计算，每个N/2点的计算再进一步分解为N/4点的计算，以此类推。根据DFT的定义式，将信号xn根据采样号n分解为偶采样点和奇采样点。设偶采样序列为yn=x2n，奇采样序列为zn=x2n+1。上式中的为旋转因子。下式则为yn与zn的表达式： 2. 蝶形变换的原理： 下图给出了蝶形变换的运算流图，可由两个N/2点的FFT（Yk和Zk得出N点FFT Xk）。同理，每个N/2点的FFT可以由两个N/4点的FFT求得。按这种方法，该过程可延迟后推到2点的FFT。下图为N=8的分解过程。图中最右边的为8个时域采样点的8点FFTXk，由偶编号采样点的4点FFT和奇编号采样点的4点得到。这4点偶编号又由偶编号的偶采样点的2点FFT和奇编号的偶采样点的2点FFT产生。相同的4点奇编号也是如此。依次往左都可以用相同的方法算出，最后由偶编号的奇采样点和奇编号的偶采样点的2点FFT算出。图中没2点FFT成为蝶形，第一级需要每组一个蝶形的4组，第二级有每组两个蝶形的两组，最后一级需要一组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=3for 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时有不等于零的结果，但是其值为数量级，误差极小）。5、 实验结论和讨论1. 通过自行编写FFT的运行程序，更加深入的了解了FFT的运算方法。如果设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”，那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多时，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+（N2）/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，点数越多，运算量的节约就越大，这就是FFT的优越性。2. 本次实验中的myfft函数由两个子函数：disbutterfly和rader以及两个for循环组成。这两个for循环分别实现了对于输入矩阵xn需要进行h次的蝶形运算和用于调用每次蝶形运算所需要的xn序列的第某号数。3.这次实验是在老师的一步一步指导下完成