bianhuan.m

实现FFT变换的函数。调用形式为y=bianhuan(xr,n) 其中

%利用快速算法计算离散时间序列的傅立叶变换
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