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

按时间抽取的基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蝶形运算流图2、 DIT-FFT算法的运算规律及编程思想1.原位计算 对N=点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M级运算后，原来存放输入序列数据的N个存储单元中可依次存放X(k)的N个值，这种原位(址)计算的方法可节省大量内存。2.旋转因子的变化规律N点DITFFT运算流图中，每个蝶形都要乘以旋转因子，p称为旋转因子的指数。例如N8 时各级的旋转因子： 第一级：L=1， 有1个旋转因子：= J=0 第二级：L=2，有2个旋转因子：= = J=0,1 第三级：L=3，有4个旋转因子：= = J=0,1,2,3对于N的一般情况，第L级共有个不同的旋转因子： = J=0,1,2, ,1 = N故： 按照上面两式可以确定第L级运算的旋转因子3、同一级中，同一旋转因子对应蝶形数目 第L级FFT运算中，同一旋转因子用在个蝶形中；4、 同一级中，蝶形运算使用相同旋转因子之间相隔的“距离” 第L级中，蝶距：D=；5、同一蝶形运算两输入数据的距离 在输入倒序，输出原序的FFT变换中，第L级的每一个蝶形的2个输入数据相距：B=。 6、码位颠倒 输入序列x(n)经过M级时域奇、偶抽选后，输出序列X(k)的顺序和输入序列的顺序关系为倒位关系。将十进制顺序数用I表示，与之对应的二进制是用IB表示，十进制倒序数用J表示，与之对应的二进制是用JB表示。十进制顺序数I增加1，相当于IB最低位加1且逢2向高位进1，即相当于JB最高位加1且逢2向低位进1。JB的变化规律反映到J的变化分为两种情况，若JB的最高位是0（JN/2），则直接由加1（JJ+N/2）得到下一个倒序值，若JB的最高位是1（JN/2），则要先将最高位变0（JJ-N/2），再在次高位加1（JJ+N/4），但次高位加1时，同样要判断0、1值，如果是0（JN/4）,则直接加1（JJ+N/4）,否则要先将次高位变0（JJ-N/4）再判断下一位，依次类推，直到完成最高位加1，逢2向右进位的运算。I=J时不需要交换，只对IJ时的情况进行数据交换即可，数据倒序程序框图如如2。7、蝶形运算的规律 序列经过时域抽选后，存入数组中，如果蝶形运算的两个输入数据相距B个点，应用原位计算，蝶形运算可表示成如下形式： XL-1(J)X L-1 (J+B)XL (J)= XL-1(J)+ WNp X L-1 (J+B)XL (J) = XL-1(J)-WNp X L-1 (J+B)p=J2M-L， J=0,1,2, ,2L-118、 DIT-FFT程序框图根据DIT-FFT原理和过程，DIT-FFT的完整程序框图如图2：(1)倒序:输入自然顺序序列x(n)，根据倒序规律，进行倒序处理；(2)循环层1:确定运算的级数，L=1M (N=)；确定一蝶形两输入数据距离B=(3)循环层2:确定L级的B=个旋转因子；旋转因子指数p=J，J=0B-1；(4)循环层3:对于同一旋转因子，用于同一级个蝶形运算中：k的取值从J到N-1，步长为 (使用同一旋转因子的蝶形相距的距离)(5)完成一个蝶形运算。 图2 数据倒序程序框图 图3 DIT-FFT的完整程序框图 3、 程序源代码设计函数myDitFFT(xn)完成一个序列的DIT-FFT运算：function y=myDitFFT(xn)M=nextpow2(length(xn);N=2M;disp(调用fft函数运算的结果:),fftxn=fft(xn,N);if length(xn)N xn=xn,zeros(1,N-length(xn); endfor 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=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:2L: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 列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.0000i 8 列 -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.0000i 8 列 -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.0000i 8 列 -4.0000 - 9.6569i经对比可知DIT-FFT与直接DFT的运行结果完全相同。4、 总结经过验证可发现DIT-FFT较直接DFT运算有着明显的优势，我们可以将这个函数运用在多个领域以简化运算，例如计算离散时间序列的卷积或计算IDFT时都可以应用到DIT-FFT算法，我感受到数字信号处理中科学思想的魅力。由于对设计思路的缺乏，我在设计程序时，在网络上查找了很多有关DIT-FFT的资料，经过学习他人的解决思路最后才整理出DIT-FFT的程序，在有些地方我自己理解的还不是很透彻，比如在实现数据倒序的程序我认为比较困难；当然即使自己想不到能学习一下别人的思路也是很好的，