并行FFT实现文档格式.docx
《并行FFT实现文档格式.docx》由会员分享,可在线阅读,更多相关《并行FFT实现文档格式.docx(13页珍藏版)》请在冰豆网上搜索。
由式(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,K<
N-1,n,k为整数,N为偶数),x0(i)=x(2i),x1(i)=x(2i+1)(0<
=i<
=N/2-1,i为整数)。
若X0(k)=DFT[x0(i)],若X1(k)=DFT[x1(i)],0<
=k<
=N/2-l,则有:
x(k)=x0(k)+X1(k)*wkN,X(k+N/2)=X0(k)-x1(k)*wkN(0<
=N/2-1,k为整数)。
2.1并行可行性分析
由定理可知,基二时分FFT算法在开始时就可以将要运算的数据分成对等的两部分,每部分都可以独立运算,这符合并行程序上的数据分解;
且FFT运算内部不用建堆,运算后的数据可以替换先前的数据,这同样也满足了多线程的特点。
若是必须要在线程内部建堆,那么当堆到一定程度时,就会出现内存读写错误。
而此算法可以很好的避开线程内部建堆。
即本算法适用于并行计算。
2.2难点分析
本算法的难点有两处,一是数据分解的问题,一般情况下可以将数据平均分为两部分,但是基二时分FFT算法不行,因为此算法的后一步会用到上一步的结果,如果直接分成两部分得出的结果一定错误。
根据算法本身的特点,再结合并行程序的特性,将数据以序列的形式分为奇偶两部分,将这两部分分别以基二时分FFT算法进行运算,得出的结果再根据定理进行汇总,以得出最终正确的结果。
第二个难点在于基二时分FFT算法本身,算法中的蝶形变换本身就是一个三重循环体,每次内部循环的结果都会在下次循环被用到。
图1-11为八点蝶形变换情况:
图1-1八点蝶形变换图
2.3FFT算法的实现
2.3.1计算旋转因子
旋转因子的计算应分为两部分,一是先计算以数据量为整体数据量的一半的旋转因子,二是计算整体数据量的旋转因子的头半部分。
依照数据的顺序将数据分解成为奇偶两部分,且这两部分的数据量相等,因为FFT的数据量要求是2的n次方(n为整数)。
2.3.2创建线程
创建两个线程,每个线程都以FFT的方法运行上面的一个数据块,因这两个数据块的大小相同,所以这两个线程运行的时间也基本相同,这样可以减少线程间等待的时间,创建线程函数如下:
HANDLEhThread[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(intpk=0;
pk<
n分解为二进制的位数;
pk++)
得出后两层所需数据:
nl=2*n1,n2=n1/2,d=pk
第二层循环为:
for(j=0;
j<
n2;
j++)
第三层循环将其分为两个线程,首先得出第三层循环的两个界限为:
实现
过程与以奇偶形式划分数据集的FFT方式同。
4实验结果分析
4.1FFT算法效率分析
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可知,随着数据量的增加,并行算法的优势越来越明显,当达到时并行运算较串行运算效率提升39%左右。
随着数据量的进一步增加,并行算法的效率提升的幅度趋于平稳,基本上在40%左右。
出现这种情况的原因是线程的创建、关闭、线程间的切换、线程间的同步需要一定的时间,同时数据的合并也需要一部分时间。
因此当很小时,算法的运行时间比线程维护时间小的多,使得并行算法的运行时间比串行算法的运行时间长几十甚至几百倍。
当N到达一定的数量级时,算法的运行时间远远大于线程维护时间,且两个线程可以并行运行,此时并行算法的优越性才体现出来。
但是由于数据集合并还需要一部分时间,所以随着数据量的增加并行算法的效率可以提升40%左右。
图4-1FFT算法性能对比图
图4-2FFT算法性能对比图
当处理大批数据时,并行快速傅里叶变换算法在多核处理器上的运行效率较串行快速傅里叶变换算法有明显的提高,从而证明了并行傅立叶变换应用的可行性。
由于多核的普及,采用并行编程思想开发基于多核的应用程序变得非常重要。
多核技术的发展必将使软件研发领域发生基础性变化,特别对那些面向一般应用、运行在PC和低端服务器上的应用软件更有非同一般的意义。
4.1.2蝶形变换内部多线程效率分析
为了比较将蝶形变换内部分解为多个线程的并行傅立叶变换算法的效率,将和串行傅立叶变换的算法在采用不同个数的数据与其进行比较。
结果如表4-2所示,其中N为实验数据个数,T1为串行程序运行速度,T2为并行程序运行速度。
表4-2蝶形变换内部多线程FFT性能对比数据表
由表4-2可知,当N很小时,并行算法的运行速度比串行要慢,当N逐渐增大时,并行算法的效率和串行的相比是反而更差。
出现这种现象主要有以下两种原因:
(1)在循环体内频繁的创建线程和关闭线程,创建线程虽然没有创建进程所需资源多,但是频繁的创建线程,对系统资源的消耗是很严重的,因此会出现创建多线程后速度不增反降的现象。
(2)线程同步:
因为下次循环会用到上次循环的结果,所以当一次循环没有结束时,虽然建立了两个线程;
即便一个线程运行结束,它也不能立即进入下一层循环,它必须等到另外一个线程运行结束,才可以一起进入下一层循环,这就是两个线程同步的问题,两个线程相互等待会使硬件资源处于空闲状态,不被充分利用,这和多核多线程设计的思想相违背。
基于以上两点的分析,若不能正确运用多线程编程技术,不但达不到预期的效果,反而使程序运行的效率降低,只有正确的采用并行编程思想,才会使应用程序性能大幅度提升。
5原程序
#include"
stdafx.h"
/*intmain(intargc,char*argv[])
{
printf("
HelloWorld!
\n"
);
return0;
}*/
/*****************mainprograme********************/
#include<
time.h>
math.h>
stdio.h>
stdlib.h>
#definenumber65537
structcompx
doublereal,imag;
};
doubleresult[number];
structcompxs[number];
//intNum=256;
constfloatpp=3.14159;
structcompxEE(structcompxb1,structcompxb2)
structcompxb3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}
voidFFT(structcompx*xin,longN)
longf,m,nv2,nm1,i,k,j=1,l;
/*intf,m,nv2,nm1,i,k,j=N/2,l;
*/
structcompxv,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<
j){t=xin[j];
xin[j]=xin[i];
xin[i]=t;
k=nv2;
while(k<
j){j=j-k;
k=k/2;
j=j+k;
}
longle,lei,ip;
floatpi;
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/