myfft1.m

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

function y=myfft1(x_ori,M)
%---------------------------------------输入原始信号，正弦函数频率为5--------------------------------%
% t=linspace(1,100,512)/100;                                                     %取了512个点
% x_ori=sin(2*pi*t*5)';                                                              %抽样得到原始函数的离散形式
%--------------------------------------将频域顺序求码位倒置后的时域顺序-----------------------------%
% M=9;
N=2^M;
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第一列
x=rand(N,M+1);
a=zeros(1,N);
for ll=1:N
x(ll,1)=x_ori(label_ori(ll)+1);
end                                                                                       %至此获得可利用蝶形单元的原始信号
%-------------------------------------------此时开始FIT变换-------------------------------------------------%
for m=2:(M+1)                                                                     %m指要得到的蝶形单元m-2右侧结果在x中的列数为m
%列举常用变量
    zushu=2^(M-m+1);                                                          %这是在第m列的左边的蝶形单元的组数
    jg=2^(m-1);                                                                     %这是组与组之间的间隔
    jg2=2^(m-2);                                                                   %这四每组一半的间隔
    
%首先求出每一级的W因子来，用a(p)顺序代表，每级a(p)都更新
p=1;
for aa=1:jg2                                                                       %总共有jg2个W，比如m=3，jg2=2，两个W               
    a(p)=W(aa-1,jg);                                                           %得到exp(-2*pi*j*(aa-1)/(2^(m-1)))
    p=p+1;                                                                         %标号加一，循环得到下一个W
end
a'
cc=1;
%接着可以利用W来求fft了
for zu=1:zushu                                                   %第zu组
%++++++方法一++++++%
%         p=1;
%          for zu_in=1:jg2
% %          k=1+(zu-1)*jg+zu_in-1;                             %计算过程
%              k=(zu-1)*jg+zu_in;                                    %这是每一组首地址+zu_in
%              q=k+2^(m-2);                                           %这是一个蝶形单元两个之间的距离，是前一级的，所以为m-2
%              x(k,m)=x(k,m-1) +a(p)*x(q,m-1);              %每一组的前一半元素
%              p=p+1;
%          end
%          p=1;
%          for zu_in=(jg2+1):jg
%  %         k=1+(zu-1)*jg+zu_in-1;                            %计算过程
%              k=(zu-1)*jg+zu_in;                                   %这是绝对位置坐标
%              p=k-2^(m-2);
%              x(k,m)=x(p,m-1)-a(p)*x(k,m-1) ;         %每一组的后一半元素
%          end
%          p=1;

%++++++方法二++++++%

p=1; 
         for h=cc:(cc+(jg2-1))                              %实际是每组开头到结尾
             x(h,m)=x(h,m-1)+a(p)*x(h+jg2,m-1);   %第k列第h个元素是k-1列h个元素和因子乘以。。。
             x(h+jg2,m)=x(h,m-1)-a(p)*x(h+jg2,m-1); 
             p=p+1; 
         end 
         cc=jg*zu+1; 
end
end
%-------------------------------------------此时画出结果-------------------------------------------------------%
S=abs(x(:,M+1)); 
plot((0:49)/50*50,abs(S(1:50)),'linewidth',2)