myfft2.m

采用DIT和DIF两种蝶形算法实现的FFT变换

function y=myfft2(x_ori,M)
%---------------------------------------输入原始信号，正弦函数频率为5--------------------------------%
% t=linspace(1,100,512)/100;                                                     %取了512个点
% x_ori=sin(2*pi*t*5)';                                                             %抽样得到原始函数的离散形式
% M=9;
N=2^M;
x=zeros(N,M+1);
x(:,1)=x_ori;
%-------------------------------------------此时开始FIT变换-------------------------------------------------%
for m=2:(M+1)                                                                     %m指要得到的蝶形单元m-2右侧结果在x中的列数为m
%列举常用变量
    zushu=2^(m-2)                                                                 %这是在第m列的左边的蝶形单元的组数
    jg=2^(M-m+2)                                                                  %这是组与组之间的间隔
    jg2=2^(M-m+1)                                                                %这四每组一半的间隔
    a=zeros(1,N);    
%首先求出每一级的W因子来，用a(p)顺序代表，每级a(p)都更新
p=1;
for aa=1:jg2                                                                        %总共有jg2个W，比如m=3，jg2=8，8个W               
    a(p)=W(aa-1,jg);                                                             %得到exp(-2*pi*j*(aa-1)/(2^(M-m+2)))
    p=p+1;                                                                           %标号加一，循环得到下一个W
end
%接着可以利用W来求fft了
cc=1;    
for zu=1:zushu                                                                    %第zu组
%++++++方法一++++++%
        pp=1;
         for zu_in=1:jg2
 %        k=1+(zu-1)*jg+zu_in-1;                                         %计算过程
             k=(zu-1)*jg+zu_in;                                               %这是每一组首地址+zu_in
             q=k+jg2;                                                              %这是一个蝶形单元两个之间的距离，是前一级的，所以为m-2
             x(k,m)=x(k,m-1) +x(q,m-1);                                 %每一组的前一半元素
          end
         for zu_in=(jg2+1):jg
%         k=1+(zu-1)*jg+zu_in-1;                                        %计算过程
             k=(zu-1)*jg+zu_in;                                              %这是绝对位置坐标
             p=k-jg2;
              x(k,m)=a(pp)*(x(p,m-1)-x(k,m-1)) ;                    %每一组的后一半元素,如果是M+1列,就不用乘了
            pp=pp+1;
         end
         p=1;
%++++++方法二++++++%
% p=1; 
%          for h=cc:(cc+(jg2-1))                              %实际是每组开头到结尾
%              x(h,m)=x(h,m-1)+x(h+jg2,m-1);          %第k列第h个元素是k-1列h个元素和因子乘以。。。
%              x(h+jg2,m)=a(p)*(x(h,m-1)-x(h+jg2,m-1)); 
%              p=p+1; 
%          end 
%          cc=jg*zu+1; 
end
end
%至此得到了次序不对的X(k)
%--------------------------------------将频域顺序求码位倒置后的频域顺序-----------------------------%
label_last=[0:N-1];                                                                 %输出端的标号顺序
label_binary=dec2bin(label_last);                                            %取二进制数
for j=1:N
    for jj=1:(size(label_binary,2)-1*rem(N,2))/2                         %如果N是奇数，就是0:(N-1)/2,如果偶数，就是0:N/2
        mid=label_binary(j,jj);
        label_binary(j,jj)=label_binary(j,size(label_binary,2)-jj+1);
        label_binary(j,size(label_binary,2)-jj+1)=mid;
    end                                                                                    %至此得到了时域信号的顺序label_ori
end
label_ori=bin2dec(label_binary);                                            %进行转换后原始时域信号放到x第一列
xk=rand(N,1);
for ll=1:N
xk(label_ori(ll)+1)=x(ll,M+1);
end                                                                                       %至此获得可利用蝶形单元的原始信号
%-------------------------------------------此时画出结果-------------------------------------------------------%
S=abs(xk); 
plot((0:49)/50*50,abs(S(1:50)),'linewidth',2)