以前从没手动实现过FFT
参考
% 2012.3.10
function test
clear;
clc;
close all;
nFFT=2^10;
N=nFFT;
t=linspace(0,1,N)';
fs=N;
x=sin(2*pi*300*t);
freq=fs*((1:nFFT)-1)/nFFT; % bug fixed 1:nFFT
figure;
disp('time[Matlab fft]:');
tic
y0=fft(x,nFFT);
toc
plot(freq,abs(y0));
title('fft');
disp('time[recursive]:');
tic
y=my_fft(x,nFFT);
toc
figure;
plot(freq,abs(y));
title('my-fft');
figure(3);
plot(freq, abs(y-y0));
title('my-fft-error');
disp('time[inplace]:');
tic
y_inplace=my_ffy_inplace(x,nFFT);
toc
figure(4);
plot(freq,abs(y_inplace-y0));
title('inplace-error');
figure(5);
plot(freq,abs(y_inplace));
function x=my_ffy_inplace(x,n)
J=zeros(n,1);
for k=1:n
J(k)=bit_reverse(n,k-1);
end
J=J+1; % 1-base
x=x(J);
t=log2(n);
for q=1:t
L=2^q; r=n/L; L_2=L/2;
for j=0:L_2-1
w=cos(2*pi*j/L)-1i*sin(2*pi*j/L);
for k=0:r-1
tau=w*x((k)*L+j+1+L_2);
x((k)*L+j+1+L_2)=x((k)*L+j+1)-tau;
x((k)*L+j+1)= x((k)*L+j+1)+tau;
end
end
end
end
function j=bit_reverse(n,k)
% 0-base index
t=log2(n);
j=0;
for q=1:t
s=floor(k/2);
b_q=k-2*s;
% horner
j=2*j+b_q;
k=s;
end
end
function y=my_fft(x,n)
if(n==1)
y=x;
else
m=n/2;
w=exp(-2*pi*1i/n);
Omega=zeros(m,m);
for j=1:m
Omega(j,j)=w^(j-1);
end
Zr=my_fft(x(1:2:n),m);
Zb=Omega*my_fft(x(2:2:n),m);
Im=eye(m);
y=[Im, Im; Im,-Im]*[Zr;Zb];
end
end
end