%該程序?qū)崿F(xiàn)快速傅立葉變換
%輸入?yún)?shù)可以是行或列向量
%函數(shù)返回值為x的傅立葉變換值及其做變換所用的點(diǎn)數(shù)N
function [y1,N]=ditfft(x)

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

N=2^p;

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

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

handle1=@flip1;%調(diào)用flip1倒序函數(shù)
y1=feval(handle1,x1,N);

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