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

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

二、

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

1.原位计算

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

2.旋转因子的变化规律

N点DIT―FFT运算流图中，每个蝶形都要乘以旋转因子p W

N

，p称为旋转因子的指数。例如N＝8 ＝32时各级的旋转因子：

第一级：L=1，有1个旋转因子：p W

N =J

/4

W

N

=J

2L

W

J=0

第二级：L=2，有2个旋转因子：p W

N =J

/2

W

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＝M2的一般情况，第L级共有1-L2个不同的旋转因子：

p W

N =J

2L

W J=0,1,2,… ,1-L2－1

L

2=M2×M-L2= N·M-L2

故：按照上面两式可以确定第L级运算的旋转因子

3、同一级中，同一旋转因子对应蝶形数目

第L级FFT运算中，同一旋转因子用在L-M2个蝶形中；

4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离”

第L级中，蝶距：D=L2；

5、同一蝶形运算两输入数据的距离

在输入倒序，输出原序的FFT变换中，第L 级的每一个蝶形的2个输入数据相距：B=1-L2。

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

7、蝶形运算的规律

序列经过时域抽选后，存入数组中，如果蝶形运算的两个输入数据相距B 个点，应用原位计算，蝶形运算可表示成如下形式：

8、 DIT-FFT 程序框图

根据DIT-FFT 原理和过程，DIT-FFT 的完整程序框图如图2：

(1)倒序:输入自然顺序序列x(n)，根据倒序规律，进行倒序处理；

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 W N

p=J ×2M-L ， J=0,1,2,… ,2L-1－1

(2)循环层1:确定运算的级数，L=1~M (N=M2)；确定一蝶形两输入数据距离B=1-L2

(3)循环层2:确定L级的B=1-L2个旋转因子；旋转因子指数p=J×L-M2，J=0~B-1；

(4)循环层3:对于同一旋转因子，用于同一级L-M2个蝶形运算中：k的取值从J到N-1，步长为L2 (使用同一旋转因子的蝶形相距的距离)

(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)

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

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;%分级计算