串行FFT递归算法(蝶式递归计算原理)求傅里叶变换

串行FFT递归算法（蝶式递归计算原理）求傅里叶变换

摘要

FFT，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）^2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2（N/2）^2=N+N^2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog(2)(N)次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。

关键字：FFT 蝶式计算傅里叶变换

目录

一．题目及要求 (1)

1.1题目 (1)

二．设计算法、算法原理 (1)

2.1算法原理与设计 (1)

2.2设计步骤 (2)

三．算法描述、设计流程 (4)

3.1算法描述 (4)

3.2流程图 (5)

四．源程序代码及运行结果 (7)

4.1源程序代码 (7)

4.2运行结果 (12)

五．算法分析、优缺点 (13)

5.1算法分析 (13)

5.2优缺点 (14)

六．总结 (15)

七．参考文献 (16)

一．题目及要求

1.1题目

对给定的)23,27,16,8,30,74,22,21(--=α，利用串行FFT 递归算法（蝶式递归计算原理）计算其傅里叶变换的结果。 1.2要求

利用串行递归与蝶式递归原理，对给定的向量求解傅里叶变换的结果。

二．设计算法、算法原理

2.1算法原理与设计

蝶式递归计算原理：令 为n/2次单位元根，则有 ， 将b 向量的偶数项 和奇数项 分别记 为 和 。 注意推导中反复使用： 。

图2.1 公式图形

)2//(2~n i e πω＝2~ωω

＝T n b b b ),...,,(220-T n b b b ),...,,(1

3

1

-T

n b b b ),...,,(1

10-'''T

n b b b ),...,,(1102

-''''''p p sn n n ωωωωω==-==+,1,1,1ln 2

/

2.2设计步骤

DFT a a a a a a b b b T n T n n n n 的是因此，向量),...,,(),...,,(111102202

2

2

--+-+++

DFT a a a a a a b b b T n T

n n n n n 的是因此，向量))(),...,(),((),...,,(111

1101312

2

2

2

---+----ω

ω

对于以上的分析可画出如图 2.2所示的离散傅里叶变换递归计算流图。图2.3就是一个按此递归方法计算的n=8的FFT 蝶式计算图。

()

()

()()1

,,1,0)(~)

(~)(~)(~)()

()()()(210

1112

22110111222411201

122

412112241201

222

2

2

22

2

2

2

2

22

2

2

-=+=++++++++=++++++++=+++++

++++==='∑∑-=+---++---++--++---=n k k

k l k n l l l n l

l l n l

l

l l

l

l

n k k lk l l l a a a a

a a a a a a a a a a a a a a a a a a a a a a a

b b n

n n

n n n n n n n n n n n n n n n ωωωωω

ωωωωωωωωω偶数时：()()()()()()()1

,,1,0)(~)

(~)(~)(~)()

()()()(210

1111222211011112222

411201

11212111222

41201

)12(11)12(1)12(1)12(12)

12(211

201

)12(122

2

2

2

2

2

2

2

22

2

2

2

22

2

2

2

2

2

2

-=-=-++-+-+-=-++-+-+-=----

++++=++++

++++===''∑∑-=+----++----++---+----+-++++-+-++-=++n k k

k k l k n l l l n l

l

l

n l

l

l

l

l

n l n l l l l l n k k k l l l l a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a b b n n n n n n n n n n n n n n n n n n n n n n n n n

n n ωωωωωωωωωωωωωωωωωωωωωωωωωω

ω

ω

ω

ω

ω奇数时：

FFT 的蝶式递归计算图（有计算原理推出）：

图2.2 递归计算流图

特别的，ｎ＝８的FFT 蝶式计算图（展开的）：

图2.3 蝶式计算图

按输入元素k a 展开，前面将输出序列之元素j b 按其偶下标（l j 2=）和（12+=l j ）

展开，导出 ∑-=++10

2

)(~n

n k k

k l k a a ω和∑-=+-10

2

)(~n

n k k

k k l k a a ωω递归计算式，按此构造出了如图1所示的FFT 递归计算流程图。事实上，我们也可以将输入序列之元素k a 按其偶下标（k 2）

a 0

a

1

a n -1

a n 2

-1

a n 2a

n 2

+1

和几下标（12+k ）展开，则导出另一种形式的FFT 递归计算式 ∑-==1

n k k kj j a w b 。

三．算法描述、设计流程

3.1算法描述

SISD 上的FFT 分治递归算法：

输入: a=(a 0,a 1,…,a n-1); 输出: B=(b 0,b 1,…,b n-1) Procedure RFFT(a,b) begin

if n=1 then b 0=a 0 else

(1)RFFT(a 0,a 2,…,a n-2, u 0,u 1,…,u n/2-1) (2)RFFT(a 1,a 3,…,a n-1, v 0,v 1,…,v n/2-1) (3)z=1

(4)for j=0 to n-1 do (4.1)b j =u j mod n/2+zv j mod n/2 (4.2)z=z ω

endfor endif end

注: (1)算法时间复杂度t(n)=2t(n/2)+O(n) t(n)=O(nlogn)

n=8的FFT蝶式计算图:

图3.1 FFT蝶式计算图n=6的FFT递归计算流程图:

图3.2 FFT递归计算流程图

3.2流程图

四．源程序代码及运行结果

4.1源程序代码

/************FFT***********/

#include //整个程序输入和输出利用同一个空间x[N]，节约空间#include

#include

#define N 1000 //定义输入或者输出空间的最大长度typedefstruct

{

double real;

doubleimg;

}complex; //定义复数型变量的结构体

void fft(); //快速傅里叶变换函数声明

void initW(); //计算W(0)~W(size_x-1)的值函数声明void change(); //码元位置倒置函数函数声明

void add(complex,complex,complex *); /*复数加法*/

void mul(complex,complex,complex *); /*复数乘法*/

void sub(complex,complex,complex *); /*复数减法*/

void divi(complex,complex,complex *); /*复数除法*/

void output(); /*输出结果*/

complex x[N],*W; /*输出序列的值*/

intsize_x=0; /*输入序列的长度，只限2的N次方*/

double PI; //pi的值

int main()

{

inti;

system("cls");

PI=atan(1)*4;

printf("Please input the size of x:\n"); /*输入序列的长度*/ scanf("%d",&size_x);

printf("Please input the data in x[N]:(such as:5 6)\n");

for(i=0;i

initW(); //计算W(0)~W(size_x-1)的值

fft(); //利用fft快速算法进行DFT变化output(); //顺序输出size_x个fft的结果

return 0;

}

/*进行基-2 FFT运算,蝶形算法。这个算法的思路就是，先把计算过程分为

log(size_x)/log(2)-1级（用i控制级数）；

然后把每一级蝶形单元分组（用j控制组的第一个元素起始下标）；

最后算出某一级某一组每一个蝶形单元（用k控制个数，共l个）。*/

voidfft()

{

inti=0,j=0,k=0,l=0;

complexup,down,product;

change(); //实现对码位的倒置

for(i=0;i

{

l=1<

for(j=0;j

{

for(k=0;k

{ //算出j组中第k个蝶形单元

mul(x[j+k+l],W[(size_x/2/l)*k],&product); /*size/2/l是该级W的相邻上标差，l

是该级该组取的W总个数*/

add(x[j+k],product,&up);

sub(x[j+k],product,&down);

x[j+k]=up; //up为蝶形单元右上方的值

x[j+k+l]=down; //down为蝶形单元右下方的值

}

}

}

}

void initW() //计算W的实现函数

{

inti;

W=(complex *)malloc(sizeof(complex) * size_x);

/*申请size_x个复数W的空间(这部申请的空间有点多，实际上只要申请size_x/2个即可)*/

for(i=0;i<(size_x/2);i++)

/*预先计算出size_x/2个W的值，存放,由于蝶形算法只需要前size_x/2个值即可*/ {

W[i].real=cos(2*PI/size_x*i); //计算W的实部

W[i].img=-1*sin(2*PI/size_x*i); //计算W的虚部

}

}

void change() //输入的码组码位倒置实现函数{

complex temp;

unsigned short i=0,j=0,k=0;

double t;

for(i=0;i

{

k=i;

j=0;

t=(log(size_x)/log(2));

while((t--)>0)

{

j=j<<1;

j|=(k & 1);

k=k>>1;

}

if(j>i)

{

temp=x[i];

x[i]=x[j];

x[j]=temp;

}

}

}

void output() //输出结果实现函数{

inti;

printf("The result are as follows\n");

for(i=0;i

{

printf("%.4f",x[i].real); //输出实部if(x[i].img>=0.0001) //如果虚部的值大于0.0001，输出+jx.img的形式printf("+j%.4f\n",x[i].img);

else if(fabs(x[i].img)<0.0001)

printf("\n");

else

printf("-j%.4f\n",fabs(x[i].img));

//如果虚部的值小于-0.0001，输出-jx.img的形式

}

}

void add(complex a,complexb,complex *c) //复数加法实现函数{

c->real = a.real + b.real; //复数实部相加

c->img = a.img + b.img; //复数虚部相加

}

void mul(complex a,complexb,complex *c) //复数乘法实现函数{

c->real = a.real*b.real - a.img*b.img; //获取相乘结果的实部

c->img = a.real*b.img + a.img*b.real; //获取相乘结果的虚部

}

void sub(complex a,complexb,complex *c) //复数减法实现函数{

c->real = a.real - b.real; //复数实部相减

c->img = a.img - b.img; //复数虚部相减

}

void divi(complex a,complexb,complex *c) //复数除法实现函数{

c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real+b.img*b.img);

//获取相除结果的实部c->img = (a.img*b.real - a.real*b.img) / (b.real*b.real+b.img*b.img);

//获取相除结果的虚部}

4.2运行结果 （1）处理器p=8：

图4.1当)8,7,6,5,4,3,2,1(=α时串行FFT 输出结果

（2）处理器p=8：

当)23,27,16,8,30,74,22,21(--=α时输出结果与计算结果相符如图4.2所示

图4.2运行图

五．算法分析、优缺点

5.1算法分析

（1）FFT算法的基本原理是把长序列的DFT逐次分解为较短序列的DFT。按照抽取方式的不同可分为DIT-FFT（按时间抽取）和DIF-FFT（按频率抽取）算法。按照蝶形运算的构成不同可分为基2、基4、基8以及任意因子（2n,n为大于1的整数），基2、基4算法较为常用。

（2）总体结构说明

输入数据为串行的数据流，故在第一级蝶形运算模块前加入串并转换模块，将串行数据流转换为并行的两列数据流以适应基2蝶形运算模块的输入信号要求。

由于每级蝶形运算一次处理的两个输入数据不能直接由前一级蝶形运算一次性输出，故在两个蝶形运算单元之间插入延时对齐模块，将前一级蝶形运算的结果（两列并行的数据流）作适当的延时并通过转接器对齐，形成后一级蝶形运算模块所需要的2列输入序列。

在最后一级蝶形运算后加入串并转换模块，将2列并行的数据流合成为1列。最后加入倒序模块将DIF-FFT得到的倒序输出序列整理为顺序输出。

旋转因子产生模块产生各级基2蝶形运算所需的旋转因子。由运算流图可以看出最后一级的旋转因子其实是1，故可省略最后一级蝶形运算单元中的旋转因子乘法器。因此用一个双口ROM将两组数据分别输出到第一级和第二级的蝶形运算单元即可。

基2蝶形运算模块由两个复数加法器和一个复数乘法器构成。旋转因子由ROM产生后，作为复数乘法器的输入之一，与前面复数加法器得到的结果相乘完成一次蝶形运算。为提高系统的运行速度可在蝶形运算单元中插入流水线寄存中间结果。

（3）蝶形运算单元如下所示：

图5.1运算单元

5.2优缺点

优点：结构清晰，可读性强，而且容易用数学归纳法来证明算法的正确性，因此它为设计算法、调试程序带来很大方便。

缺点：递归算法的运行效率较低，无论是耗费的计算时间还是占用的存储空间都比非递归算法要多。

六．总结

经过一周的课程设计，使我更加深刻的学习课本知识，又复习巩固了以前学过的知识。虽然这次课程是那么短暂的一周时间，但是我感觉到这些天我的收获远胜过我在课堂上所学，虽然这次任务主要上是设计，其实就是让你利用课本所学与实践相结合，是让我对所学算法和原来学的编程的熟悉和应用。

众所周知，通过课程设计能使我们了解到一些实际与理论之间的差异，对我们是很重要的。通过课程设计不仅可以巩固专业知识，为以后的工作打下了坚实的基础，而其还可以培养和熟练使用资料，运用工具书的能力，把我们所学的课本知识与实践结合起来，起到温故而知新的作用。课程设计诚然是一门专业课，给我很多专业知识以及专业技能上的提升，同时又是一门讲道课，一门设计课，给了我许多道，给了我很多思，给了我莫大的空间。同时，设计让我感触很深。使我对抽象的理论有了一些具体的认识。在课程设计过程中。我们要比较系统的了解串行算法FFT的设计中的每一个环节，包括从总体设计原则，以设计要求和设计题目的指导思想为中心，参照有关资料，有计划有头绪、有逻辑地把这次设计搞好！

总而言之，这次课程设计使我收获良多，以往的知识也捡起了不少。真挚的感谢我们的辅导老师和同学，在设计过程中，老师和同学们的帮助才使得我以顺利的完成这次设计。

七．参考文献

[1] 陈国良, 吴俊敏, 章锋等, 并行计算机体系结构[M], 北京: 高等教育出版社,

2002.

[2] 陈国良, 并行算法的设计与分析[M] 第三版, 北京: 高等教育出版社, 2009.

[3] 陈国良等, 并行算法实践[M], 北京: 高等教育出版社, 2003.

[4] 陈国良, 并行计算----结构·算法·编程[M], 第三版北京：高等教育出版社, 2011.