[筆記]以Matlab 實現cooley-tukey快速傅立葉演算法

這學期在學matlab
大學完全沒用過...
因此以下若有錯誤可以提出糾正

日前
在網路上都找不太到用matlab寫的cooley-tukey快速傅立葉演算法(其實是作業寫不出來上網查)

搞了半天就寫出來了

這邊是設定輸入序列的長度需為2的次方(因為沒做零填補運算)

並以時域間取法來做計算
具體的流程與公式可參考wiki

程式碼如下:
(輸入序列A輸出轉換過的序列X)

clear all;

A=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15];
% Fa=fft(A); %檢查用
L=length(A);

%先對序列做排序
inx=swap_index(L);
for i=1:L
    X(i)=A(inx(i)+1);
end

K=0;
index=1;
for d=1:log2(L) %d代表流程數 
    for e=1:L/(2^d) %每流程的對稱點數
        for g=1:2^(d-1) %上群
            W(index)=X(index)+exp((-1i*2*pi*K)/L)*X(index+(2^(d-1)));
            index=index+1;
            K=K+(L/(2^d));
        end
        for g=1:2^(d-1) %下群
            W(index)=X(index-(2^(d-1)))+exp((-1i*2*pi*K)/L)*X(index);
            index=index+1;
            K=K+(L/(2^d));
        end
        K=0;
    end
    for t=1:L %更新
        X(t)=W(t);
    end
    index=1;
end

swap_index()

function [ inx ] = swap_index( L )
    dep=log2(L);
    for i=0:L-1
        temp=dec2bin(i,dep)-48;
        a=0;
        for j=dep-1:-1:0
            a=a+((2^j)*temp(j+1));
        end
        inx(i+1)=a;    
    end
end

簡單的說

就是先將輸入的序列做排序
排序的方式為
先將索引轉成2進制然後反轉
再轉回十進制
產生新的索引序列
接著再依照公式做利用規則與對稱性做運算