FFT算法（精选7篇）
FFT算法 篇1
1 FFT处理器的设计与实现
1.1 FFT处理器的总体设计
FFT数据处理器主要由三个模块组成, 各个模块的作用如下:
第一, in_ram的作用在于将按照一定顺序输入的每个点保存起来, 包括实部数据和虚部数据。在所有的16个点的数据一一输入之后, 以并行的方式传递给FFT运算单元。
第二, fft_16_8是FFT运算单元的最关键部分, 其作用在于运用复数FFT并行迭代处理16个点的8位数据;其运算过程是, 首先通过设计, 完成一级FFT运算的gout模块, 在此基础上经过四次迭代调用, 完成16个点的FFT运算。
第三, out_ram的作用在于, 将FFT运算的16个点的运算结果进行保存。在此基础上, 按照输入顺序将每个点的实部和虚部的数据以串行的方式输出。
1.2 各运算模块的设计和实现
1.2.1 加法器的设计
FFT运算的加法器是指一个有时序的, 8位带进位的二进制补码输入加法行波。其运算过程是:通过数据输入, 在2个时钟的运算周期内实现一次二进制补码加法运算。
实现1位全加器的逻辑表达式如式 (1-1) :
其中Ai与Bi表示两个不同的二进制数字;Ci表示一个进位输入数字;Si表示两个不同的二进制的数字输出和。
1.2.2 乘法器的设计
FFT运算的乘法器是指一个有时序的, 8位带进位的二进制补码运用带符号位的阵列的输入行波乘法算法。其具体的逻辑结构为 (n+1) 位× (n+1) 位带求补级阵列。
1.2.3 蝶形运算单元设计
FFT运算的蝶形运算是指首先分别输入2个点的实部数据和虚部数据 (表示为xrl、xil和xr2、xi2) , 这两个点的实部和虚部数据均用8位补码数表示;其次分别输入旋转因子的实部数据和虚部数据 (表示为w nr和w ni) , 这两个点的实部和虚部数据用8位带符号位的原码表示。其各个时钟周期的运算公式如下所示:
1) 四次乘法运算
从上述算法中可以看出, 第二个点的实部数据和虚部数据, 以及旋转因子的实部数据和虚部数据经过运算后, 相比之前扩大了26倍, 并用8位补码数呈现。四个乘法器的输出:Mullout、mul2out、m ul3out、m ul4out也同样扩大了26倍, 并用16位补码数呈现。
2) 一次中间减法运算和一次中间加法运算
上述公式中, 中间加法器与减法器的输出结果分别用midsubout、m idaddout表示, 以8位数的补码数呈现。并在运算过程中对数据进行了相关处理, 使其结果固定缩小2倍。
3) 最后输出前的两次加法运算和两次减法运算
上述公式中, 蝶形运算的两个序列的实部数据和虚部数据分别用xkrl、xkil、xkr2、xki2表示, 用9位补码数表示, 其目的在于防止数据溢出。并在运算过程中, 运用相关处理, 使输入数据固定缩小2倍。由此, 经过蝶形运算后的结果也将固定缩小2被, 最后输出时截取码值的低7位和符号位。
根据补码加法和补码减法的运算公式, 如式 (1-12) 和式 (1-13) 所示:
[x]补+[y]补=[x+y]补式 (1-12)
[x-y]=[x]补-[y]补=[x]补+[-y]补式 (1-13)
从上述公式中可以看出, 负数的加法需要运用补码化的加法运算;同理, 减法可以将其先转化为式 (1-13) , 再进行加法运算。
1.2.4 一级FFT单元设计
由于本次FFT运算是16个点的基-2FFT, 通过上述加法器运算、乘法器运算、蝶形运算等多算法的讨论得出, 完成一级FFT运算需要8个蝶形运算单元, 而且这8个蝶形运算单元是并行运行的。
1.2.5 FFT运算单元设计
综上多运算的讨论, 在这一步骤中, 首先将16个点的实部数据和虚部数据并行输入, 然后通过四次迭代运算, 调用gout模块, 根据基-2FFT, 按照FFT的时序序列算法特点, 在四次迭代运算后并行输出FFT的实部数据和虚部数据结果。
2 综合及结果分析
通过上述多算法和多单元模块的设计与运算, 运用Xilinx的Spartan-3E器件, 遵循FPGA开发的方法和步骤, 首先运用Ve rilog HDL进行代码编写, 并完成语法检查, 即Che ck Syntax。在此基础上, 调用“P1an Ahead”工具, 自动分配I/O引脚。最后运用Xilinx ISE自带的综合工具Synthe s ize-XST进行编译、分析、综合, 执行Ge ne rate Pos t-Synthe s is Sim ulation Mode l, 生成后仿模型, 同时产生综合报告文件 (.syr) 。
3 总结
本文以基于FPGA的FFT算法的设计与实现作为选题, 掌握了FPGA设计的基本流程。同时, 用Matlab工具对本设计的算法流程进行了验证, 并且借助Model Sim和Matlab工具对设计结果进行了联合仿真。
参考文献
[1]王旭东, 靳雁霞.MATALB及其在FPGA中的应用[M].国防工业出版社, 2008.
[2]高西全, 丁玉美.数字信号处理[M].西安电子科技大学出版社, 2008.
[3]李全利, 李静.基于FPGA的FFT算法研究[J].自动化技术与应用, 2007.
FFT算法的FPGA实现 篇2
离散傅里叶变换(Discrete-time Fourier Transform,DFT)是信号分析与处理中的一种重要变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。1965年由J.W.Cooley和J.W.Tukey提出了FTT算法，它把计算N点DFT的乘法运算量从N2次下降到N/2 log2N次，使DFT的运算效率提高了1 2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了良好的条件，大大推动了数字信号处理技术的发展。FFT算法的出现也因此被看作是数字信号处理发展史上的里程碑。目前，FFT在图像处理、石油勘探、地震和雷达等领域中都得到了广泛的应用。对于FFT的工程实现，目前比较通用的方法有三种：数字信号处理器(DSP)、专用FFT处理芯片和FPGA。其中DSP适用于流程复杂的算法，利用DSP进行FFT运算将占用DSP大量的运算时间，降低整个系统的数据吞吐率，也无法发挥DSP的灵活性;若用专用FFT处理芯片，尽管速度尚可，但外围电路复杂，可扩展性差，并且价格昂贵;FPGA具有可重构的特点，适合于算法结构固定、运算量大的前端数字信号处理。采用FPGA技术可以提高元器件的优质利用性，降低设计成本，并能够并行处理数据，容易实现流水线结构，而且升级简便，提高了设计的灵活性，所有这些都非常适合实现FFT算法。因此，本文主要研究基于FPGA芯片的FFT，把FFT对实时性的要求和FPGA芯片设计的灵活性结合起来，实现并行算法与硬件结构的优化配置，提高FFT处理速度，满足现代信号处理的高速、高可靠性要求。
1 算法原理
DFT是数字信号处理中常用的一种算法，用来对信号频谱进行分析，其基本公式是：
式中x(n)表示时域信号，X(k)表示频域信号，为旋转因子。由式(1)可以看出，对于N点的DFT需做N2次复数乘法运算和N(N-1)次复数加法运算，乘法和加法的运算量均与N2成正比，当N较大时，运算量相当可观。FFT算法就是不断地把长序列的DFT分解成几个短序列的DFT，并利用的周期性和对称性来减少DFT的运算次数。
设序列x(n)的长度为N(N=2M,M为任意正整数)，按n的奇偶把x(n)分解成两个N/2点的子序列：
则x(n)的N点DFT为
由于
所以
其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点的DFT。由于X1(k)和X2(k)均以N/2为周期,且,所以X(k)又可表示为
这样就将N点的DFT分解为两个N/2点的DFT以及上面两式的运算,这种运算可用如图1所示的流图符号表示,根据其形状称之为蝶形运算符号。其中
按照这种分解运算的思想，将X1(k)和X2(k)继续向下分解，直到最后变为一组2点DFT运算的组合。以8点的FFT运算为例，其运算过程如图2所示[1]。
2 设计思路
2.1 数的定标
在蝶形运算单元中，，即在蝶形单元运算中存在复数，且它的实部和虚部均为小数。而在FPGA的设计中，由于输入输出均为固定数目的引脚，操作数一般采用整型数来表示，所以必须首先通过数的定标，用整型数来合理地表示小数。数的定标常用的有Q表示法。在小数和整数的互相表示中，存在以下关系[2]：小数(x)转换为定点的整数(xq)：;定点的整数(xq)转换为小数(x)：
2.2整体设计
整体设计可以分为输入、输出两个部分，输入部分主要由数据接收处理模块组成，输出部分主要由运算处理部分组成。两个模块之间的数据传递通过地址生成及数据查找模块完成。实现FFT算法各个模块的实现框图如图3所示。
2.3 各功能模块的实现
2.3.1 数据接收处理模块
如图2所示，进行第0级运算时并不是按地址顺序进行取数，而是需要数据接收处理单元把存放在接收RAM中的每个采样点按地址读出，按码位倒置后再存入运算RAM中准备进行运算。
其中运算RAM的作用是为FFT提供数据缓存，并且交替作为数据读出和运算结果写入单元，直到第9级蝶型运算完成。运算RAM采用双口RAM模式，允许同时将蝶型运算所需要的两个入口数据读出，也可以同时将两个运算结果写入，这样就比单口RAM节省1/2的读写周期。
2.3.2 蝶型运算模块
蝶型运算模块的功能是完成对输入数据A、B的蝶型运算，采用标准的加法器和乘法器来完成，以保证运算的精确性。由于输入数据A、B和旋转因子均为32位，为了保证没有溢出错误，乘法器输出应为64位，加法器输出应为65位，每级蝶型运算结果动态范围最大需要扩展33位;因此需要对蝶型运算结果进行检测，判断当前结果动态范围扩展位数，计算出当前级的最大扩展位数。每级FFT运算结果的动态范围最多需要扩展5位，则中间运算结果的位宽应为37位。在送回到运算RAM之前，按照动态范围将这37位数据扩展标准的块浮点格式，以完成定阶。最后，截取其中的32位送回运算RAM为下一级运算做准备。
2.3.3 地址生成及数据查找模块
在FFT运算中，中间运算数据和旋转因子的寻址是实现设计的关键，在设计中利用额外的两块片上RAM存储FFT中间运算结果的实部和虚部，利用两块片上ROM(宽度为32位，深度为10位)存储旋转因子的实部与虚部，这种方法消耗大量的片上资源。一般情况下，N点FFT只用log2N位地址码、N/2旋转因子仅需要位地址码来建立查找表，而在设计中选择的RAM数据宽度为32位，深度为10位，ROM数据宽度为32位，深度为10位，比地址查找表的方法使用了更多的片上资源，但是这样做大大加快了寻址速度。
3 仿真分析
本文采用Alter公司的CycloneⅡ系列FPGA芯片EP2C35F672C8，用VHDL语言编程，并在QuartusⅡ6.0平台上进行逻辑综合与时序分析。在设计中选择初相为0，频率为f1(100Hz)/f2(200Hz)的正弦信号作为输入信号，输入数据是对此输入信号的采样，采样数为1024点，数据格式为32bit补码，存储在两块片上单口RAM中。经FFT运算的中间结果和最后结果输出分为实部、虚部两路，数据格式为每路32位补码，存储在另外两块片上双口RAM中。旋转因子的实部、虚部和输入的数据位相同，都为32位，存储在两块ROM上。图4为QuartusⅡ6.0的仿真结果，
仿真结果显示共占用了FPGA中34%的逻辑单元，50%的片上RAM和2 5%的全局时钟，时钟最快可达102.353MHz。完成1024点FFT，需要时钟约49μs。
将FFT的计算结果依次读出后，利用Matlab软件描点，得到的波形图如图5所示。由图5可以看出，如果排除随机噪声的影响，计算结果和Matlab软件计算得到的结果(如图6)基本一致。
4 结束语
通过比较可以看出仿真结果与理论值吻合得很好，所采用的方案能达到设计的要求，并且控制实现简单，具有较高的运算速度。随着可编程器件规模、速度的提高，利用FPGA实现FFT运算的优势越来越明显，并且这种方案实现的FFT很容易扩展，具有很大的可行性和实用性。
摘要：在信号处理中,FFT占有很重要的位置,其运算时间影响整个系统的性能。传统的实现方法速度很慢,难以满足信号处理的实时性要求。针对这个问题,本文研究了基于FPGA芯片的FFT算法,把FFT算法对实时性的要求和FPGA芯片设计的灵活性结合起来,采用Alter公司的CycloneⅡ系列FPGA芯片EP2C35F672C8,用VHDL语言编程,最后分别使用Quartus Ⅱ和Matlab软件开发工具验证实现。
关键词：数字信号处理,快速傅立叶变换,现场可编程门阵列
参考文献
[1]丁玉梅,高西全.数字信号处理(第二版)[M].西安:西安电子科技大学出版社,2001.
[2]侯伯亨,顾新.VHDL硬件描述语言与数字逻辑电路设计[M].西安:西安电子科技大学出版社,2000.
[3]李静梅,郑超峰,金玉苹.基于FPGA的FFT算法实现[J].应用科技,2009,36(2):38-41.
FFT算法 篇3
随着电力电子技术的飞速发展及电力电子设备的广泛应用,电能质量污染问题日益严重,其中谐波问题尤为突出。谐波对电力系统安全、稳定运行构成了严重威胁,对电气环境造成了严重影响,因此谐波补偿至关重要,而谐波检测是设计谐波补偿装置的前提。电网谐波检测的方法主要有使用模拟低通滤波器并进行离散傅里叶变换的方法及小波分解方法等。低通滤波器结构简单,但误差较大,已逐渐被淘汰。小波分解方法精度较高,抗干扰能力强,但其运算速度慢,不易实现。因此,本文采用离散傅里叶变换进行谐波检测,并采用先进的定点DSP芯片TMS320F2812实现FFT算法。
1 DSP芯片定点运算
TMS320F2812具有32位硬件乘法器和累加器,处理能力可达150 MIPS,但它是一款定点处理芯片,用定点数进行数值运算,其操作数一般采用整型数来表示,不能直接处理小数。因此本文用Q格式[3]表示操作数,并进行小数运算。在进行FFT运算时,把输入数据归一化为Q15格式来表示,既可使运算结果达到一定精度,又可以提高芯片TMS320F2812处理浮点数的速度。Q15格式数据范围为[-65 536,65 535]。计算完成后,将Q15格式的所有结果转换为浮点数即得到正确结果。
浮点数xf转换为定点数xq:
定点数xq转换为浮点数xf:
2 电网谐波检测方法实现
2.1 FFT谐波检测一般步骤
用FFT进行谐波检测的步骤:(1)采样电网电压信号并将其变换为离散序列。(2)建立数据窗,忽略数据窗前后信号波形。(3)进行FFT运算得到结果。在进行以上步骤时,必须满足以下要求:首先要满足采样定理,以免引起频谱混叠;其次采样频率必须与信号频率同步,整周期采样。
对电网信号x(t)以频率fs采样N个数据,得到一组序列x(n)(n=0,1,…N-1)。序列x(n)实质上可理解为对连续信号x(t)离散采样后,用长度为N的矩形窗截断的结果。由时域卷积定理可知,时域相乘在频域对应卷积,所以对原始信号进行FFT后得到的是原始信号的频谱与矩形窗频谱的卷积,这将造成长范围频谱泄漏[5]。由于电网信号频率存在小范围内的漂移,所以在采样时不可能做到严格同步采样,对此时的采样结果进行频谱分析存在短范围频谱泄漏。为减小频谱泄漏对检测结果造成的影响,采用加窗插值FFT算法[6]进行谐波分析。
2.2 数据采集模块设计
在实际测量中,必须对输入电压及电流进行必要的处理才能作为采样单元的输入。本文选择电压互感器SPT204A来获得精度高、线性度好的输出交流电压。SPT204A是一款微安级精密电压互感器,其输入、输出额定电流均为2mA。它具有精度高、轻便、能直接焊于印刷线路板上等特点。电压互感及调理电路如图1所示,其中UIN为电网电压信号。SPT204A输出信号经过两个反向并联二极管D1、D2输入运算放大器[4]。
由于电压互感及调理电路输出信号是峰值约为3.3V的双极性交流信号,而所采用的AD转换器要求输入0~3.3V的单极性信号,因此还需设计电压偏置电路来进行信号转换。电压偏置电路如图2所示,其主要器件是运算放大器。
2.3 加窗插值FFT算法原理
加窗插值FFT算法常用的窗函数有汉宁(Hanning)窗、Blackman-Harris窗等余弦窗。本文采用汉宁窗:
对采样序列x(n)进行加窗处理后得到x′(n):
对x′(n)进行FFT得到X′(k):
汉宁窗对应的FFT插值公式[7]为
式中:Ak和Φk分别对应各次谐波的振幅和相位;
2.4 谐波检测流程
基于FFT算法的谐波检测流程如图3所示。首先初始化DSP及各模块;然后利用内部定时器产生中断,以固定采样频率采集电压信号,经AD转换得到离散序列x(n);将x(n)转换为Q15格式进行加窗、FFT及插值处理,得到Q15格式的各次谐波幅值和相位;最后利用式(2)将结果转换为浮点数。
3 仿真分析
对普通电网基波频率为50Hz的周期信号进行谐波检测,固定采样频率为12.8kHz,采集点数N=1 024。
图4为采集到的电压原始波形。利用Matlab对采集到的数据进行算法仿真,前20次谐波幅值结果如图5所示,其中第1个点为直流分量,第2个点为基波幅值,以此类推。
采用DSP计算各次谐波幅值和相位,结果如表1所示。从表1可看出,DSP计算得到的基波分量为309V,由于该信号为正弦信号,可计算出其有效值为218V,符合实际情况。对比DSP计算结果与Matlab仿真结果可知,利用DSP实现的加窗插值FFT算法能准确地计算出各次谐波,误差较小。
4 结语
介绍了基于FFT算法的电网谐波检测方法的具体实现及仿真分析。该方法在高速定点DSP芯片上实现FFT算法,将数据归一化为Q15格式进行计算,并采用加窗插值FFT算法进行谐波分析,既保证了检测速度,又可使结果满足一定精度。
摘要：针对传统电网谐波检测方法误差较大、运算速度慢的问题,提出了基于FFT算法的电网谐波检测方法;介绍了DSP芯片的定点运算、基于FFT算法的电网谐波检测方法的具体实现,并进行了仿真分析。该方法在高速定点DSP芯片TMS320F2812上实现FFT算法,将数据归一化为Q15格式进行计算;为减小频谱泄漏对检测结果造成的影响,采用加窗插值FFT算法进行谐波分析。该方法实现了对电网电压各次谐波的检测,保证了检测速度和精度。
关键词：电网,谐波检测,FFT,加窗插值,DSP定点运算
参考文献
[1]程佩清.数字信号处理教程[M].2版.北京:清华大学出版社,2001.
[2]王艳芬,王刚,张晓光,等.数字信号处理原理及实现[M].北京:清华大学出版社,2008.
[3]王潞钢,陈林康,曾岳南.DSP C2000程序员高手进阶[M].北京:机械工业出版社,2005.
[4]远坂俊昭.测量电子电路设计——模拟篇[M].彭军,译.北京:科学出版社,2006:17-23.
[5]庞浩,李东霞,俎云霄,等.应用FFT进行电力系统谐波分析的改进算法[J].中国电机工程学报,2003(6):50-54.
[6]谢明,丁康,莫克斌.频谱校正时谱线干涉的影响及判定方法[J].振动工程学报,1998(1):55-60.
FFT算法 篇4
1 实序列FFT实现原理
FFT是离散傅里叶变换(DFT)的快速算法。DFT是将时域序列变换为长度相同的频域序列,从理论上讲,输入序列为复序列。对于N点离散有限长时间序列x(n),离散傅里叶变换为[4,5]
但在实际工程实践中,输入序列x(n)通常是实序列。如果直接按照FFT运算流图计算,就是把x(n)看成一个虚部为0的复序列进行计算,这必将增加存储量和运算时间。根据数字信号处理理论[6,7],采用一种更加优良的方法实现32点实序列的FFT运算。把32点实序列DFT计算转换为一个16点复序列DFT计算,然后将16点复序列的输出进行适当的运算组合就可得到原始序列的DFT。x(n)为N=32点实序列,取x(n)的偶数点f(n)和奇数点g(n)分别作为序列h(n)的实部和虚部构造新序列,有下式
则h(n)就是一个16点复数序列,用FFT计算h(n)的16点离散傅里叶变换H(k),就可通过H(k)推导得出F(k)和G(k)分别为
求出F(k)和G(k)后就可利用下面的公式计算出x(n)的离散傅里叶变换X(k)为
这种算法比相同长度的复序列FFT算法减少约一半的运算量,运算效率提高了近一倍。
2 实序列FFT硬件实现
2.1 系统总体设计
目前FFT硬件实现结构主要分为三类:递归结构、流水线结构以及并行迭代结构。递归结构又叫作内存共享结构,这种结构的主要优点是只有一个运算单元,占用硬件资源少,缺点是运算时间较长。流水线结构每一级均采用一个运算单元,前一级运算结果直接用于下一级,因此速度上有所提高。并行迭代结构每级都采用N/2个蝶形单元进行并行运算,每一列中的运算并行进行,一列迭代至另一列是顺序进行的。这种方法处理速度快,但设备量较大需要采用大规模集成电路实现。考虑到本设计点数较少,对处理速度要求较高,因此采用并行迭代结构实现。设计主要包括以下几个模块:实/复序列转换模块、输入数据串并转换模块、基2-DIT算法实现模块、输出数据并串转换模块以及还原实序列DFT模块。实/复序列转换模块主要完成的功能为:将32点实序列转换为16点复序列,然后再进行FFT运算。还原实序列DFT模块是将16点复序列h(n)的FFT运算结果根据式(3)、式(4)求出32点实序列x(n)的FFT运算结果。基2-DIT算法实现模块是整个系统实现的核心模块,因此主要介绍基2-DIT算法实现模块。因为输入的32点实序列通过实/复转换模块将变成16点复序列,所以基2-DIT实现模块主要是针对16点复序列设计完成的。
2.2 基2-DIT算法硬件实现
基2-DIT算法实现的核心是蝶型运算单元的实现,蝶型运算单元的结构框图如图1所示。蝶型运算单元的设计主要包括以下几个模块:复数乘法模块、延时D触发模块、截取模块、缩小模块以及实数加减模块,文中主要介绍复数乘法模块、截取模块、缩小模块。
复数乘法模快采用IP核ALTMULT-COMPLEX实现,主要完成旋转因子与一个输入数据的复数乘法。考虑到16点FFT设计只需8个旋转因子,因此可以先通过Matlab计算出所需的旋转因子并存储在ROM里,按照每级运算的需求直接读取降低运算量。由于旋转因子取值范围在[-1,1]之间,为了方便定点运算,将旋转因子扩大26并取整数,这样复数乘法器的输出结果也扩大26。由于FPGA内部运算采用定点处理,因此每次复乘之后输出结果的位宽都要扩大一倍。并且,相乘之后的数据还要相加减,为了保证不溢出,输出结果位宽还要扩展一位。所以随着FFT运算点数的增加,输出结果的位宽将以2的指数次幂增加,占用的资源太多以致无法实现。因此,必须增加截取模块。
截取模块实现的功能是将复数乘法器的16位输出crout[15:0]、ciout[15:0]截取8位。通过比较,多种截取方式找到一种误差比较低的截取方式:首先判断crout[6]、ciout[6]是否为1,如果是1,则输出为{crout[15],crout[13:7]}+1、{ciout[15],ciout[13:7]}+1;如果不是1,则输出为{crout[15],crout[13:7]}、{ciout[15],ciout[13:7]}。截取后数据输出结果缩小了27倍,因为做复数乘法之前旋转因子扩大26倍,所以复数乘法器的最终输出结果缩小了2倍。因此也应该将另一个输入数据缩小2倍,才能进行后续的加减操作。因此每级蝶型运算的输出都缩小两倍,这样可以有效的避免输出结果的溢出,FFT最后的输出结果缩小了16倍。因为对于FFT运算,所关心的主要是输出数据的相对幅度,所以这种固定的缩小比例对FFT运算结果影响不大[8]。
缩小模块实现的功能是将不与旋转因子相乘的输入数据缩小2倍。通过移位的方式可以近似实现缩小两倍的功能,这种方式比除法占用更少的逻辑资源。具体的实现方式是:首先判断输入数据datain的符号位datain[7]是0还是1(即判断输入数据的正负),如果是0则输出为{1'b0,datain[7:1]};如果是1则输出为{1'b1,datain[7:1]}。
整个基2-DIT算法的硬件实现分为4级,每级用8个蝶形运算单元实现。前一级运算结束的标志作为下一级运算以及读取旋转因子的启动信号,上一级运算结果直接送到下一级作为输入数据。基2-DIT算法实现模块的输出经过并串转换后送入还原实序列DFT模块进行组合运算,输出结果即为32点实序列FFT运算的最终结果。基2-DIT模块输出结果缩小16倍,还原实序列DIT模块输出结果缩小2倍,所以整个系统的输出缩小了32倍。
3 仿真与分析
系统设计完成后用Signal Tap观察仿真结果。当串行输入32点实序列x(n)=[21 39 58 76 97 119 126117 96 78 53 36 18-10-35-66-38-8 26 38 62 4023-14-30-52-68-76-89-101-118-125]时,经过实/复转换模块输出为16点复序列y(n)=[21+39i 58+76i 97+119i 126+117i 96+78i 53+36i18-10i-35-66i-38-8i26+38i62+40i23-14i-30-52i-68-76i-89-101i-118-125i],将y(n)串并转换后送入基2-DIT实现模块得到16点复序列y(n)的FFT运算结果,仿真波形如图2所示。
最后通过输出数据并串转换模块以及还原实序列DFT模块得到最终的FFT运算结果,Signal Tap仿真结果如图3所示。
为了验证系统设计的准确性,用Matlab计算x(n)的FFT理论值并缩小32倍后的结果如图4所示。
Matlab计算理论结果与Signal Tap输出实际结果幅值比较如图5所示。通过对比可以发现,本设计仿真结果与Matlab计算得到的理论值大致相同,误差较小。FFT运算过程主要的误差有舍入误差、溢出误差、截取误差,由于本设计采用定点运算,所以主要误差为截取误差。如果想要提高精度,可以采用浮点或块浮点运算,但是硬件实现的复杂度也会相对提高。
4 结论
设计采用并行迭代处理方案实现了32点8位实序列FFT,通过Signal Tap与Matlab数据对比发现,该设计可以实现预期的功能。从串行数据全部输入开始到最后开始输出数据总共需要42个时钟周期,本设计在一些数字信号处理领域(比如信道化接收机)具有良好的应用前景。
参考文献
[1]Sanchez M A,Garrido M,Lopez-Vallejo M,et al.Imple-menting FFT based digital channelized receivers on FP-GA platforms[J].Aerospace and Electronic Systems,2008,44(4):1567-1585.
[2]刘美容.FFT算法的DSP实现[J].微电子学与计算机,2015,32(1):76-79,84.
[3]万红星,陈禾,韩月秋.一种高速并行FFT处理器的VL-SI结构设计[J].电子技术应用,2005,31(50):45-48.
[4]高西全,丁玉美.数字信号处理[M].3版.西安:西安电子科技大学出版社,2008:110-123.
[5]吴大正.信号与线性系统[M].4版.北京:高等教育出版社,2009:115-187.
[6]郑南宁.数字信号处理[M].西安:西安交通大学出版社,1991:67-69.
[7]Alan V O,Ronald W S.Discrete-time digital signal pro-cessing[M].Beijing:Tsinghua University Press,2005.
[8]王旭东,刘渝.全并行结构FFT的FPGA实现[J].南京航空航天大学学报,2006,38(2):96-100.
FFT算法 篇5
快速傅里叶变换 (FFT) 作为计算和分析工具, 在众多学科领域 (如信号处理、图像处理、生物信息学、计算物理、应用数学等) 有着广泛的应用。在高速数字信号处理领域, 如雷达信号处理, FFT的处理速度往往是整个系统设计性能的关键所在。
针对高速实时信号处理的要求, 软件实现方法显然满足不了其需要。近年来现场可编程门阵列 (FPGA) 以其高性能、高灵活性、友好的开发环境、在线可编程等特点, 使得基于FPGA的设计可以满足实时数字信号处理的要求, 在市场竞争中具有很大的优势。
在FFT算法中, 数据的宽度通常都是固定的宽度。然而, 在FFT的运算过程中, 特别是乘法运算中, 运算的结果将不可避免地带来误差。因此, 为了保证结果的准确性, 采用定点分析是非常必要的。
1 FFT算法原理
FFT算法的基本思想就是利用权函数的周期性、对称性、特殊性及周期N的可互换性, 将较长序列的DFT运算逐次分解为较短序列的DFT运算。针对N=2的整数次幂, FFT算法有基-2算法、基-4算法、实因子算法和分裂基算法等。这里, 从处理速度和占用资源的角度考虑, 选用基-4按时间抽取FFT算法 (DIT) 。对于N=4γ, 基-4 DIT具有log4N=γ次迭代运算, 每次迭代包含N/4个蝶形单元。蝶形单元的运算表达式为:
其信号流如图1。式中:A, B, C, D和A′, B′, C′, D′均为复数据;W=e-j2π/N。进行1次蝶形运算共需3次复乘和8次复加运算。N=64点的基-4DIT信号流其输入数据序列是按自然顺序排列的, 输出结果需经过整序。64点数据只需进行3次迭代运算, 每次迭代运算含有N/4=16个蝶形单元。
2 FFT算法的硬件实现
2.1 流水线方式FFT算法的实现
为了提高FFT工作频率和节省FPGA资源, 采用3级流水线结构实现64点的FFT运算。流水线处理器的结构如图2所示。
每级均由延时单元、转接器 (SW) 、蝶形运算和旋转因子乘法4个模块组成, 延时节拍由方框中的数字表示。各级转接器和延时单元起到对序列进行码位抽取并将数据拉齐的作用。每级延时在FPGA内部用FIFO实现, 不需要对序列进行寻址即可实现延时功能。数据串行输入, 经过3级流水处理后, 串行输出。
转接器有一定的工作规律。例如, 当第0级变换做完进入转接器SW1前, 先对后三路数据进行一定节拍的延时, 延迟节拍分别为4, 8, 12。为了说明规律, 把输入转接器的四路数据按照前后次序进行分组, 每4个时钟节拍为1组, 共16组, 如图3 (左) 所示。在数据流串行经过转接器SW1时, 第0组中的数据保持不变, 第1组中的数据与第4组中的数据交换;5不变, 2和8交换, 3和12交换, 6和9交换;10不变, 7和13交换, 11和14交换, 15不变。交换完毕后, 前三路数据经过延迟节拍分别为12, 8, 4的FIFO存储器输出, 位置关系如图3所示。
上述转换规律对于SW2也是适用的, 只是转接器前后的延时节拍和分组的大小有所不同。
2.2 存储单元
为了实现算法的流水线设计, 存储器RAM设计为64×16 b的双端口RAM, 即在时钟信号和写控制信号同时为低电平时, 从输入总线写入RAM;在时钟信号和读控制信号同时为高电平时, 从 RAM输出数据。
ROM为17×16 b的ROM, 储存经过量化后的旋转因子, 旋转因子为正弦函数和余弦函数的组合。根据旋转因子的对称性和周期性, 在利用ROM存储旋转因子时, 可以只存储旋转因子的一部分。
2.3 运算结构
Radix-4蝶形运算单元是整个 FFT处理器中的核心部件。在用Radix-4运算器计算时需要并行输入数据, 如果能以并发数据输入的话, 则同步性和控制度较好, 但实际上常要进行串并之间的转换。存储RAM按单节拍输出16 b位宽数据, 选择器不停旋转送入到确定的位置, 每4点全部到位后R-4使能有效;然后4个时钟节拍得到有效结果数据, 再通过选择器旋转送入到对应存储RAM中。
复数运算中, 对应复数的实部和虚部RAM用同一个地址发生器。地址发生器在进行RAM地址发生时采用两套地址, 第一套是计数器按时钟节拍顺序产生的, 用于输入数据的存储;第二套是由数据宽度为16 b的ROM产生的, ROM中存放的数据为下级运算所需倒序的序列地址, 发生地址给RAM, 然后RAM按倒序地址输出下级需要进行运算的数据。
2.4 块浮点结构
数字信号处理系统可分为定点制、浮点制和块浮点制, 它们在实现时对系统资源的要求不同, 工作速度也不同, 有着不同的适用范围。定点制算法简单, 速度快, 但动态范围有限, 需要用合适的溢出控制规则 (如定比例法) 适当压缩输入信号的动态范围。浮点表示法动态范围大, 可避免溢出, 但系统实现复杂, 硬件需求量大, 速度慢。
为了提高精度, 并减少复杂度和存储量, 采用块浮点结构。块浮点算法是以上两种表示法的结合。这种表示方法是, 一组数共用同一个阶码, 这个阶码是这组数中最大数的阶码。块浮点算法无需进行额外的指数运算, 仅对尾数进行运算即可, 其与定点运算一样方便, 但需要在每级运算结束后进行本级运算溢出最大位数判断, 以对数据块进行块指数调整。在调整时仅保留一位符号位, 因而能够充分利用有限位长。这样处理比定点方法扩大了动态范围, 并且提高了精度, 比浮点运算在速度上有了提高。块浮点结构如图4所示。
3 结 语
着重讨论基于FPGA的64点高速FFT算法的实现方法。采用高基数结构和流水线结构, 大大提高了FFT处理器的运行速度。同时块浮点结构的引入, 也大幅减少了浮点操作占用FPGA器件的资源数目, 兼顾了FPGA高精度、低资源、低功耗的特点。从实验结果看, 该方法可以满足高速实时处理数字信号的要求。
参考文献
[1]朱冰莲, 刘学刚.FPGA实现流水线结构的FFT处理器[J].重庆大学学报, 2004, 27 (9) :33-36.
[2]胡广书.数字信号处理——理论、算法与实现[M].2版.北京:清华大学出版社, 2005.
[3]陈丽安, 张培铭.定点DSP块浮点算法及其实现技术[J].福州大学学报:自然科学版, 2004;32 (6) :689-693.
[4]求是科技.VHDL应用开发技术与工程实践[M].北京:人民邮电出版社, 2005.
[5]孔利东.基于FPGA到数据采集与处理技术的研究[D].武汉:武汉理工大学, 2007.
[6]谭征, 张晓林.一种基于FPGA的超高速FFT处理器设计[J].遥测遥控, 2005, 26 (6) :46-49.
[7]任淑艳, 关丛荣.应用VHDL语言的FFT算法实现[J].哈尔滨理工大学学报, 2003, 8 (6) :24-26.
[8]田丰, 邓建国.FFT算法的一种FPGA实现[J].现代电子技术, 2005, 28 (8) :97-100.
[9]孙志坚, 刘学梅.在FPGA中实现高速FFT算法的研究[J].青岛建筑工程学院学报, 2005, 26 (2) :84-86.
FFT算法 篇6
当前GPU的体系结构为高性能计算提供了良好的可编程性[1,2,3,4,5,6]。一个GPU包含多个多处理器(MP),一个多处理器又包含若干个流处理器(SP)。在GPU上具有全局设备存储器可以被所有多处理器访问;每一个多处理器又具有一定量的共享存储器(Shared Memory)供流处理器共享使用和一定量的寄存器供流处理器独占使用。CUDA编程模型为这样的硬件提供了易用的编程方式[7]。它的运行模型包括若干不同的并行层次:线程块网格(Thread Block Gird)、线程块(Thread Block)、Warp块以及线程,它们被调度到GPU上执行。一个线程块只会被调度到一个多处理器上执行,因此,一个线程块中的线程可以共享访问一部分共享存储器。一个线程块网格包括很多线程块。因为线程块的数量可以大大超过多处理器的数量,一个线程块网格中的线程块被动态的调度到多处理器上运行。线程块中的线程被组织成Warp块,一个Warp块以SIMD的方式执行。线程块中的线程以排他的方式被调度到多处理器上执行。
CUDA提供了对众核GPU进行程序设计的底层方式,但是要想得到顶级的性能对程序员来讲确实是一个挑战。程序员必须理解GPU硬件平台的各种制约因素和限制来设计算法,并在影响性能的各种因素之间进行权衡。例如,早期硬件生产商提供的GPU上的数学库:
(1)基本线性代数运算的CUBLAS库;
(2)FFT计算的CUFFT库。
其中的一些基本运算函数(版本为CUBLAS 1.0和CUFFT 1.1)就具有比较差的性能。在本文中,对两个基本的运算函数:计算密集型的矩阵乘法和带宽/通信密集型的1维FFT进行性能优化,并以此说明GPU程序优化技术。
本文优化的代码的性能远远优于早期硬件生产商提供的GPU上的数学库。这说明,对于一个程序员而言,使用CUDA在GPU上得到顶级性能的代码是十分困难的。实际上,程序员不仅要精确了解GPU的硬件体系结构,还要具有丰富的实践经验,来帮助他们达到顶级的存储器带宽、峰值速度等性能指标。在本文中,通过一些基本的测试,得到了相应的编程原则。本文的代码性能为单精度浮点数在NVIDIA Ge Force GTX280上得到的性能。
1 基本优化原则
一系列的GTX200 GPU程序设计优化原则可以在文献[2-4]中找到。在本节中,主要讨论根据所进行的基础实验,得到的一些优化原则。
1.1 以MIMD方式控制线程Warp
在每一个时钟周期,一个多处理器以SIMD的方式运行:即一个多处理器中的每一个微处理器执行相同的指令运算不同的数据[8]。基本的SIMD的单位是Warp,每一个Warp具有相同的线程个数。在GTX200 GPU上,Warp的大小是32。Warp块是按照线程的自然顺序形成的:前32个线程是第1个Warp,而后第2个等。活动的Warp以时间片的方式占用多处理器的计算资源。
在大多数的CUDA程序中,每一个线程块中的线程具有相同的代码结构,只是根据线程ID区分处理不同的数据元素。尽管没有提倡,但实际上不同的Warp块可以根据不同的Warp ID执行不同的代码。这通过在条件语句中测试Warp ID来实现。只要同一个Warp块中的线程执行同样的条件分支,在性能上就不会有明显的损失。这使得代码可以运行在MIMD的模式下并执行不同的指令。
图1演示了这种MIMD机制的实现方式。在这种情况下,使用这种机制可以减少一个线程块的共享存储器的使用量,从而使得更多的线程块可以被加载到一个多处理器上。在FFT计算中,各个线程需要进行全局的矩阵转置:4个Warp块中的各个线程将自己的数据放入到共享存储器中并同步,之后再从共享存储器中将自己需要的数据取回,如图1(a)所示。如果按照图1(b)所示的方式进行,共享存储器的用量可以被减半:前2个Warp块中的线程先将自己的数据放入到共享存储器中并同步,各个线程将自己需要的数据取回,而后后2个Warp块中的线程再将自己的数据放入到共享存储器中并同步,各个线程将自己需要的数据取回。根据实验,具有Warp间条件指令转移的代码会带来性能10%~30%的损失。但是,由于一个多处理器可以具有更多的活动的线程块,在本文FFT的代码实现中,损失的性能可以被由于具有更多的活动线程块掩盖指令流水线延迟得到的性能所掩盖。
1.2 峰值指令吞吐量
通过测试代码可以测试一个GPU硬件所能达到的峰值指令吞吐量,通过比较该值和代码达到的性能之间的差异,可以推测该代码进一步优化的空间。Ge Force GTX280具有240个微处理器,运行在1.295 GHz的频率下,单精度浮点数的理论性能为933 Gflop/s。通过在一个展开的循环中运行寄存器间的MAD指令a=b*a+b和b=a*b+a,可以达到GPU的峰值指令吞吐量。在完全流水线的情况下,可以达到622 Gflop/s的性能。
1.3 流水线延迟
在文献[2-4]中指出,每一个多处理器必须要有192个活动的线程,才能够掩盖流水线的延迟。本文的实验表明,必须至少要有256个活动的线程才能够掩盖流水线的延迟。只有在每一个多处理器有超过256个活动线程,而不论这256个活动线程分成多少个线程块,峰值指令吞吐量620 Gflop/s才可能被达到。
1.4 增加活动线程块以掩盖同步指令
表1示意了不同同步和计算指令比例以及在一个多处理器上具有1个或2个活动线程块的情况下,代码可以达到的指令吞吐量。如果一个多处理器上只有一个活动的线程块,代码中的同步指令会导致一个多处理器成为空闲的状态。如果一个多处理器上有多个活动的线程块,在一个线程块执行同步指令的时候,另一个活动的线程块可以被调度器自动调度执行。这种使用多个活动线程块掩盖线程块同步指令的方法也被运用于FFT的代码优化中。
1.5 设备存储器带宽
设备存储器不具有缓存,所以必须使用凝聚的访问方式以得到最大的设备存储器带宽[8]。GPU可以使用一个指令将32 b,64 b或128 b的字由设备存储器读入到寄存器中,而半个warp块中的线程必须以凝聚的访问方式读写设备存储器。
通过测试设备存储器峰值带宽,凝聚的64 b的访问性能略强于32 b字的访问,而在每个线程块具有128个线程的情况下,32 b字和64 b字访问的峰值带宽分别是121.6 GB/s和124.8 GB/s(注意增加_syncthreads()指令可以略微增加带宽性能)。在矩阵乘法的实现中采用每个线程块具有128个线程的设计,这样基本上可以保证设备存储器访问的最高性能。
2 SGEMM实现和性能测试
上述程序优化原则可以被用于单精度矩阵乘法(SGEMM)实现。SGEMM属于计算密集型的任务,但在设计不当的情况下,会受到设备存储器带宽的约束。
计算矩阵A和矩阵B的乘积得到矩阵C的计算可以按照下述分块的方式进行:每一个线程块负责计算矩阵C的一个子块Csub,而该线程块中的每一个线程负责计算子块Csub中的若干个元素。Csub是两个长方形矩阵的乘积:矩阵A的和Csub具有相同行数的长方形矩阵乘以矩阵B的和Csub具有相同列数的长方形矩阵,如图2所示。
由于GPU硬件存储的限制,这两个长方形的矩阵也被分割为相应的子块,称为Asub和Bsub,Csub为对应Asub和Bsub子块矩阵乘积的和。Csub的大小决定了每轮运算必须被加载到GPU芯片的数据量大小。如果Csub过小,则SGEMM会变成一个带宽受限的问题。
假设Asub,Bsub和Csub的大小分别为m×k,k×n和m×n,而矩阵A,B和C分别被分为M×K,K×N和M×N个这样的子块。计算一个Csub需要从设备存储器中取K个Asub和Bsub。由于在矩阵C中有M×N个Csub,则整个计算需要从设备存储器中取M×m×N×n×K×k×(1 m+1 n)个浮点数。
该算法以行优先的方式设计。线程块的线程数目为128,被组织成4×32的格式。Csub的大小为16×256。每一个线程块负责计算一个Csub,而一个线程负责计算Csub中的2列。图3显示了每一个线程的流程图。
首先,线程块申请16×32的共享存储器空间用于存储矩阵A的子块,每一个线程申请寄存器用于存储Csub中的16个float2即32个float。在使用线程ID计算出访问矩阵A,B和C的指针位置之后,进入一个循环。在每一轮循环中,一个线程块从设备存储器中读入16×32的Asub的数据到共享存储器中,之后,又一个内层的循环被执行:在每一轮内层循环中,一个线程从矩阵B中读入一个float2,把这两个float和共享存储器中相应的数据做计算,将结果累加到Csub对应的寄存器中。最后,每个线程将其计算的2列Csub的数据写回到设备存储器中。
图4对本文的代码和CUBLAS代码的性能进行了比较。当矩阵大小为8 192×8 192(在GTX285 GPU上可以进行计算的最大尺寸)时,本文的代码达到了393 Gflop/s的性能,比文献[8]中实现的性能大约高5%。而使用共享存储器进行MAD计算的理论指令峰值吞吐量为410 Gflop/s。本文的实现和文献[8]中的SGEMM实现具有类似的代码结构,下列因素导致本文的代码取得了更优的性能:
(1)在本文的代码中,Csub的尺寸为16×256,而非16×64,从而减少了从设备存储器访问数据的数量;
(2)在本文的代码中,一个线程块包含128个线程,而每一个线程以凝聚方式访问一个float2即64 b字,这样可以取得较高的设备存储器带宽;
(3)在本文的代码中,一个Asub的大小为16×32,而非16×16,从而使得整个算法的同步指令次数被大大减小了。
3 FFT实现和性能测试
另一个代码优化的示例为批量的1维复数FFT,根据向量的大小不同,1维复数FFT的实现可以被分为三种情况:当向量的长度为8或16时,不需要进行线程之间的数据交换,一个向量的计算可以在寄存器中被一个线程完成;当向量的长度为64,256,512,1 024,2 048或4 096时,线程之间需要进行数据交换,且数据交换可以在共享存储器中完成;当向量的长度超过8 192时,必须通过设备存储器进行矩阵数据的转置,这可以通过调用两次内核代码完成。因此,当向量的长度为8或16时,此计算为存储带宽受限型;当向量的长度为2 048或4096时,此计算为通信受限型。在此,取向量的长度为16,2 048或4 096两种情况来说明程序优化的方法。
3.1 向量长度为16
当向量长度为16时,这是一个存储带宽受限型的问题。为了得到良好的性能,必须要使得对设备存储器的访问达到一个良好的带宽。在进行计算之前,原始数据被首先进行整理,使得线程可以以凝聚的方式,读取一个向量中需要计算的16个元素。尽管这种数据整理可以通过在共享存储器中的数据转置来代替,但在本文不涉及到此问题。所以在这个实现中,一个线程需要做的就是从设备存储器中读取16个元素,进行1维FFT计算,再把数据写回设备存储器。由于每一个元素为复数,所以每一个线程自然的就在读取64 b字的数据。此外,根据前面提到的优化原则,在每一个线程块中包含256个线程,并且在计算和读写设备存储器之间加入同步指令,以得到最优的设备存储器访问带宽。
3.2 向量长度为2 048或4 096
当向量长度为2 048时,每一个线程块包含128个线程,每一个多处理器上可以具有2个活动线程块。每一个线程负责计算2 048个复数元素,这通过在寄存器中计算长度为16,16和8的向量完成,在3次计算之间,通过共享存储器交换数据。为了在共享存储器中交换2 048个单精度复数数据,至少需要8 KB的空间用来交换复数的实部或虚部,使得一个多处理器上只能存在1个活动线程块。前面的优化原则提到,一个多处理器上至少需要2个活动线程块以掩盖同步指令的开销。本文的解决方法是,通过以下操作将共享存储器的使用量减少一半:首先前一半的线程先将自己的数据放入到共享存储器中并同步,各个线程将自己需要的数据取回,而后后一半的线程再将自己的数据放入到共享存储器中并同步,各个线程将自己需要的数据取回。这样,一个多处理器上可以同时具有两个活动的线程块。在此实现中,通过控制Warp ID来达到此种操作效果。当向量长度为4 096时,实现的思路与向量长度为2 048是类似的,只不过此时一个多处理器上只能具有一个活动的线程块。
3.3 性能测试
对64 MB的数据进行批量1维FFT的性能测试:当向量大小为1 024时,对应的信号数目为8 129;当向量大小为2 048时,对应的信号数目为4 096,等。图5比较了本文的代码和CUFFT 1.1,Volkov和Kazian的代码[5,6]的性能。注意,在本文的实现中,当向量大小为2 048和4 096时,该计算通过一次内核调用完成,从而使得数据不需要进行额外的整理,而其他一些实现则要求数据在GPU计算完成之后使用CPU在主存中对数据进行整理。
4 结语
书写高性能的GPU代码是十分困难的,本文介绍的工作并不是为了得到特定算法的高性能代码,而是为了总结帮助程序员得到高性能GPU代码的程序设计原则。程序员在进行算法设计时,往往需要进行自顶向下和自下向上的反复优化。本文中提到的优化原则可以作为对CUDA程序设计指南和其他一些工作的补充。在硬件体系结构发生变化的情况下,原先非最优的代码可能又会取得顶级的性能。取得特定算法的高性能代码本身并不是本文所看重的,而在进行代码优化过程中取得的经验,有助于对GPU和GPU集群程序设计问题进行进一步的研究。
参考文献
[1]吴恩华,柳有权.基于图形处理器(GPU)的通用计算[J].计算机辅助设计与图形学学报,2004,16(5):601-612.
[2]吴恩华.图形处理器用于通用计算的技术、现状及其挑战[J].软件学报,2004,15(10):101-112.
[3]苏超轼,赵明昌,张向文.GPU加速的八叉树体绘制加速算法[J].计算机应用,2008(10):51-62.
[4]储璟骏,杨新,高艳.使用GPU编程的光线投射体绘制算法[J].计算机辅助设计与图形学学报,2007,19(5):257-262.
[5]朱亚平,杨慧珠,董渊,等.OpenGL技术在地震数据可视化中的应用[J].石油地球物理勘探,2000(8):35-42.
[6]徐品,蓝善祯,刘兰兰.利用GPU进行通用数值计算的研究[J].中国传媒大学学报:自然科学版,2009(4):101-112.
[7]NVIDIA Corporation.NV1DIA CUDA programming guide Ver sion 0.1[EB/OL].[2009-04 28].http://developer.nvidia.con/cuda.
[8]COPPERSMITH D,WINOGRAD S.Matrix multipli cation viaarithmetic progressions[C]//Proceedings of the nineteenth an nual ACM symposium on Theory of computing.New York,NY,USA:ACM,1987:1-6.
[9]RYOO S.Optimization principles and application performanceevaluation of a multithreaded GPU using CUDA[C]//Pro ceedings of 13th ACM SIGPLAN Symposium on Principles andPractice of Parallel Programming.[S.l.]:ACM,2008:73-82.
[10]VOLKOV V,DEMMEL J W.Benchmarking GPUs to tunedense linear algebra SC08[C]//Proceedings of the 2008 ACM/IEEE Conference on Supercomputing.[S.l.]:IEEE Press,2008:1-11.
[11]GOTO K,VAN DE GEIJN R A.Anatomy of high performancematrix multiplication[J].ACM Transactions on MathematicalSoftware,2008,3(43):12-25.
FFT算法 篇7
设长度为N的有限序列x (n) 的DFT为:
设序列x (n) 的长度N (N=2M, M为任意整数) , 根据N的奇偶性把x (n) 分解为两个N/2点的子序列如式 (2) 所示:
则x (n) 的N点DFT为式 (3)
由于
所以
式 (5) 中, G (k) 和H (k) 分别为g (m) 和h (m) 的N/2点DFT, 表达式如下 (6) (7) 所示:
由于G (k) 和H (k) 均以N/2为周期, 考虑到对称性, 有WÂÁÁÂÁÂ=-WÂÁ, X (k) 就可以用 (8) (9) 表示:
通过将N点的DET分解和 (8) (9) 的运算。可以有效减少DFT的运算量。由于N是2的正整数次幂, N/2点的DET可以继续进行分解, 直到分解为2点的DET为止。这种运算称为蝶形运算符。
2 算法分析
本设计实现了在Altera FPGA上实现的基于Open CL的快速傅里叶算法 (FFT) , 该基准程序并行处理32点复杂单精度浮点数据。数据以有序输入, 以倒位序输出。
本设计能够每时钟周期处理8个数据点单基数2的FFT引擎。根据现有设备资源, 设计还可以额外例化模块, 以实现更高的性能, 也可以定制变换的大小。
FFT算法模块被配置为单个工作项内核, 它使用一个滑动窗口来表示延迟元件。该应用程序从全局内存中读取输入数据, 并将数据馈送到FFT单元然后将存储结果保存到全局内存。性能将受到内存带宽的影响。
3 性能分析与比较
本设计基于64位Windows 7的操作系统, Intel酷睿i3四核四线程的处理器, 以及DE1-So C FPGA开发板。对于处理不同的数据量在不同硬件环境下的运行情况如表1所示:
数据是10次测试结果的平均值, 从表中的结果可以看到, 对于相同大小的数据量, FPGA的加速效果比Cortex A9有很大的提升, 并且, 随着数据量的加大, 效果也越来越明显。此外, 在保证处理不同规模问题质量的前提下, 基于FPGA平台的Open CL设计也很大程度的降低了系统的功耗。
4 结论
本设计使用的SDK为Altera发布的支持Open CL并行编程的SDK, 本设计简单有效的解析了Open CL和FPGA的并行体系架构, 使用DE1-So C平台和FFT算法为实验对象, 对并行的加速算法Open CL在FPGA上实现的原理和方法进行研究。与传统的串行平台作对比, 测试加速性能和功耗情况。经过对数据分析对比, 由FFT算法的实验数据结果可见, 在FPGA上实现开放计算语言Open CL标准并且与Open CL并行编程模型和拥有强大并行体系结构的FP-GA相结合, 不仅可以更大程度的提高性能, 也能更大程度的有效降低系统功耗。
参考文献
[1]陈国良, 孙广中, 徐云, 龙柏.并行计算的一体化研究现状与发展趋势[J].中国科学.2009, 54 (8) :1043-1049.
[2]陈国良.并行计算:结构·算法·编程[M].北京:高等教育出版社, 2003:368-402.
[3]詹云, 赵新灿, 谭同德.基于Open CL的异构系统并行编程[J].计算机工程与设计, 2012, 33 (11) :4191-4195.
[4]苏金国, 李璜, 杨健康等译, Open CL编程指南[M].北京:机械工业出版社, 2013:75-77.