%該程序實現快速傅立葉變換
%輸入參數可以是行或列向量
%函數返回值為x的傅立葉變換值及其做變換所用的點數N
function [y1,N]=ditfft(x)

for p=1:inf  %找x的長度對應的2的最低冪次
    if 2^p>=length(x);break;end
end 

N=2^p;

if size(x,1)>1, x=x';end %若為列向量則轉換為行向量

if length(x<N)
  x1=[x,zeros(1,N-length(x))];%補零
end

handle1=@flip1;%調用flip1倒序函數
y1=feval(handle1,x1,N);

for p1=1:p     %將DFT作P次分解，依次計算
 
 nm=2^p1;u1=1;
 wn=exp(-i*2*pi/nm); %本次分解的基本因子wn
  
 for j=1:nm/2       %本次跨越間隔內的各碟形運算
       
     for k=j:nm:N  %本次碟形運算的跨越間隔為nm=2^p1
           kp=k+nm/2; %k,kp為碟形運算的對應單元下標
           t=y1(kp)*u1;
           y1(kp)=y1(k)-t;
           y1(k)=y1(k)+t;
       end
    
       u1=u1*wn;
       
   end
   
end