并行計算MPI程序設計

1、*實踐教學*蘭州理工大學理學院 2016年春季學期并行計算課程設計專業(yè)班級： 2013級信息與計算科學姓名： 學號：13550418 指導教師： 成績：摘要FFT，即為快速傅氏變換，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機系統(tǒng)或者說數(shù)字系統(tǒng)中應用離散傅立葉變換，可以說是進了一大步。   設x(n)為N項的復數(shù)序列，由DFT變換，任一X（m）的計算都需要N次復數(shù)乘法和N-1次復數(shù)加法，而一次復數(shù)乘法等于四次實數(shù)乘法和兩次實數(shù)加法，一次復數(shù)加法等于兩次

2、實數(shù)加法，即使把一次復數(shù)乘法和一次復數(shù)加法定義成一次“運算”（四次實數(shù)乘法和四次實數(shù)加法），那么求出N項復數(shù)序列的X（m）,即N點DFT變換大約就需要N2次運算。當N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列（設N=2k,k為正整數(shù)），分為兩個N/2項的子序列，每個N/2點DFT變換需要（N/2）2次運算，再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以后，總的運算次數(shù)就變成N+2（N/2）2=N+N2/2。繼續(xù)上面的例子，N=1024時，總的運算次數(shù)就變成了525312次，節(jié)省了大約50%的運算

3、量。而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那么N點的DFT變換就只需要Nlog(2)(N)次的運算，N在1024點時，運算量僅有10240次，是先前的直接算法的1%，點數(shù)越多，運算量的節(jié)約就越大，這就是FFT的優(yōu)越性關鍵字：FFT 蝶式計算 傅里葉變換目錄摘要2目錄3一、題目及要求41.1題目41.2要求4二、算法設計與算法原理52.1 算法原理與設計52.2設計求解步驟6三、算法描述與算法流程73.1算法描述73.2 流程圖9四、源程序代碼與運行結果104.1源程序104.2運行結果16五、算法分析及其優(yōu)缺點175.1 算法分析175.2優(yōu)缺點18

4、六、總結19七、參考文獻20一、題目及要求1.1題目對給定的=（1,2，4，3，8，6，7，2），利用串行FFT遞歸算法(蝶式遞歸計算原理)計算其傅里葉變換的結果1.2要求利用串行遞歸與蝶式遞歸原理，對給定的向量求解傅里葉變換的結果二、算法設計與算法原理2.1 算法原理與設計 令 為n/2次單位元根，則有 . 將b向量的偶數(shù)項和奇數(shù)項分別記為 和 注意推導中反復使用圖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飛級數(shù)i>=？級數(shù)i加1 是 輸出fft結果序列結束 否 該級該組起始下標j>=？計算出該級需要的W的個數(shù)l 是 否組起始下標加2*l該級該組元素序數(shù)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，節(jié)約空間 #include <stdio.h> #include <math.h> #include <stdlib.h> #define

6、 N 1000 /定義輸入或者輸出空間的最大長度typedefstruct double real;doubleimg; complex; /定義復數(shù)型變量的結構體 void fft(); /快速傅里葉變換函數(shù)聲明 void initW(); /計算W(0)W(size_x-1)的值函數(shù)聲明 void change(); /碼元位置倒置函數(shù)函數(shù)聲明 void add(complex,complex,complex *); /*復數(shù)加法*/ void mul(complex,complex,complex *); /*復數(shù)乘法*/ void sub(complex,complex,complex

7、 *); /*復數(shù)減法*/ void divi(complex,complex,complex *); /*復數(shù)除法*/ 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控制級數(shù)）； 然后把每一級蝶形單元分組（用j控制組的第一個元素起始下標）； 最后算出某一級某一組每一個蝶形單元（用k控制個數(shù)，共l個）。*/voidfft() inti=0,j=0,k=0,l=0; complexup,down,product; change(); /實現(xiàn)對碼位的倒置 for(i=0;i<log(size_x)/log(2);i+) /循環(huán)算出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級內(nèi)j組蝶形單元的結果 /算出j組中第k個蝶形單元mul(xj+k+l,W(size_x/2/l)*k,&product); /*size/2/l是該級W的相鄰上標差，l是該級該組取的W總個數(shù)*/add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; /up為蝶形單元右上方的值 xj+k+l=down; /down為蝶形單元右下方的值 void initW() /計算W的實現(xiàn)函數(shù) inti; W=(complex *)malloc(sizeof(complex) * size_x); /*

11、申請size_x個復數(shù)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() /輸入的碼組碼位倒置實現(xiàn)函數(shù) 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() /輸出結果實現(xiàn)函數(shù) 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) /復數(shù)加法實現(xiàn)函數(shù) c->real = a.real + b.real; /復數(shù)實部相

14、加 c->img = a.img + b.img; /復數(shù)虛部相加 void mul(complex a,complexb,complex *c) /復數(shù)乘法實現(xiàn)函數(shù) 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) /復數(shù)減法實現(xiàn)函數(shù) c->real = a.real - b.real; /復數(shù)實部相減 c->img = a.img -

15、b.img; /復數(shù)虛部相減 void divi(complex a,complexb,complex *c) /復數(shù)除法實現(xiàn)函數(shù) 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五、算法分析及其優(yōu)缺點5.1 算法分析（1）FFT算法的

16、基本原理是把長序列的DFT逐次分解為較短序列的DFT。按照抽取方式的不同可分為DIT-FFT（按時間抽?。┖虳IF-FFT（按頻率抽?。┧惴?。按照蝶形運算的構成不同可分為基2、基4、基8以及任意因子（2n,n為大于1的整數(shù)），基2、基4算法較為常用。（2）總體結構說明 輸入數(shù)據(jù)為串行的數(shù)據(jù)流，故在第一級蝶形運算模塊前加入串并轉換模塊，將串行數(shù)據(jù)流轉換為并行的兩列數(shù)據(jù)流以適應基2蝶形運算模塊的輸入信號要求。 由于每級蝶形運算一次處理的兩個輸入數(shù)據(jù)不能直接由前一級蝶形運算一次性輸出，故在兩個蝶形運算單元之間插入延時對齊模塊，將前一級蝶形運算的結果（兩列并行的數(shù)據(jù)流）作適當?shù)难訒r

17、并通過轉接器對齊，形成后一級蝶形運算模塊所需要的2列輸入序列。 在最后一級蝶形運算后加入串并轉換模塊，將2列并行的數(shù)據(jù)流合成為1列。最后加入倒序模塊將DIF-FFT得到的倒序輸出序列整理為順序輸出。 旋轉因子產(chǎn)生模塊產(chǎn)生各級基2蝶形運算所需的旋轉因子。由運算流圖可以看出最后一級的旋轉因子其實是1，故可省略最后一級蝶形運算單元中的旋轉因子乘法器。因此用一個雙口ROM將兩組數(shù)據(jù)分別輸出到第一級和第二級的蝶形運算單元即可。 基2蝶形運算模塊由兩個復數(shù)加法器和一個復數(shù)乘法器構成。旋轉因子由ROM產(chǎn)生后，作為復數(shù)乘法器的輸入之一，與前面復數(shù)加法器得到的結果相乘完成一次蝶形運

18、算。為提高系統(tǒng)的運行速度可在蝶形運算單元中插入流水線寄存中間結果（3）蝶形運算單元如下所示：5.2優(yōu)缺點（1）優(yōu)點：結構清晰，可讀性強，而且容易用數(shù)學歸納法來證明算法的正確性，因此它為設計算法、調(diào)試程序帶來很大方便。 （2）缺點：遞歸算法的運行效率較低，無論是耗費的計算時間還是占用的存儲空間都比非遞歸算法要多。六、總結通過這次課程設計，讓我更加深刻了解課本知識，和以往對知識的疏忽得以補充。雖然這次課程是那么短暫的1周時間，我感覺到這些天我的所學勝過我這一學期所學，這次任務原則上是設計，其實就是一次大的作業(yè)，是讓我對課本知識的鞏固和對基本公式的熟悉和應用課程設計是一個重要的教學環(huán)節(jié)，通過課程設計使我們了解到一些實際與理論之間的差異。通過課程設計不僅可以鞏固專業(yè)知識，為以后的工作打下了堅實的基礎，而其還可以培養(yǎng)和熟練使用資料，運用工具書的能力，把我們所學的課本知識與實踐結合起來，起到溫故而知新的作用。課程設計誠然是一門專業(yè)課，給我很多專業(yè)知識以及專業(yè)技能上的提升，同時又是一門講道課，一門