并行计算MPI程序设计

*******************实践教学*******************兰州理工大学理学院2016年春季学期并行计算课程设计专业班级：2013级信息与计算科学姓名：学号：13550418指导教师：成绩：摘要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蝶式计算傅里叶变换目录摘要...............................................................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六、总结..........................................................19七、参考文献......................................................20一、题目及要求1.1题目对给定的α=（1,2，4，3，8，6，7，2），利用串行FFT递归算法(蝶式递归计算原理)计算其傅里叶变换的结果1.2要求利用串行递归与蝶式递归原理，对给定的向量求解傅里叶变换的结果二、算法设计与算法原理2.1算法原理与设计令为n/2次单位元根，则有.将b向量的偶数项和奇数项分别记为和注意推导中反复使用图2.1)2//(2~nie＝2~＝Tnbbb),...,,(220Tnbbb),...,,(131Tnbbb),...,,(1102Tnbbb),...,,(1102,,1,1,1ln2/ppsnnn2.2设计求解步骤DFTaaaaaabbblaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaabbTnTnnkkkklknlllnlllnlllllnlnlllllnkkklllnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn的是因此，向量奇数时：))(),...,(),((),...,,(1,,1,0)(~)(~)(~)(~)()()()()(11111013121011112222110111122224112011121211122241201)12(11)12(1)12(1)12(12)12(2112010)12(122222222222222222222222222222222DFTaaaaaabbblaaaaaaaaaaaaaaaaaaaaaaaaaaabbTnTnnkkklknlllnlllnllllllnkklkllnnnnnnnnnnnnnnnnnnnnn的是向量因此，偶数时：),...,,(),...,,(1,,1,0)(~)(~)(~)(~)()()()()(11110220210111222110111222411201122412112241201022222222222222222222222三、算法描述与算法流程3.1算法描述n=8的FFT蝶式计算图图3.1图3.2FFT递归计算流程图3.2流程图飞是否是否是否图3.3输入序列长度size_x输入序列对应值（例如5+j3，输入53）计算出前size_x/2个exp(-j*2π*k/size_x)个值，即W的值级数i=？输出fft结果序列开始结束计算出该级需要的W的个数l该级该组起始下标j=？级数i加1组起始下标加2*l该级该组元素序数k=？X[j+k]X[j+k]lX[j+k+l]*W[(size_x/2/l)*k]X[j+k+l]-1K加1四、源程序代码与运行结果4.1源程序/************FFT***********///整个程序输入和输出利用同一个空间x[N]，节约空间#includestdio.h#includemath.h#includestdlib.h#defineN1000//定义输入或者输出空间的最大长度typedefstruct{doublereal;doubleimg;}complex;//定义复数型变量的结构体voidfft();//快速傅里叶变换函数声明voidinitW();//计算W(0)~W(size_x-1)的值函数声明voidchange();//码元位置倒置函数函数声明voidadd(complex,complex,complex*);/*复数加法*/voidmul(complex,complex,complex*);/*复数乘法*/voidsub(complex,complex,complex*);/*复数减法*/voiddivi(complex,complex,complex*);/*复数除法*/voidoutput();/*输出结果*/complexx[N],*W;/*输出序列的值*/intsize_x=0;/*输入序列的长度，只限2的N次方*/doublePI;//pi的值intmain(){inti;system(cls);PI=atan(1)*4;printf(Pleaseinputthesizeofx:\n);/*输入序列的长度*/scanf(%d,&size_x);printf(Pleaseinputthedatainx[N]:(suchas:56)\n);/*输入序列对应的值*/for(i=0;isize_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的结果return0;}/*进行基-2FFT运算,蝶形算法。这个算法的思路就是，先把计算过程分为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;ilog(size_x)/log(2);i++)//循环算出fft的结果{l=1i;for(j=0;jsize_x;j+=2*l)//算出第m=i级的结果【i从0到（log(size_x)/log(2)）-1】{for(k=0;kl;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为蝶形单元右下方的值}}}}voidinitW()//计算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的虚部}}voidchange()//输入的码组码位倒置实现函数{complextemp;unsignedshorti=0,j=0,k=0;double