FFT变换的Matlab程序

1.Main Program

clear all;

A1=1;f1=100;fi1=0.0;fi1=fi1*pi/180;

A2=1;f2=150;fi2=0.0;fi2=fi2*pi/180;

N=128;fs=400;

n=[0:N-1];

dt=n/fs;

x=A1*cos(2*pi*f1.*dt+fi1);

y=A2*cos(2*pi*f2.*dt+fi2);

x=x+i*y;

y=fft(x);

z=fft_lyn(x);

figure(1)

subplot(2,1,1), plot(dt,real(x));

subplot(2,1,2), plot(dt,imag(x));

figure(2)

subplot(2,1,1), plot(n,real(y));

subplot(2,1,2), plot(n,real(z));

figure(3)

subplot(2,1,1), plot(n,imag(y));

subplot(2,1,2), plot(n,imag(z));

%clear all;

2.Sub Program

function X=fft_lyn(x)

[tmp,N]=size(x);

X=x;

m=log2(N);

p=1;

%Prepare w;

w0=exp(-i*2*pi/N);

l=N/2-1;

for n=1:l

w(n)=w0^(n);

end

p=1;

Nk=N;

for k=1:m

hNk=Nk/2;hNk1=hNk+1;hNk2=hNk1+1;

%Butterfly;

for Nk0=0:Nk:N-1

for n=1:hNk

j=Nk0+n;

j1=hNk+j;

t=X(j);

X(j)=X(j)+X(j1);

X(j1)=-X(j1)+t;

end

end

%Twiddle;

for Nk0=0:Nk:N-1

l=1;

for n=hNk2:Nk

j=Nk0+n;

X(j)=X(j)*w(p*l);

l=l+1;

end

end

Nk=hNk;

p=p+p;

end

%Bit reverse order;

hN=N/2;hN1=hN+1;hN_1=hN-1; for n=1:N

Nk=hN;

l=n;

j=1;

Nj=1;

for k=1:m

if(l>Nk)

l=l-Nk;

j=j+Nj;

end

Nk=Nk/2;

Nj=Nj+Nj;

end

if(n

t=X(n);

X(n)=X(j);

X(j)=t;

end

end