并行计算课程设计任务书分析

*******************实践教学*******************兰州理工大学理学院2016年春季学期并行计算课程设计专业班级：2013级信息与计算科学2班姓名：学号：13540216指导教师：成绩：基于串行FFT蝶式递归计算摘要本文主要设计快速傅氏变换两种经典的串行算法之一的递归算法（蝶式）。

它是利用分治的思想来推导递归的FFT计算算法，将b元素或a元素之下标分别按偶下标与奇下标展开而推导出的递归式。

算法主要过程如下：首先主进程对输入的元素按其下标进行奇偶划分，然后分配给两个处理器，分别在各自的处理器在进行同样的划分、并计算其结果，将结果带入到上一层再进行计算，最后将结果带入到主进程合并进行计算，输出结果。

关键字：蝶式计算傅里叶变换递归算法目录一．题目及要求 (1)1.1设计题目 (1)1.2设计参数 (1)1.3测试实例 (1)二．设计算法及要求 (1)2.1 算法设计原理 (1)2.2 算法设计 (2)三．算法描述与设计流程 (3)3.1算法描述 (3)3.2 流程图 (4)四．原程序代码与运行结果 (5)4.1 源程序 (5)4.2 运行结果 (10)五．算法分析及其优缺点 (11)5.1 算法分析 (11)5.2 优缺点 (12)六．总结 (13)七．参考文献 (14)一．题目及要求1.1设计题目串行FFT 递归算法（蝶式递归计算原理）求傅里叶变换1.2设计参数蝶式递归计算原理、算法性能和误差分析1.3测试实例对给定的)3,7,6,78,30,4,32,201(---=α，利用串行FFT 递归算法 （蝶式递归计算原理）计算其傅里叶变换的结果。

二．设计算法及要求2.1 算法设计原理令 为n/2次单位元根，则有 . 将b 向量的偶数项 和奇数项 分别记为 和 注意推导中反复使用()()()()DFTa a a a a ab b b 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 T n T n n k kk l k n l l l n ll l n ll l ll l n k klk l l n n n nn n n n n n n n n n n n n n n n n 的是向量因此，：),...,,(),...,,(1,,1,0)(~)(~)(~)(~)()()()()(偶数时1111022021011122211011122241120112241211224120122222222222222222--+--=+---++---++--++---=+++-=+=++++++++=++++++++=+++++++++==='∑∑ ωωωωωωωωωωωωωωTn b b b ),...,,(220-Tn b b b ),...,,(131-T n b b b ),...,,(1102-'''T n b b b ),...,,(1102-'''''')2//(2＝~n i e πω2~ωω＝,,1,1,1ln 2/p p sn n nωωωωω==-==+()()()()()()()DFTa a a a a ab b b l a a a a a a a a a a a aa a a a a a a a a aa a a a a a a a a a ab b Tn T n n k kk 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 kk l l l n nn n nn n n n n n n nnn n n n nn n n nnn n nn nn n 的是因此，向量奇数时：))(),...,(),((),...,,(1,,1,0)(~)(~)(~)(~)()()()()(11111013121011112222110111122224112011121211122241201)12(11)12(1)12(1)12(12)12(211201)12(1222222222222222222---+--=+----++----++---+----+-++++-+-++-=++----=-=-++-+-+-=-++-+-+-=----++++=++++++++===''∑∑ωωωωωωωωωωωωωωωωωωωωωωωωωωωωωωωωω2.2 算法设计对于以上的分析可画出如图 1所示的离散傅里叶变换递归计算流图。

图1 n=8的FFT 蝶式计算图三．算法描述与设计流程3.1算法描述3.2图四．原程序代码与运行结果4.1 源程序/************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间#include <stdio.h>#include <math.h>#include <stdlib.h>#define N 1000 //定义输入或者输出空间的最大长度typedef struct{double real;double img;}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; /*输出序列的值*/int size_x=0; /*输入序列的长度，只限2的N次方*/ double PI; //pi的值int main(){int i;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<size_x;i++)scanf("%lf %lf",&x[i].real,&x[i].img);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个）。

*/void fft(){int i=0,j=0,k=0,l=0;complex up,down,product;change(); //实现对码位的倒置for(i=0;i<log(size_x)/log(2);i++) //循环算出fft的结果{l=1<<i;for(j=0;j<size_x;j+=2*l)//算出第m=i级的结果【i从0到（log(size_x)/log(2)）-1】{for(k=0;k<l;k++) //算出第i级内j组蝶形单元的结果{ //算出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的实现函数{int i;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<size_x;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() //输出结果实现函数{int i;printf("The result are as follows\n");for(i=0;i<size_x;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) //如果虚部的绝对值小于0.0001，无需输出printf("\n");elseprintf("-j%.4f\n",fabs(x[i].img));//如果虚部的值小于-0.0001，输出-jx.img的形式}}void add(complex a,complex b,complex *c) //复数加法实现函数{c->real = a.real + b.real; //复数实部相加c->img = a.img + b.img; //复数虚部相加}void mul(complex a,complex b,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,complex b,complex *c) //复数减法实现函数{c->real = a.real - b.real; //复数实部相减c->img = a.img - b.img; //复数虚部相减}void divi(complex a,complex b,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 运行结果图3 运行结果图五．算法分析及其优缺点5.1 算法分析（1）FFT算法的基本原理是把长序列的DFT逐次分解为较短序列的DFT。