【Matlab】时域基2FFT

时域基2FFT

程序较为简陋，需要进一步编写加强灵活性。
要求：输入的序列和N的点数必须为2的幂。
Matlab：2019b

function [y_n] = MyFFT(x_n,N)
%MYFFT 时域基2的快速傅里叶变换
%   x_n:需要进行fft的序列
%   N  :做N点DFT，应当为2的指数幂
%   y_n:输出序列

M = log2(N);    %   M级蝶形运算
L = length(x_n);%   输入序列长度
x_input = 1:L;  %   倒序后的实际输入序列,建议直接赋值，用zeros之类的函数会定义为矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   倒序
last_index = 0;
x_input(1) = x_n(1);
for k=1:L-1       %   进行L-1次运算
    index = last_index;%   序列下标
    div = 2;    %   分母
    while 1
        if index<(N/div)
            index = index + N/div;
            x_input(k+1) = x_n(index+1);
            last_index = index;
            break
        else
            index = index - N/div;
            div = div * 2;
            continue
        end 
    end
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   蝶形运算
for k=1:M        %  进行M次蝶形运算,第k级运算
   J = 0:2^(k-1)-1;
   ro = exp(-j .* J .* 2^(M-k)*2*pi/N);   %   每级的旋转因子
   ro_index = 0;
   x_input_index = 1;   
   
   for conter = 1:N/2 %   每级应当进行N/2次蝶形计算
       stable = x_input(x_input_index);
       change = x_input(x_input_index+2^(k-1));
       temp1 = stable + change .* ro(rem(ro_index,2^(k-1))+1);
       temp2 = stable - change .* ro(rem(ro_index,2^(k-1))+1);
       x_input(x_input_index) = temp1;
       x_input(x_input_index+2^(k-1)) = temp2;
       
       if rem(conter,2^(k-1))==0
          x_input_index = x_input_index +2^(k-1)+1; % 蝶形运算满2^(k-1)次后加2^(k-1)+1个，否则加1
       else
          x_input_index = x_input_index + 1; 
       end
       ro_index = ro_index + 1;
   end
end

y_n = x_input;
end