按时间抽取的基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 运算流图中,每个蝶形都要乘以旋转因子pW N ,p 称为旋转因子的指数。
例如N =8 =32 时各级的旋转因子:第一级:L=1, 有1个旋转因子:pW N =J/4W N =J2L W J=0第二级:L=2,有2个旋转因子:p W N =J /2W N =J2L W J=0,1第三级:L=3,有4个旋转因子:p W N =J W N =J2L W J=0,1,2,3 对于N =M 2的一般情况,第L 级共有1-L 2个不同的旋转因子:pW N =J 2L W J=0,1,2,… ,1-L 2-1L2=M2×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=M2);确定一蝶形两输入数据距离B=1-L 2XL -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(3)循环层2:确定L 级的B=1-L 2个旋转因子;旋转因子指数p=J ×L -M 2,J=0~B-1; (4)循环层3:对于同一旋转因子,用于同一级L -M 2个蝶形运算中:k 的取值从J 到N-1,步长为L 2 (使用同一旋转因子的蝶形相距的距离) (5)完成一个蝶形运算。
开 始送入x (n ),MN =2M 倒 序L =1 , M J=0 , B - 1P =2M -L J k = J , N -1 , 2L pN 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)<Nxn=[xn,zeros(1,N-length(xn))]; endfor m=0:N/2-1;%旋转因子指数范围WN(m+1)=exp(-j*2*pi/N)^m;%计算旋转因子enddisp('输入到各存储单元的数据:'),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;endJ=J+K;enddisp('倒序后各存储单元的数据:'),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;endenddisp('本级运算后各存储单元的数据:'),disp(xn);end在主函数中调用myDitFFT(xn)函数实现DIT-FFT并和直接DFT运算结果做对比:xn=[0,1,2,3,4,5,6,7];myDitFFT(xn);调用fft函数运算的结果:1 至 7 列28.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 + 0.0000i -4.0000 - 1.6569i -4.0000 - 4.0000i8 列-4.0000 - 9.6569i调用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 列12.0000 + 0.0000i -4.0000 + 4.0000i -4.0000 + 0.0000i -4.0000 - 4.0000i 16.0000 + 0.0000i -4.0000 + 4.0000i -4.0000 + 0.0000i8 列-4.0000 - 4.0000i运算级次:3本级运算后各存储单元的数据:1 至 7 列28.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 + 0.0000i -4.0000 - 1.6569i -4.0000 - 4.0000i8 列-4.0000 - 9.6569i经对比可知DIT-FFT与直接DFT的运行结果完全相同。
四、总结经过验证可发现DIT-FFT较直接DFT运算有着明显的优势,我们可以将这个函数运用在多个领域以简化运算,例如计算离散时间序列的卷积或计算IDFT时都可以应用到DIT-FFT算法,我感受到数字信号处理中科学思想的魅力。
由于对设计思路的缺乏,我在设计程序时,在网络上查找了很多有关DIT-FFT的资料,经过学习他人的解决思路最后才整理出DIT-FFT的程序,在有些地方我自己理解的还不是很透彻,比如在实现数据倒序的程序我认为比较困难;当然即使自己想不到能学习一下别人的思路也是很好的,这个程序的代码量并不大,我自身的能力还很低,要在以后的学习中不断进步才能完成更加复杂的任务。
这次课程设计让我对快速傅里叶变换有了更多的了解,也认识到了科学计算方法的重要性,我感到很充实。
参考文献——百度百科;按时间抽取的基2FFT算法分析及MATLAB实现[J].电子技术,2011(2)数字信号处理(第四版)西安电子科技大学出版社高希全丁玉美编。