fft.m

数值分析,通过自身递归调用实现奇偶分离,用于一维快速傅立叶程序变换代码的实现

tic
p=10;
N=2^p;
a1=zeros(1,N);a2=zeros(1,N);
w=zeros(1,N/2);f=zeros(1,N);
xk=zeros(1,N);ck=zeros(1,N);
h=2*pi/N;
xi=[0.0:h:(N-1)*h];
for k=1:N
   f(k)=xi(k)*exp(-xi(k));
   xk(k)=f(k)/N;
end
for k=1:N
   a1(k)=xk(k);
end
for m=0:N/2-1
   w(m+1)=cos(-m*h)+sin(-m*h)*i;
end
for q=1:p
   if(mod(q,2)~=0)
      for k=0:2^(p-q)-1
         for j=0:2^(q-1)-1
            a2(k*2^q+j+1)=a1(k*2^(q-1)+j+1)+a1(k*2^(q-1)+j+2^(p-1)+1);
            a2(k*2^q+j+2^(q-1)+1)=(a1(k*2^(q-1)+j+1)-a1(k*2^(q-1)+j+2^(p-1)+1))*w(k*2^(q-1)+1);
         end
      end
  else
      for k=0:2^(p-q)-1
         for j=0:2^(q-1)-1
            a1(k*2^q+j+1)=a2(k*2^(q-1)+j+1)+a2(k*2^(q-1)+j+2^(p-1)+1);
            a1(k*2^q+j+2^(q-1)+1)=(a2(k*2^(q-1)+j+1)-a2(k*2^(q-1)+j+2^(p-1)+1))*w(k*2^(q-1)+1);
         end 
      end
   end
end
if(mod(p,2)==0)
   ck=a1;
else
   ck=a2;
end
toc
 xq=[-2*pi:h:2*pi];
 NN=length(xq);
 s=zeros(1,NN); re=zeros(1,NN);
for m=1:NN
    for k=1:N
        s(m)=s(m)+ck(k)*(cos((k-1)*xq(m))+i*sin((k-1)*xq(m)));
    end
 end
re=real(s);
plot(xq,re,'r')