利用matlab编写fft快速傅里叶变换.doc

实验目的
利用MATLAB编写FFT快速傅里叶变换。
比较编写的myfft程序运算结果与MATLAB中的FFT的有无误差。
二、实验条件
PC机，
实验原理
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个蝶形。
实验内容
，程序根据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;
，纠正输出序列的输出顺序。
function y=rader(x,N)
n=[0:N-1];
bn=dec2bin(n);
rbn=fliplr(bn);
rn=bin2dec(rbn);
y=x(rn+1);
，程序中套了两个循环。
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);
，将myfft与fft相减，比较之间的误差