radix-2fft.txt

基2FFT,是DFT(离散傅立叶变换)的快速算法

N=input('输入2的幂N='),n=0:N-1,r=log2(N);                   %输入N  
xn=...                                                                                %输入x[n]

k=0:N-1;
subplot(221),stem(n,real(xn)),title('x[n]实部'),xlabel('n');
subplot(222),stem(n,imag(xn)),title('x[n]虚部'),xlabel('n');         %画图

X=xn;                                                                             %先将X初始化，便于接下来的运算
for m=1:r;                                                                        %m表示第m级的蝶形流图计算
xn=X;                                                                             %在第m级蝶形流图计算前先将xn变为上一级蝶形流图计算的结果，便于接下来的一级蝶形流图计算
for k=0:2^(r-m)-1;
X(k+1)=xn(k+1)+xn(k+1+2^(r-m)),X(k+1+2^(r-m))=(xn(k+1)-xn(k+1+2^(r-m)))*exp(-j*2*pi*(k*2^(m-1))/N);
end;                                                                                  %此循环计算第一个交叉在一起的蝶形
while k+2^(r-m)~=N-1;g=k+2^(r-m)+1;                          %判断是否计算出了该级蝶形流图计算的所有结果，若没有则接着计算下一个交叉在一起的蝶形
for k=g:g+2^(r-m)-1;
X(k+1)=xn(k+1)+xn(k+1+2^(r-m)),X(k+1+2^(r-m))=(xn(k+1)-xn(k+1+2^(r-m)))*exp(-j*2*pi*(k*2^(m-1))/N);
end;
end;
end;

Xk=zeros(1,N);                                          %由于X为倒序位排列，所以要对其进行重新排序，并附给Xk，Xk便是最后的结果
for g=1:N;
a=zeros(1,r);
b=1;
c=g-1;
while b<=r;
a(b)=rem(c,2),b=b+1,c=fix(c/2);                %rem()为取余数，fix()为向0方向取整数
end;
a=fliplr(a),e=0;                                           %fliplr()为矩阵作左右翻转
for f=1:r;
e=e+a(f)*(2^(f-1));
end;
Xk(e+1)=X(g);                                            %MATLAB中序列的序号只能从1开始，所以这里要加1
end;


k=0:N-1;                                                                                          %因为此时的k不是数组，为了画图需要将k变为数组
subplot(223),stem(k,real(Xk)),title('Xk实部'),xlabel('k');
subplot(224),stem(k,imag(Xk)),title('Xk虚部'),xlabel('k');                                               %画图