并行FFT实现

并行FFT算法的分析与实现

目录

1 FFT算法分析与实现 ............................................................................................................. 2 1.1 FFT算法的分析 ........................................................................................................... 2 1.1.1 FFT算法的意义 ....................................................................................................... 2 1.2 DFT算法和FFT算法在计算量上的对比 ...................................................................... 2 1.2.1 DFT正变换 ............................................................................................................... 2 1.2.2 变换后以2为底的FFT算法 ................................................................................. 2 2 基于二时分FFT算法的蝶式运算定理................................................................................ 3 2.1 并行可行性分析 ............................................................................................................ 3 2.2 难点分析 ........................................................................................................................ 4 2.3 FFT算法的实现 ........................................................................................................... 4 2.3.1 计算旋转因子 ......................................................................................................... 4 2.3.2 创建线程 ................................................................................................................. 5 2.3.3 关闭线程 ................................................................................................................. 5 3 蝶形变换 ............................................................................................................................... 5 4 实验结果分析 ....................................................................................................................... 6 4.1 FFT算法效率分析 .......................................................................................................... 6 4.1.1 奇偶两部分的算法效率分析 ................................................................................. 6 4.1.2 蝶形变换内部多线程效率分析 ............................................................................. 8 5 原程序 ................................................................................................................................... 9

1 FFT算法分析与实现

1.1 FFT算法的分析

1.1.1 FFT算法的意义

在离散傅里叶变换中，若要处理的数据量很大，工作量也变得非常巨大。正是由于这个原因，使得离散傅里叶变换的应用范围在过去很长的时间里受到了严格的。而快速傅里叶变换有效的解决了运算量的问题，使得它的理论与方法在数理方程、线性系统分析、信号处理与仿真等很多学科领域都有着广泛的应用。由于计算机只能处理有限长度的离散的序列，所以真正在计算机上运算的是一种离散傅里叶变换。若把此算法应用到并行程序设计中，它的效率将会有进一步的提高，将能得到更广泛的应用。

FFT算法有多个不同的变形：基二时分FFT算法、基二频分FFT算法、基四时分FFT算法、基四频分FFT算法。其中基二时分FFT算法是最重要的一种，它是J.W.Cooley和J.W.Tukey在1965年发表的一篇论文中提出。

1.2 DFT算法和FFT算法在计算量上的对比

1.2.1 DFT正变换

(1-1)

公式(1-1)可以看出，对每个n，计算X(n)须作N次复数乘法和N-1次复数加法，要全部计算完成需要N2次复数乘法和N*(N-1)次加法运算，当N非常大时，计算量将非常大。例如N=1024，则要运算2*(1024*1023)次。

1.2.2 变换后以2为底的FFT算法

由式(1-2)可以看出，共有r重和式，每重和式N个方程，每个方程仅需一次乘法运算，因此FFT算法的总计算量为Nr乘法和Nr次加法。如果在考虑WNK=-WNk+N/2，乘法的运算量还可以减少一半。如果在多核环境下用多线程技术让算法很好的并行运行，那么并行部分运

行的效率有可能再提升一倍。

（1-2）

由公式(1-1)和(1-2)可知：FFT的时间复杂度为O(N2)，DFT算法的时间复杂度为O(Nlog2N)，当N很大时，FFT算法的效率可以明显的显示出来。

2 基于二时分FFT算法的蝶式运算定理

设X(k)=DFT[x(n)](0<=N,K2.1 并行可行性分析

由定理可知，基二时分FFT算法在开始时就可以将要运算的数据分成对等的两部分，每部分都可以运算，这符合并行程序上的数据分解；且FFT运算内部不用建堆，运算后的数据可以替换先前的数据，这同样也满足了多线程的特点。若是必须要在线程内部建堆，那么当堆到一定程度时，就会出现内存读写错误。而此算法可以很好的避开线程内部建堆。即本算法适用于并行计算。

2.2 难点分析

本算法的难点有两处，一是数据分解的问题，一般情况下可以将数据平均分为两部分，但是基二时分FFT算法不行，因为此算法的后一步会用到上一步的结果，如果直接分成两部分得出的结果一定错误。根据算法本身的特点，再结合并行程序的特性，将数据以序列的形式分为奇偶两部分，将这两部分分别以基二时分FFT算法进行运算，得出的结果再根据定理进行汇总，以得出最终正确的结果。第二个难点在于基二时分FFT算法本身，算法中的蝶形变换本身就是一个三重循环体，每次内部循环的结果都会在下次循环被用到。图1-11为八点蝶形变换情况：

图1-1 八点蝶形变换图

2.3 FFT算法的实现

2.3.1 计算旋转因子

旋转因子的计算应分为两部分，一是先计算以数据量为整体数据量的一半的旋转因子，二是计算整体数据量的旋转因子的头半部分。

依照数据的顺序将数据分解成为奇偶两部分，且这两部分的数据量相等，因为FFT的数据量要求是2的n次方(n为整数)。

2.3.2 创建线程

创建两个线程，每个线程都以FFT的方法运行上面的一个数据块，因这两个数据块的大小相同，所以这两个线程运行的时间也基本相同，这样可以减少线程间等待的时间，创建线程函数如下：

HANDLE hThread[2]；

hThread[0]=CreateThread(NULL，0，FtmlProc，thread[O]，0,NULL)； hThread[1]=CreateThread(NULL，0，FunlProc,thread[1]，0,NULL)；

一个线程运行结束后，必须等待另外的一个线程运行结束，只有当两个线程都运行结束后才能进行下一步，阻塞等待如下：

WaitForMultipleObjects(NUM_THREADS，handle，true，INFINITE)；

2.3.3 关闭线程

关闭创建的线程以释放其所占的系统资源，关闭线程句柄如下： CloseHandle(Thread[0])； CloseHandle(Thread[1])；

将运行后的数据合并，当数据量很大时同样可以将合并部分分为多线程来进行。

3 蝶形变换

作为对比，先不考虑算法的特性，直接看蝶形变换的三层循环，将第三层循

环体分解为多个线程，若按照这个思路，整个程序的运行过程如下：先经过一次循环变换将数据反序输入，确保数据正序输出，再将数据输入至UFFT算法的核心(三层循环体)，由于上两层的循环需要用到以前的结果，不能出来，所以只能对第三层循环体进行线程划分。

第一层循环为： For(int pk=0；pk第三层循环将其分为两个线程，首先得出第三层循环的两个界限为：实现

过程与以奇偶形式划分数据集的FFT方式同。

4 实验结果分析

4.1 FFT算法效率分析

4.1.1 奇偶两部分的算法效率分析

分别使用并行傅立叶变换算法和串行傅立叶变换算法对不同个数的数据进行处理，比较两种算法的运行时间。实验结果如表4-1所示，其中N为输入数据总数，T1为FFT串行程序运行时间，T2为FFT并行程序运行时间，h为效率提升的幅度，时间单位是10-7s。

表4-1奇偶两部分FFT算法性能对比数据表

根据表4-1做图4-1。在双核机器上运行此算法，趋势相同，取时间平均值。为了使算

法性能的对比效果更明显，以指数的形式表示数据，其中横坐标为数据量，N代表2N个，纵坐标为时间T，T代表2T*10-7秒，做图4-2。

由图4-2可知，当N小于212(4096)时，并行算法的运算速度没有串行算法速度快，当数据总数N达到212时，并行算法与串行算法的效率相当。当N大于212时，并行算法的优势就表现出来。由图4-1可知，随着数据量的增加，并行算法的优势越来越明显，当达到4194304时并行运算较串行运算效率提升39％左右。随着数据量的进一步增加，并行算法的效率提升的幅度趋于平稳，基本上在40％左右。出现这种情况的原因是线程的创建、关闭、线程间的切换、线程间的同步需要一定的时间，同时数据的合并也需要一部分时间。因此当很小时，算法的运行时间比线程维护时间小的多，使得并行算法的运行时间比串行算法的运行时间长几十甚至几百倍。当N到达一定的数量级时，算法的运行时间远远大于线程维护时间，且两个线程可以并行运行，此时并行算法的优越性才体现出来。但是由于数据集合并还需要一部分时间，所以随着数据量的增加并行算法的效率可以提升40％左右。

图4-1 FFT算法性能对比图

图4-2 FFT算法性能对比图

当处理大批数据时，并行快速傅里叶变换算法在多核处理器上的运行效率较串行快速傅里叶变换算法有明显的提高，从而证明了并行傅立叶变换应用的可行性。由于多核的普及，采用并行编程思想开发基于多核的应用程序变得非常重要。多核技术的发展必将使软件研发领域发生基础性变化，特别对那些面向一般应用、运行在PC和低端服务器上的应用软件更有非同一般的意义。

4.1.2 蝶形变换内部多线程效率分析

为了比较将蝶形变换内部分解为多个线程的并行傅立叶变换算法的效率，将和串行傅立叶变换的算法在采用不同个数的数据与其进行比较。结果如表4-2所示，其中N为实验数据个数，T1为串行程序运行速度，T2为并行程序运行速度。

表4-2蝶形变换内部多线程FFT性能对比数据表

由表4-2可知，当N很小时，并行算法的运行速度比串行要慢，当N逐渐增大时，并行算法的效率和串行的相比是反而更差。出现这种现象主要有以下两种原因：

(1)在循环体内频繁的创建线程和关闭线程，创建线程虽然没有创建进程所需资源多，

但是频繁的创建线程，对系统资源的消耗是很严重的，因此会出现创建多线程后速度不增反降的现象。

(2)线程同步：因为下次循环会用到上次循环的结果，所以当一次循环没有结束时，虽然建立了两个线程；即便一个线程运行结束，它也不能立即进入下一层循环，它必须等到另外一个线程运行结束，才可以一起进入下一层循环，这就是两个线程同步的问题，两个线程相互等待会使硬件资源处于空闲状态，不被充分利用，这和多核多线程设计的思想相违背。基于以上两点的分析，若不能正确运用多线程编程技术，不但达不到预期的效果，反而使程序运行的效率降低，只有正确的采用并行编程思想，才会使应用程序性能大幅度提升。

5 原程序

#include \"stdafx.h\"

/*int main(int argc, char* argv[]) { }*/

/*****************main programe********************/ #include #include #include #include #define number 65537 struct compx

printf(\"Hello World!\\n\"); return 0;

{ };

double result[number]; struct compx s[number]; //int Num=256; const float pp=3.14159;

struct compx EE(struct compx b1,struct compx b2) {

struct compx b3;

b3.real=b1.real*b2.real-b1.imag*b2.imag; b3.imag=b1.real*b2.imag+b1.imag*b2.real; return(b3); }

void FFT(struct compx *xin,long N) {

double real,imag;

long f,m,nv2,nm1,i,k,j=1,l; /*int f,m,nv2,nm1,i,k,j=N/2,l;*/ struct compx v,w,t; nv2=N/2; f=N;

for(m=1;(f=f/2)!=1;m++){;} nm1=N-1;

/*变址运算*/ for(i=1;i <=nm1;i++) {

if(i while(k }

{

long le,lei,ip; float pi;

for(l=1;l <=m;l++) { le=pow(2,l);// 这里用的是L而不是1 lei =le/2; pi=3.14159; v.real=1; v.imag=0; w.real=cos(pi/lei); w.imag=-sin(pi/lei); for(j=1;j <=lei;j++)

{

/*double p=pow(2,m-l)*j; double ps=2*pi/N*p; w.real=cos(ps); w.imag=-sin(ps);*/ for(i=j;i <=N;i=i+le) { /* w.real=cos(ps); w.imag=-sin(ps);*/ ip=i+lei;

t=EE(xin[ip],v);

xin[ip].real=xin[i].real-t.real;

!!!!

}

}

xin[ip].imag=xin[i].imag-t.imag; xin[i].real=xin[i].real+t.real; xin[i].imag=xin[i].imag+t.imag; }

v=EE(v,w); }

} return;

///////////////////////////////////////////////////////// main() {

clock_t c_start, c_end,time;

unsigned int i=1;

int j;

c_start = clock(); for(j=0;j<100;j++) {

for(i=1;i s[i].real=sin(pp*i/32); s[i].imag=0; }

FFT(s,number-1);

for(i=1;i void FFT(double *dataR, double *dataI) {

int i,j,k; int r,step;

double sin_tap,cos_tap;

double real1,image1,real2,image2; double ansR[NUM], ansI[NUM];

int temp = (int)(log(NUM*1.0)/log(2.0)+0.5); time = c_end - c_start; printf(\"time = %d\\n\}

c_end = clock();

{

result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2)); }

/***************************following code FFT*********************************/

for(i=1;i<=temp;i++) {

step=NUM >> i;

for(k=0; k}

}

}

dataR[r+k]=real1+real2; dataI[r+k]=image1+image2;

dataR[r+k+step]=(real1-real2)*cos_tap-(image1-image2)*sin_tap; dataI[r+k+step]=(real1-real2)*sin_tap+(image1-image2)*cos_tap;

real1=dataR[r+k]; image1=dataI[r+k]; real2=dataR[r+k+step]; image2=dataI[r+k+step];

for(r=0;r<=step-1;r++) {

cos_tap=cos(pi*r/step); sin_tap=sin(pi*r/step);

//开始整序 for(i = 0; i< NUM; i++)

{

ansR[i]=0, ansI[i]=0; int ttemp = 0, k = i; for(j = 0; j <= temp-1; j++) { }

ansR[ttemp] = dataR[i];

ttemp*=2;//左移（本句等价于ttemp <<=1;） ttemp += k%2; k = k >> 1;//右移

}

ansI[ttemp] = dataI[i];

for(i = 0; i < NUM; i++) {

dataR[i] = ansR[i]; dataI[i] = ansI[i];

}/*整序完成*/

}/*end FFT*/

/**********************************end FFT***************************************/

int main(){

double dataR[NUM], dataI[NUM]; int i;

for(i = 0; i < NUM; i++) dataR[i] = dataI[i] = 0; dataR[0] = 1; dataR[1] = 1; dataR[2] = 1; dataR[3] = 1;

dataI[0] = 1; dataI[1] = 1; dataI[2] = 1; dataI[3] = 1;

double s = clock(); FFT(dataR, dataI);

}

double e = clock();

for(i = 0; i < NUM; i++) { }

dataR[i]=dataR[i]/NUM; dataI[i]=dataI[i]/NUM;

cout <printf(\"%lf\cout<<\"ms\"<