matlab怎么测FFT,FFT算法MATLAB实现与测试 | 学步园

以前从没手动实现过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