1.编写的fft函数
function X=FFT(x,N)
b=length(x); %x(n)的长度
%若x(n)的长度小于N,则补零
if(b<N)
x=[x,zeros(1,N-b)];
end
b=length(x); %补零后x(n)的长度
xm=zeros(1,b); %x(n)转换数组
M=log2(N);
%将x(n)的序号码位倒置
nxd=fliplr(dec2bin([1:N]-1,M));%将系数取二进制并对调
nxd=bin2dec(nxd)+1; %取十进制 (matlab的下标从1开始,需要加1)
x=x(nxd); %将x(n)倒序排列作为输入端的初始值
% N个1点DFT,即x(n),赋值到数组xm
for k=1:N
xm(k)=x(k);
end
%将N点DFT作M次基2分解,从左到右,对每次分解作DFT运算,共M级蝶形运算
for L=1:M
B=2^(L-1); %第L级中,每个蝶形单元两个输入数据的距离,共有B个旋转因子
for J=0:B-1 %第L级中不同的旋转因子
R=J*2^(M-L); %旋转因子的指数
for p=(J+1):2^L:N %本次蝶形运算的跨越间隔为2^L
W=exp(-j*2*pi*R/N); %计算本次运算的旋转因子
t=xm(p)+xm(p+B)*W; %进行蝶形运算
xm(p+B)=xm(p)-xm(p+B)*W;
xm(p)=t;
end
end
end
X=xm; %输出X为xm
2.与MATLAB内部的fft函数对比
clear
x=[1 1 1 1 1 1 1 1];
N=128;
X1=fft(x,N); %MATLAB内部的fft函数
X=FFT(x,N); %调用自己编写的fft函数
subplot(211);
stem(abs(X1),'.');grid on;
title('MATLAB内部的FFT变换图像');
subplot(212);
stem(abs(X),'.');grid on;
title('自己编写的FFT变换图像');
3.运行结果