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

一、实验目的

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

2.比较编写的myfft 程序运算结果与MATLAB 中的FFT 的有无误差。

二、实验条件

PC 机，MATLAB7.0

三、实验原理

1. FFT （快速傅里叶变换）原理：

将一个N 点的计算分解为两个N/2点的计算，每个N/2点的计算再进一步分解为N/4点的计算，以此类推。根据DFT 的定义式，将信号x[n]根据采样号n 分解为偶采样点和奇采样点。设偶采样序列为y[n]=x[2n]，奇采样序列为z[n]=x[2n+1]。

上式中的k N W -为旋转因子N k j e /2π-。下式则为y[n]与z[n]的表达式：

2. 蝶形变换的原理：

下图给出了蝶形变换的运算流图，可由两个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 的定义：]2[][][N n x n x n y +

+=、n N W N n x n x n z -+-=])2

[][(][，将序列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时有不等于零的结果，但是其值为15-10数量级，误差极小）。

五、实验结论和讨论

1.通过自行编写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=k2,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+（N^2）/2。继续上

面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，点数越多，运算量的节约就越大，这就是FFT的优越性。

2.本次实验中的myfft函数由两个子函数：disbutterfly和rader以及两个for循环组成。这两个for循环分别实现了对于输入矩阵x[n]需要进行h次的蝶形运算和用于调用每次蝶形运算所需要的x[n]序列的第某号数。

3.这次实验是在老师的一步一步指导下完成的，在跟着老师的思路的完成程序的时候，我意识到自己虽然能意识到应该写的程序的轮廓以及程序要实现的功能，却不知道该如何将这种思想变为计算机的逻辑思想。这让我知道自己对计算机思维转换的不熟练，需要更多的练习来巩固和完善自己。