并行计算MPI程序设计

1、*实践教学*兰州理工大学理学院 2016年春季学期并行计算课程设计专业班级： 2013级信息与计算科学姓名： 学号：13550418 指导教师： 成绩：摘要FFT，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。   设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次

2、实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。当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+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算

3、量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog(2)(N)次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性关键字：FFT 蝶式计算 傅里叶变换目录摘要2目录3一、题目及要求41.1题目41.2要求4二、算法设计与算法原理52.1 算法原理与设计52.2设计求解步骤6三、算法描述与算法流程73.1算法描述73.2 流程图9四、源程序代码与运行结果104.1源程序104.2运行结果16五、算法分析及其优缺点175.1 算法分析175.2优缺点18

4、六、总结19七、参考文献20一、题目及要求1.1题目对给定的=（1,2，4，3，8，6，7，2），利用串行FFT递归算法(蝶式递归计算原理)计算其傅里叶变换的结果1.2要求利用串行递归与蝶式递归原理，对给定的向量求解傅里叶变换的结果二、算法设计与算法原理2.1 算法原理与设计 令 为n/2次单位元根，则有 . 将b向量的偶数项和奇数项分别记为 和 注意推导中反复使用图2.12.2设计求解步骤三、算法描述与算法流程3.1算法描述n=8的FFT蝶式计算图图3.1图3.2 FFT递归计算流程图3.2 流程图开始计算出前size_x/2个exp(-j*2*k/size_x)个值，即W的值输入序列对应值

5、（例如5+j3，输入5 3）输入序列长度size_x飞级数i>=？级数i加1 是 输出fft结果序列结束 否 该级该组起始下标j>=？计算出该级需要的W的个数l 是 否组起始下标加2*l该级该组元素序数k>=？K加1 是Xj+k Xj+klXj+k+l*W(size_x/2/l)*k Xj+k+l -1 否图3.3四、源程序代码与运行结果4.1源程序/*FFT*/ /整个程序输入和输出利用同一个空间xN，节约空间 #include <stdio.h> #include <math.h> #include <stdlib.h> #define

6、 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

7、 *); /*复数减法*/ void divi(complex,complex,complex *); /*复数除法*/ void output(); /*输出结果*/ complex xN,*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",&

8、size_x);printf("Please input the data in xN:(such as:5 6)n"); /*输入序列对应的值*/for(i=0;i<size_x;i+)scanf("%lf %lf",&xi.real,&xi.img);initW(); /计算W(0)W(size_x-1)的值 fft(); /利用fft快速算法进行DFT变化 output(); /顺序输出size_x个fft的结果return 0; /*进行基-2 FFT运算,蝶形算法。这个算法的思路就是， 先把计算过程分为log(size_x

9、)/log(2)-1级（用i控制级数）； 然后把每一级蝶形单元分组（用j控制组的第一个元素起始下标）； 最后算出某一级某一组每一个蝶形单元（用k控制个数，共l个）。*/voidfft() inti=0,j=0,k=0,l=0; complexup,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<

10、;l;k+) /算出第i级内j组蝶形单元的结果 /算出j组中第k个蝶形单元mul(xj+k+l,W(size_x/2/l)*k,&product); /*size/2/l是该级W的相邻上标差，l是该级该组取的W总个数*/add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; /up为蝶形单元右上方的值 xj+k+l=down; /down为蝶形单元右下方的值 void initW() /计算W的实现函数 inti; W=(complex *)malloc(sizeof(complex) * size_x); /*

11、申请size_x个复数W的空间(这部申请的空间有点多，实际上只要申请size_x/2个即可)*/ for(i=0;i<(size_x/2);i+) /*预先计算出size_x/2个W的值，存放,由于蝶形算法只需要前size_x/2个值即可*/ Wi.real=cos(2*PI/size_x*i); /计算W的实部 Wi.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;

12、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=xi; xi=xj; xj=temp; void output() /输出结果实现函数 inti; printf("The result are as followsn"); for(i=0;i<size_x;i+) printf("%.4f",xi.real); /输出实部 if(xi.img>=0.0001) /如果

13、虚部的值大于0.0001，输出+jx.img的形式printf("+j%.4fn",xi.img); else if(fabs(xi.img)<0.0001) /如果虚部的绝对值小于0.0001，无需输出 printf("n"); elseprintf("-j%.4fn",fabs(xi.img); /如果虚部的值小于-0.0001，输出-jx.img的形式 void add(complex a,complexb,complex *c) /复数加法实现函数 c->real = a.real + b.real; /复数实部相

14、加 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 -

15、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=4图4.2.1（2）处理器p=2图4.2.2五、算法分析及其优缺点5.1 算法分析（1）FFT算法的

16、基本原理是把长序列的DFT逐次分解为较短序列的DFT。按照抽取方式的不同可分为DIT-FFT（按时间抽取）和DIF-FFT（按频率抽取）算法。按照蝶形运算的构成不同可分为基2、基4、基8以及任意因子（2n,n为大于1的整数），基2、基4算法较为常用。（2）总体结构说明 输入数据为串行的数据流，故在第一级蝶形运算模块前加入串并转换模块，将串行数据流转换为并行的两列数据流以适应基2蝶形运算模块的输入信号要求。 由于每级蝶形运算一次处理的两个输入数据不能直接由前一级蝶形运算一次性输出，故在两个蝶形运算单元之间插入延时对齐模块，将前一级蝶形运算的结果（两列并行的数据流）作适当的延时

17、并通过转接器对齐，形成后一级蝶形运算模块所需要的2列输入序列。 在最后一级蝶形运算后加入串并转换模块，将2列并行的数据流合成为1列。最后加入倒序模块将DIF-FFT得到的倒序输出序列整理为顺序输出。 旋转因子产生模块产生各级基2蝶形运算所需的旋转因子。由运算流图可以看出最后一级的旋转因子其实是1，故可省略最后一级蝶形运算单元中的旋转因子乘法器。因此用一个双口ROM将两组数据分别输出到第一级和第二级的蝶形运算单元即可。 基2蝶形运算模块由两个复数加法器和一个复数乘法器构成。旋转因子由ROM产生后，作为复数乘法器的输入之一，与前面复数加法器得到的结果相乘完成一次蝶形运

18、算。为提高系统的运行速度可在蝶形运算单元中插入流水线寄存中间结果（3）蝶形运算单元如下所示：5.2优缺点（1）优点：结构清晰，可读性强，而且容易用数学归纳法来证明算法的正确性，因此它为设计算法、调试程序带来很大方便。 （2）缺点：递归算法的运行效率较低，无论是耗费的计算时间还是占用的存储空间都比非递归算法要多。六、总结通过这次课程设计，让我更加深刻了解课本知识，和以往对知识的疏忽得以补充。虽然这次课程是那么短暂的1周时间，我感觉到这些天我的所学胜过我这一学期所学，这次任务原则上是设计，其实就是一次大的作业，是让我对课本知识的巩固和对基本公式的熟悉和应用课程设计是一个重要的教学环节，通过课程设计使我们了解到一些实际与理论之间的差异。通过课程设计不仅可以巩固专业知识，为以后的工作打下了坚实的基础，而其还可以培养和熟练使用资料，运用工具书的能力，把我们所学的课本知识与实践结合起来，起到温故而知新的作用。课程设计诚然是一门专业课，给我很多专业知识以及专业技能上的提升，同时又是一门讲道课，一门