%利用快速算法計算離散時間序列的傅立葉變換
function y=bianhuan(xr,n)
%參數說明 xr為單精度實型一維數組，長度為n;
%n為整型變量，是數據長度，也是進行FFT變換的點數，必須為2的整數次冪
for i=1:n
    xi(n)=0;%設定虛部為零
end
 j=0;      
 for i=0:n-2;
     %實現數組元素的交換
     if i<j
         t=xr(i+1);
         xr(i+1)=xr(j+1);
         xr(j+1)=t;
         t=xi(i+1);
         xi(i+1)=xi(j+1);
         xi(j+1)=t;
    end
     k=n/2;
     while (k<=j)
         j=j-k;
         k=k/2;
     end
     j=j+k;
 end
 %求出蝶形運算的級數
 ep=0;
 i=n;
 while i>1
     ep=ep+1;
     i=i/2;
 end
 ph=2*pi/n;
 
 for i=0:n-1;
     s(i+1)=sin(ph*i);
     c(i+1)=cos(ph*i);
 end
 a=2;
 b=1;
 for mt=1:ep;%共進行ep級運算，每級N/2個蝶形運算，N=2的ep次冪
     s0=n/a;
     s1=1;
     for k=0:b-1;
         i=k;
         while i<n;
             arg=i+b;
             if k==0;
                 pr=xr(arg+1);
                 ppi=xi(arg+1);
             else%確定pr、ppi，為迭代運算做準備
                 pr=xr(arg+1)*c(s1+1)-xi(arg+1)*s(s1+1);
                 ppi=xr(arg+1)*s(s1+1)+xi(arg+1)*c(s1+1);
             end%每個蝶形結構完成下列基本迭代運算
             
             xr(arg+1)=xr(i+1)-pr;
             xi(arg+1)=xi(i+1)-ppi;
             xr(i+1)=xr(i+1)+pr;
             xi(i+1)=xi(i+1)+ppi;
             i=i+a;
         end
         s1=s1+s0;
     end
     a=2*a;
     b=b*2;
 end

for k=1:n;
y(k)=sqrt(xr(k)*xr(k)+xi(k)*xi(k));  %求模
end