fft.txt

基2－FFT算法的Matlab代码的实现。其中有一个子程序可以用Matlab内部的程序来实现

function X=MYFFT(x,N)
%x:series;N,the number of fft-complex required,which must be a power of2
%X:output the fft of x
%-----------------------------------------------------------
 nargchk(1,2,nargin);
 if nargin<2%
     N=length(x);
 elseif N<=0%input must be a nonnegtive integer
     error('N must be a nonnegtive integer');
 end
num=log2(N);%N is a power of 2 ???
if (num/ceil(num)~=1)%case:No
    disp('N isnot a power of 2!');%message
    N=(2^ceil(num));%modify N
%     x=[x,zeros(1,N-length(x))];%zero-padding
end
if length(x)<N
    disp('Zero-padding was used!');
    x=[x,zeros(1,N-length(x))];%zero-padding
end
[m,n]=size(x);
if m<n
    x=x.';
end
%-----------------------------------------------------------
id=bitinv(N);%call subprogram
x_temp=x(id+1);%reshaped series
num=log2(N);%the No. of recycle
NUM=N;%the total No. of x that want to be compute
for K=1:num
    n=2^K;%xlh's or xfh's length
    idpublic=1:2^K:N;%public id
    Xbin=[];%initial
    for I=1:NUM/2%divide x to NUM/2 parts
        idbegin=idpublic(I);%2*(I-1)+2^(K-1);%first id
        tn=n/2;%id1=;id2=id1+n/2;
        id1=idbegin+(0:(tn-1));id2=idbegin+((tn):2*tn-1);
        xfh=x_temp(id1);
        xlh=x_temp(id2);
        Xbin(1:id2(end),1)=[Xbin;fft2_subs(xfh,xlh)];
     end
NUM=floor(NUM/2);
x_temp=Xbin;%x_temp changed
end
%-----------------------------------------------------------------    function Xbins=fft2_subs(xfh,xlh)
        r=length(xfh);ns=2*r;% in order to decide Wn
        Wn=exp(-j*2*pi/ns);
        k=(0:length(xlh)-1)';
        Xbins(:,1)=[xfh+xlh.*Wn.^k;xfh-xlh.*Wn.^k];
    end
%----------------------------------------------------------------    function id=bitinv(N)
%该程序也可以由Matlab中的bitrevord()代替
        %N:number,it must be a power of 2;
        %id,output the index of the sery to caculate fft
        if N<0
            error('input must be a nonnegtive integer!');
        end
        n=ceil(log2(N));%No. of bits needed to represent the index
        index=0:2^n-1;%index of the series
        binvec=zeros(1,n);
        id=binvec;
        for i=1:2^n
            binvec(i,:)=fliplr(dec2binvec(index(i),n));
            id(i)=binvec2dec(binvec(i,:));
        end
    end
%------------------------------------------------------------
X=Xbin;%
end