myfft.m

用matlab语言

function data=myfft(datat,nn,isign)
datat=mybitrevorder(datat,nn);
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;                     %也可以用直接写的
          [A(j+1),A(i)]=deal(A(i),A(j+1));  %可以用deal function%%%
    end
    k=N/2;
    while j>=k
        j=j-k;
        k=k/2;
    end
    j=j+k;
end