FFT_matlab算法实现与验证
一、算法代码：
DIT_FFT_algorithm:
clear,clc,
clear all;
xn=[0,1,2,3,4,5,6,7];
N=length(xn);
A=xn;
%DIT_FFT
NI=N/2;
for I=1:N-1
if I<NI
t=A(I+1);
A(I+1)=A(NI+1);
A(NI+1)=t;
end
T=N/2;
while NI>=T
NI= NI-T;
T=T/2;
end NI= NI+T;
end
disp('逆序x[n]:'),disp(A);
%butterfly
WN=exp(-i*2*pi/N);
v=floor(log2(N));
for m=1:v
for k=0:2^m:N-1
for K=0:2^(m-1)-1
p=k+K;
q=p+2^(m-1);
r=2^(v-m)*mod(p,2^m);
B(p+1)=A(p+1)+A(q+1)*WN^r;
B(q+1)=A(p+1)-A(q+1)*WN^r;
end
end A=B;
end
disp('FFT_X[k]:'),disp(A);
DIF_FFT_Algorithm:
clear,clc,
clear all;
xn=[0,1,2,3,4,5,6,7];
N=length(xn);
A=xn;
%DIF_FFT
v=floor(log2(N));
WN=exp(-i*2*pi/N);
for m=1:v
for k=0:2^(v-m+1):N-1
for K=0:2^(v-m)-1
p=k+K;
q=p+2^(v-m);
r=2^(m-1)*mod(p,2^(v-m+1));%基于DIT_FFT的修改；
B(p+1)=A(p+1)+A(q+1);
B(q+1)=(A(p+1)-A(q+1))*WN^r;
end
end
A=B;
disp(A);
end
NI=N/2;
for I=1:N-1
if I<NI
t=A(I+1);
A(I+1)=A(NI+1);
A(NI+1)=t;
end
T=N/2;
while NI>=T
NI=NI-T;
T=T/2;
end NI=NI+T;
end
disp('X[k]:'),disp(A);
我们选择xn=[0,1,2,3,4,5,6,7]，通过matlab自带的fft函数验证正确性，具有一般的代表性（因为一开始通过x[n]=[1,0,0,0,0,0,0,0]，虽然结果得到X[k]=[1,1,1,1,1,1,1,1]，但当改变输入离散信号时得到的结果却和fft函数相差甚远）：
%test_fft
xn=[0,1,2,3,4,5,6,7];
A=fft(xn);
disp(A);
最终得到的结果是相同的：
DIT_FFT:
最终Xk：
Columns 1 through 5
28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000
Columns 6 through 8
-4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
DIF_FFT:
最终Xk：
Columns 1 through 5
28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000
Columns 6 through 8
-4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
matlab_fft：
Columns 1 through 5
28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000
Columns 6 through 8
-4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
其次，根据DFT公式，我们同样可以验证结果：
%DFT
clear all;
xn=[0,1,2,3,4,5,6,7];
N=length(xn);
k=0:N-1;
n=0:N-1;
WN=exp(-i*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
disp(Xk);
同样可以得到相同的结果。
二、算法结构分析
以DIT_FFT 为例，
1. 逆序，根据bit_reverse 规则，将x[n]每一位数值按照二进制地址（log2(N)位）进行对称
逆运算，得到的逆序序列在进行butterfly 变换。
为防止在交换过程中出现重复交换，使得不到结果，规定只有地址小的与地址大的数值交换。
2. “碟式”运算，根据如下规则：
1[][][]r m m N m X p X p W X q +=+
1[][][]r m m N m X q X p W X q +=-
12m q p -=+
2[mod(2)]v m m r p -=*
其中， 2,v
N =m 表示运算级次，m=1,2…v.
总体分为三重循环：
1) 共进行v 级butterfly 运算，令m=1:v ；
2) 第m 级运算中，又分为zu=2^(v-m)组butterfly 运算，每组的运算数据长度为L=2^m ，
每一组的开始数据地址为dizhi=0+k*(2^m), k=0,1,…zu ；
3) 在每一组中，又分为2^(m-1)次运算，每一次运算采用如上图所示的计算法则，得
到该级该组该次的butterfly 运算结果。
3. 最终将FFT 变换的结果顺序输出，得到X[k]。
DIF_FFT 则是先进性butterfly 运算，然后在对错序的X ’[k]进行bit_reverse 操作。
三、图表说明
当输入离散数值x[n]=[1,1,1,1,1,1,1,1]时，易知[]0,0X k k =≠
运行程序结果如下：
x[n]:
1 1 1 1 1 1 1 1
X[k]:
8 0 0 0 0 0 0 0
四、解决任意N 问题
对于任意的N 数据，如果不是2的整数次幂，我们可以通过补零，使数据满足butterfly 运算的要求，程序如下：
%%%%补零
clear,clc,
clear all;
xn=[1,0,0,0,0,0,0,0,0,0,0]; %%%%%%
N1=length(xn);
n=floor(log2(N1));
if 2^n<N1
n=n+1;
end
N=2^n;
A=[xn,zeros(1,N-N1)];
%余下程序段相同。
