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

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

1 DIT-FFT算法的基本原理

有限长序列x(n)的N点DFT定义为：X(k) x(n)W，式中WNk nN

n 0N 1 e j2 N，其整数次幂简称为旋转因子。

直接进行DFT运算大约需要2N2次三角函数计算、4N2次实数乘法计算和2N(2N 1)次实数加法计算，且需许多索引和寻址操作［2］。本文列出了直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在实际工作中并不实用。基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点时域序列。一个完整的8点DIT-FFT运算流程如图2所示［4］。图中的输入序列不再是顺序排列但有规律可循，数组A(存储地址)用于存放输入数据和每级运算的结果。

2 DIT-FFT算法的运算规律及编程思想

为了编写DIT-FFT算法的运算程序，首先要分析其运

算规律，总结编程思想并绘出程序框图。由图2可知，

DIT-FFT算法的运算过程很有规律。

2.1 原位计算

对N 2M点的FFT共进行M级运算，每级由N/2个

蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶

有用，且输出节点与输入节点在同一水平线上，这就意味

着每算完一个蝶后，所得数据可立即存入原输入数据所占

用的数组元素(存储单元)，这种原位(址)计算的方法可节

省大量内存。

2.2 蝶形运算 图1 DIT蝶形运算流图符号

图2 8点DIT-FFT运算流程

实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。对N 2M点的FFT共有M级运算，用L表示从左到右的运算级数(L=1，2，…，M )。第L级共有B 2L 1个不同指数的旋转因子，用R表示这些不同指数旋转因子从上到下的顺序(R=0，1，…，B-1)。第R个旋转因子的指数P 2M LR，旋转因子指数为P的第一个蝶的第一节点标号k从R开始，由于本级中旋转因子指数相同的蝶共有2M L个，且这些蝶的相邻间距为2L，故旋转因子指数为P的最后一个蝶的第一节点标号k为：(2M L 1) 2L R N 2L R，本级中各蝶的第二个节点与第一个节点都相距B点。应用原位计算，蝶形运算可表示成如下形式：

PAL(k) AL 1(k) WNAL 1(k B) PAL(k B) AL 1(k) WNAL 1(k B)

总结上述运算规律，可采用如下运算方法进行DIT-FFT

运算。首先读入数据，根据数据长度确定运

算级数M，运算总点数N 2M，不足补0处理。然后对读入数据进行

数据倒序操作。数据倒序后从第1级开始逐级进行，共进行M级运算。

在进行第L级运算时，先算出该级不同旋转因子的个数B 2L 1(也是该

级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到

R=B-1结束。每个R对应的旋转因子指数P 2M LR，旋转因子指数相

同的蝶从上往下依次逐个运算，各个蝶的第一节点标号k都是从R开始，

以2L为步长，到N 2L R(可简取极值N-2)结束。考虑到蝶形运算有两

个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，

输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据

完毕后才能覆盖。这样数据倒序后的运算可用三重循环程序实现。整个

运算流程如图3所示。

2.3 序列倒序

为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序输入，即

在运算前要先对输入的序列进行位序颠倒。如果总点数为N 2M的x(n)

的顺序数是用M位二进制数表示，则倒序数只需将顺序数的二进制位

倒置即可［5］，按照这一规律用硬件电路和汇编语言很容易产生倒序数。

但用MATLAB等高级语言实现倒序时，直接倒置二进制数位的方法不可

取，还须找出产生倒序的十进制规律。将十进制顺序数用I表示，与之

对应的二进制数用IB表示。十进制倒序数用J表示，与之对应的二进

制数用JB表示。JB是IB的位倒置结果，十进制顺序数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时不需要交换，当I ≠ J时需要交换数据。另外，为了避免再次调换前面已经调换过的一对数据，只对I<J的情况进行数据交换即可实现数据倒序操作。图3中数据倒序的程序框图如图4所示。 图3 DIT-FFT运算程序框图

3 MATLAB程序实现

MATLAB提供的fft函数是一个计算DFT的智能程序，能自动选择快速算法进行DFT运算，由于它是一个内建函数，用type命令看不到程序代码。为了体现编程思想，下面结合DIT-FFT程序框图，列出MATLAB环境下的实现程序。MATLAB的数组元素按序存储，可用下标寻访，但下标是从1开始的，所以在MATLAB程序中，寻访数组中的元素(数据)时，下标要在原序号上加1。旋转因子可按指数预先计算出来并存放到数组WN中，虽然占用了一些内存，但程序运行时可直接寻访调用，无需反复计算，可进一步提高运算速度。用MATLAB实现DIT-FFT算法的程序如下:

%基2DIT-FFT运算的MATLAB程序

clc;close all;clear;format compact;

%输入数据并计算常量

xn=[0,1,2,3,4,5,6,7];%可取任意序列

M=nextpow2(length(xn)), N=2^M,

for m=0:N/2-1;%旋转因子指数范围

WN(m+1)=exp(-j*2*pi/N)^m;%计算旋转因子

end

A=[xn,zeros(1,N-length(xn))]; %数据输入

disp('输入到各存储单元的数据:'),disp(A);

%数据倒序操作

J=0;%给倒序数赋初值

for I=0:N-1;%按序交换数据和算倒序数

if I<J;%条件判断及数据交换

T=A(I+1);A(I+1)=A(J+1);A(J+1)=T;

end

%算下一个倒序数

K=N/2;

while J>=K;

J=J-K;K=K/2;

end

J=J+K;

end

disp('倒序后各存储单元的数据:'),disp(A);

%分级按序依次进行蝶形运算

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=A(K+1)+A(K+B+1)*WN(P+1);

A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1);

A(K+1)=T;

end

end

disp('本级运算后各存储单元的数据:'),disp(A);

end

disp('输出各存储单元的数据:'),Xk=A,

disp(' …… 此处隐藏：2065字，全部文档内容请下载后查看。喜欢就下载吧 ……