myfft.m

用蝶形算法和码位倒置法编写的fft变换程序

% 用matlab编的fft变换程序
% 
% 
% 用蝶形算法和码位倒置法编写的fft变换程序，并验证之，有兴趣的可以看看。
function data=myfft(datat,nn,isign)
datat=mybitrevorder(datat,nn);                %这个地方也可以用matlab自带的bitrevorder,有兴趣的还是自己编一下
for i=0:length(datat)-1
    data(2*i+1)=datat(i+1);
    data(2*i+2)=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=nn.*2;
mmax=2;
while n>mmax;
    istep=2.*mmax;
    theta=6.28318530717959/(isign*mmax);
    wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
    for m=1:2:mmax
        for i=m:istep:n
            j=i+mmax;
   tempr=wr*data(j)-wi*data(j+1);
   tempi=wr*data(j+1)+wi*data(j);
   data(j)=data(i)-tempr;
   data(j+1)=data(i+1)-tempi;
   data(i)=data(i)+tempr;
   data(i+1)=data(i+1)+tempi;
        end
        wtemp=wr;
        wr=wtemp*wpr-wi*wpi+wr;
  wi=wi*wpr+wtemp*wpi+wi;
    end
    mmax=istep;
end  
%%%%%%%%%%%%%码位倒置程序%%%%%%%%%%%%%%%%%%%%
function  A=mybitrevorder(A,N)
j=N./2;
for i=2:1:N-1
    if i<j
        temp=A(i);
        A(i)=A(j+1);
        A(j+1)=temp;
    end
    k=N/2;
    while j>=k
        j=j-k;
        k=k/2;
    end
    j=j+k;
end

%%%%%%%%%%%%%进行验证(matlab自带的fft，自己编写的fft，和函数的解析fft三者的比较%%%%%%%
clear all;
isign=1;
nn=256;
t=linspace(0,3,nn);
f=2.*exp(-3.*t);
data=myfft(f,nn,isign);
Ts=t(2)-t(1);
Ws=2.*pi./Ts;
W=Ws.*(0:nn./2)./nn;
for i=0:nn-1
    datar(i+1)=data(2.*i+1).*Ts;
    datai(i+1)=data(2.*i+2).*Ts;
end
for i=1:nn./2+1
    co(i)=datar(i);
    si(i)=datai(i);
    power(i)=(co(i).^2+si(i).^2).^0.5;
end
F=fft(f);
Fp=F(1:nn./2+1)*Ts;
Fpr=real(Fp);
Fpi=imag(Fp);
W=Ws.*(0:nn./2)./nn;
A=abs(Fp);
Fa=2./(3+j.*W);
B=abs(Fa);
plot(W,B,'o',W,A,'+r',W,power,'g');