clc;
clear;
%FFT
tic;
m=10;
N=2^m;
dt=0.01;
%將函數離散化輸入
for k=0:1:(N-1)
    xk(k+1)=k*dt;
    %xk(k+1)=2*pi*k/N;
    fxk(k+1)=exp(-abs(xk(k+1)));
end
%plot(fxk)；

%矩陣計算
T=[1,1;1,-1];
for j=1:1:m;
    W(j)=exp(-2*pi*i/2^j);
end
for j=2:1:m
    for n=1:1:2^(j-1)
        diagonal(n)=W(j)^(n-1);
    end
    K=diag(diagonal);
    I=eye(2^(j-1));
    T=[I,I;I,-I]*[I,zeros(size(I));zeros(size(I)),K]*...
    [T,zeros(size(T));zeros(size(T)),T];
end

%二進制下標換位
K=zeros(N,m);%倒序的二進制腳標
for n=0:1:(N-1)
    a=n;
    for j=1:1:m
        K(n+1,j)=mod(a,2);
        a=(a-K(n+1,j))/2;
    end
end
k_mark=zeros(1,N);%倒序后的二進制腳標轉換為十進制
for n=1:1:N
    for j=1:1:m
        k_mark(n)=k_mark(n)+K(n,j)*2^(m-j);
    end
end
for n=1:1:N%倒序
    Ak(n)=fxk(k_mark(n)+1);
end

Ck=T*Ak';%DFT;
t1=toc;%運行時間

tic;
ff=fft(fxk);
t2=toc;%matlab中fft運行時間