按时间抽取的基2FFT算法分析及MATLAB实现

按时间抽取的基2FFT 算法分析及MATLAB 实现
一、DIT-FFT 算法的基本原理
基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT ，再进行适当的组合，得到原N 点序列的DFT ，最终达到减少运算次数，提高运算速度的目的。

按时间抽取的基2FFT 算法，先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT 运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N 点序列的DFT 。

只要N 是2的整数次幂，这种分解就可一直进行下去，直到其DFT 就是本身的1点时域序列。

图1 DIT-FFT 蝶形运算流图
二、DIT-FFT 算法的运算规律及编程思想
1.原位计算
[
对N=M
2点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。

在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M 级运算后，原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值，这种原位(址)计算的方法可节省大量内存。

2.旋转因子的变化规律
N 点DIT ―FFT 运算流图中，每个蝶形都要乘以旋转因子p
W N ，p 称为旋转因子的指数。

例如N ＝8 ＝3
2 时各级的旋转因子：
第一级：L=1， 有1个旋转因子：p
W N =J
/4W N =J
2L W J=0
第二级：L=2，有2个旋转因子：p
W N =J
/2W N =J
2L W J=0,1
第三级：L=3，有4个旋转因子：p
W N =J
W N =J 2L W J=0,1,2,3
对于N ＝M 2的一般情况，第L 级共有1-L 2个不同的旋转因子：
p W N =J 2L W J=0,1,2,… ,1-L 2
－1 L 2=M 2×M -L 2= N ·M -L 2
故： 按照上面两式可以确定第L 级运算的旋转因子
\
3、同一级中，同一旋转因子对应蝶形数目 第L 级FFT 运算中，同一旋转因子用在L
-M 2
个蝶形中；
4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离” 第L 级中，蝶距：D=L
2； 5、同一蝶形运算两输入数据的距离
在输入倒序，输出原序的FFT 变换中，第L 级的每一个蝶形的2个输入数据相距：B=1
-L 2。

6、码位颠倒
输入序列x(n)经过M 级时域奇、偶抽选后，输出序列X(k)的顺序和输入序列的顺序关系为倒位关系。

'
将十进制顺序数用I 表示，与之对应的二进制是用IB 表示，十进制倒序数用J 表示，与之对应的二进制是用JB 表示。

十进制顺序数I 增加1，相当于IB 最低位加1且逢2向高位进1，即相当于JB 最高位加1且逢2向低位进1。

JB 的变化规律反映到J 的变化分为两种情况，若JB 的最高位是0（J<N/2），则直接由加1（J ←J+N/2）得到下一个倒序值，若JB 的最高位是1（J ≧N/2），则要先将最高位变0（J ←J-N/2），再在次高位加1（J ←J+N/4），但次高位加1时，同样要判断0、1值，如果是0（J<N/4）,则直接加1（J ←J+N/4）,否则要先将次高位变0（J ←J-N/4）再判断下一位，依次类推，直到完成最高位加1，逢2向右进位的运算。

I=J 时不需要交换，只对I<J 时的情况进行数据交换即可，数据倒序程序框图如如2。

7、蝶形运算的规律
序列经过时域抽选后，存入数组中，如果蝶形运算的两个输入数据相距B 个点，应用原位计算，蝶形运算可表示成如下形式：
8、 DIT-FFT 程序框图
~
根据DIT-FFT 原理和过程，DIT-FFT 的完整程序框图如图2：
(1)倒序:输入自然顺序序列x(n)，根据倒序规律，进行倒序处理；
(2)循环层1:确定运算的级数，L=1~M (N=M
2)；确定一蝶形两输入数据距离B=1
-L 2
(3)循环层2:确定L 级的B=1
-L 2
个旋转因子；旋转因子指数p=J×L
-M 2
，J=0~B-1；
XL-1(J)
X L-1 (J+B)
XL (J)= XL-1(J)+ WNp X L-1
(J+B)
X L (J) = X L-1(J)W N p X L-1
(J+B)
p=J×2M-L ， J=0,1,2,… ,2L-1－1
(4)循环层3:对于同一旋转因子，用于同一级L -M 2个蝶形运算中：k 的取值从J 到N-1，步长为L 2 (使用同一旋转因子的蝶形相距的距离) (5)完成一个蝶形运算。

开 始
送入x (n )，M
N ＝2
M 倒 序
L ＝1 , M Ｊ＝0 , B － 1
P ＝2
M －L J k ＝ J , N －1 , 2
L p
N p N W B k X k X B k X W B k X k X k X )()()()()()(+-⇐+++⇐输 出
结 束
B 2 L －1
图2 数据倒序程序框图 图3 DIT-FFT 的完整程序框图
三、程序源代码
/
设计函数myDitFFT(xn)完成一个序列的DIT-FFT 运算：
function y=myDitFFT(xn) M=nextpow2(length(xn)); N=2^M;
disp('调用fft函数运算的结果:'),
fftxn=fft(xn,N);
if length(xn)<N
xn=[xn,zeros(1,N-length(xn))];
end
for m=0:N/2-1;%旋转因子指数范围
}
WN(m+1)=exp(-j*2*pi/N)^m;%计算旋转因子end
disp('输入到各存储单元的数据:'),disp(xn);
%数据倒序操作
J=0;%给倒序数赋初值
for I=0:N-1;%按序交换数据和算倒序数
if I<J;%条件判断及数据交换
T=xn(I+1);xn(I+1)=xn(J+1);xn(J+1)=T;
end
%
%算下一个倒序数
K=N/2;
while J>=K;
J=J-K;K=K/2;
end
J=J+K;
end
disp('倒序后各存储单元的数据:'),
disp(xn);
% 分级按序依次进行蝶形运算
·
for L=1:M;%分级计算
disp('运算级次:'),disp(L);
B=2^(L-1);
for R=0:B-1;%各级按序蝶算
P=2^(M-L)*R;
for K=R:2^L:N-2;%每序依次计算
T=xn(K+1)+xn(K+B+1)*WN(P+1);
xn(K+B+1)=xn(K+1)-xn(K+B+1)*WN(P+1);
xn(K+1)=T;
end
！
end
disp('本级运算后各存储单元的数据:'),disp(xn);
end
在主函数中调用myDitFFT(xn)函数实现DIT-FFT并和直接DFT运算结果做对比：
xn=[0,1,2,3,4,5,6,7];
myDitFFT(xn);
调用fft函数运算的结果:
1 至7 列
`
+ + + + + - -
8 列
-
调用myDitFFT(xn)函数运行的结果：
输入到各存储单元的数据:
0 1 2 3 4 5 6 7 倒序后各存储单元的数据:
0 4 2 6 1 5 3 7 (
运算级次:
1
本级运算后各存储单元的数据:
4 -4 8 -4 6 -4 10 -4 运算级次:
2
本级运算后各存储单元的数据:
1 至7 列
+ + + - + + +
8 列
[
-
运算级次:
3
本级运算后各存储单元的数据:
1 至7 列
+ + + + + - -
8 列
-
经对比可知DIT-FFT与直接DFT的运行结果完全相同。

四、总结
经过验证可发现DIT-FFT较直接DFT运算有着明显的优势，我们可以将这个函数运用在多个领域以简化运算，例如计算离散时间序列的卷积或计算IDFT时都可以应用到DIT-FFT算法，我感受到数字信号处理中科学思想的魅力。

由于对设计思路的缺乏，我在设计程序时，在网络上查找了很多有关DIT-FFT的资料，经过学习他人的解决思路最后才整理出DIT-FFT的程序，在有些地方我自己理解的还不是很透彻，比如在实现数据倒序的程序我认为比较困难；当然即使自己想不到能学习一下别人的思路也是很好的，这个程序的代码量并不大，我自身的能力还很低，要在以后的学习中不断进步才能完成更加复杂的任务。

这次课程设计让我对快速傅里叶变换有了更多的了解，也认识到了科学计算方法的重要性，我感到很充实。

参考文献——
百度百科；
按时间抽取的基2FFT算法分析及MATLAB实现[J].电子技术，2011（2）
数字信号处理（第四版）西安电子科技大学出版社高希全丁玉美编)。