%自己定義的FFT函數，算法采用的就是課本上講的雷德算法

function y = FFT_self(x,M)

%--------------------------------------------------------------------------
%下面這部分是當輸入序列的長度不為2的指數冪時在序列后補零達到2的指數冪的長度
q = nextpow2(M);                          %返回第一個q使得 2^q >= abs(M)
N = 2^q;
if M<N
    x=[x,zeros(1,N-M)];                   %在序列后面補零
end
%--------------------------------------------------------------------------
%實際的FFT算法
%--------------------------------------------------------------------------
%這一部分是用于變址
nv2=N/2;
nm1=N-1;
i=0;                                      %用于表示自然順序的輸入數據
p=0;                                      %用于表示倒位序的數據
while i<nm1                        
    if i<p                                %當i<p時要進行變址，在i>p時不需要變址，可以避免重復換址
        t=x(p+1);
        x(p+1)=x(i+1);
        x(i+1)=t;
    end
    k=nv2;                                %nv2二進制的高位為1，其他位為0
%---------------------------------已知某個倒位序數，求下一個倒位序數
    while k<=p                            %P的高位為1，減去N/2將高位變為0
        p=p-k;
        k=k/2;                            %用于次高位的處理
    end
    p=p+k;                                %當P的高位為0時，直接加N/2變為1
    i=i+1;
end
%--------------------------------------------------------------------------
%L級遞推計算
l=1;                              
m=log2(N);                                %蝶形運算的級數
while l<=m                        
    le=2^l;                               %各蝶形結依次的距離
    le1=le/2;                             %同一種蝶形結構中參加運算的兩節點的間距
    u=1;                                  %開始的exp(-j*2*pi*0/N)
    w=exp(-j*2*pi/le);                    %每一級的系數
    p=0;
    while p<le1;                          %這層循環控制不同旋轉因子的碟形結
        i=p;
        while i<N                         %這層循環控制同一種蝶形結的運算
            ip=i +le1;                    %i為蝶形結中一個點，ip是另一個點
            t=u*x(ip+1);                  %下面三句就是進行相應的蝶形運算
            x(ip+1)=x(i+1)-t;
            x(i+1)=t+x(i+1);
            i=i+le;                       
        end
        u=u*w;                            %改變不同蝶形結的旋轉因子
        p=p+1;
    end
    l=l+1;
end
y = x;                                    %把經過FFT處理后的輸入序列x傳遞給輸出序列y