并行之美基于DSP的FFT算法并行分析

基于DSP的FFT算法并行分析1摘要—快速傅里叶变换FFT是一种应用广泛的实现离散傅里叶变换DFT的快速算法，研究FFT的高速实现有着非常重要的意义和价值。采用通用DSP来实现FFT，可以满足高速实时处理的要求,但在高速信号处理系统中，一些复杂的信号算法通常需要多次FFT运算才能求得结果，因此缩短FFT运算时间具有实际意义。本文讨论了FFT算法的内在并行性，如何将FFT算法进行并行化，需要对算法的基本结构进行分析，这里主要以频率抽取基2算法为例，对FFT的并行性和串行性进行分析。关键字—并行处理技术，基于DSP,FFT变换I.引言散傅里叶变换(DFT)是数字信号处理中最基本、最重要的运算，许多算法，例如卷积、滤波、波谱分析等都可以化为DFT来实现。[1]但是DFT的计算量相当大，1965年Cooley和Tukey提出了快速傅里叶变换算法(FFT)，大大提高了DFT的计算速度。快速傅立叶变换算法主要是利用了旋转因子的对称性、周期性和还原性，从而将N点的有限长序列的DFT分解成几个较小点数的DFT。[2]由于DFT的计算量与N平方成正比，所以长度N越小的序列其DFT的运算量越少，从而提高了运算速度。快速傅里叶变换FFT广泛应用于雷达、数字通信、图象处理、语音信号处理和生物医学等领域，且在高速信号处理系统中，FFT运算已成为一个基本运算，一些复杂的信号（如小波变换）算法通常需要多次FFT运算才能求得结果，所以研究FFT算法的并行性，缩短FFT算法运算时间具有很大的现实意义。[3]本文主要以频率抽取基2为例分析FFT的串行性与并行性，研究基于4TS201DSP硬件平台的并行FFT算法设计。II.FFT的并行化A.频率抽取FFT算法的并行性对于一个N点DFT(离散时间傅立叶变换)，其变换定义为),1,1,0()()(210NjNNnnkNeWNkWnxkX（1-1）FFT算法的本身，即是将N点DFT进行分组，以基2FFT算法为例，先将点DFT分成两个N/2点DFT，再分为四个N/4点DFT，进而八个N/8点DFT，至分为N/2个两个DFT。例如频率抽取(DIF)基二FFT算法，是将DFT定义公式中的x(n)按序号分为上下两部分:nkNnkNNnNNnnkNNnnkNWNnxWnxWnxWnxkX)]2/()([)()()(2/1/2012/1/20（1-2）式中k/2)-1(nkNW，分别令k=2r,k=2r+1，而r=0,1,...,N/2-1。再令)2/()()(Nnxnxng，nNWNnxnxnh)]2/()([)(，即可得到由输入数据进行分组的迭代公式:12/02/12/02/)()12()()2(NnnrNNnnrNWnhrXWngrX（1-3）从上述分析可知，对于频率抽取基2类FFT算法，其算法本身的思想即是将输入数据在每一级蝶形运算是按照序号将输入数据按照上下两部分进行分组。因此，可以按照分而治之优化原则对蝶形运算的输入数据进行分割，将运算量分配到多个计算节点。同时，将N点DFT分为相互独立的N/2点DFT并向下细分成N/2个两点蝶形运算，每一个蝶形运算和DFT都是相互独立的，可以并行计算。B.时间抽取FFT算法的并行性和频率抽取FFT算法相比，时间抽取FFT算法的数据分割原则有所不同。首先按照x(n)的序号的奇偶进行分组，有:12/0)12(1/202)12()2()(NnkrNNnrkNWrxWrxkX（1-4）令1/202/)2()(NrrkNWrxkA，1/202/)12()(NrrkNWrxkB则有)12/,1,0)(()()(NkkBWkAkXnN(1-5)其中A(k)和B(k)都是N/2点DFT，而X(k)则是N点DFT。按照奇偶分割原则继续往下划分，即可得到完整的时间抽取FFT运算蝶形图。按照以上分析，无论频率抽取FFT基2算法还是时间抽取基2算法，其基本思想都是对输入数据进行分割，类似的，基4算法和分裂基算法，都是在数据分割方法上对于最基本的基二算法进行的改进。这一类的FFT算法现在均有非常成熟的串行程序，利用分治原则和算法内部数据的可分割性，可以方便将输入数据进行分割处理，将运算量分摊到多个计算节点，从而设计出相应的并行算法。从理论上将，除了数据的分割需要耗费处理器一定计算量以外，剩余的计算量是将串行算法的计算量进行平均分配而无需额外的开销，因此按照输入数据分割方法设计的并行算法拥有优异的加速性能。C.FFT算法的串行性上文分析了FFT算法内在的并行性，可以通过分治策略将属于数据进行分割，从而达到在多片DSP中的并行算的运算负载分配。[4]虽然FFT算法天然地具备了这一并行特性，但是算法内部流程也具有无法并行化的部分，在算法并行之美——基于DSP的FFT算法并行分析离基于DSP的FFT算法并行分析2的流程层面上具备固有的串行性。在FFT算法当中，旋转因子是蝶形运算中的重要内容。因此，在蝶形运算之前必须先得到旋转因子NjNeW2。根据实际应用需求，旋转因子可以通过查表或者实时计算生成。此外对于频率抽取基2FFT算法，输入数据为顺序排列，蝶形运算后输出数据为变序输出，因此需要进行码元翻转才能成为正常的输出序列，因此计算顺序可分为3步:图1-1FFT算法串行流程示意图由上图可以看到，FFT算法在流程层面上具有很强的依赖性，旋转因子是蝶形运算的必要输入数据，蝶形运算的实质是将两个时域点数据依旋转因子加权之和。[5]对于DIFFFT算法而蝶形运算的输出则是逆序的，需要通过码元翻转才能得到正确排序的频域数据，因此这三个流程以上一级流程的输出作为输入，具有不可改变的内在串行性，在该层面上难以通过算法并行化策略进行并行化处理，只能采用块分配策略(BlockArranging)进行并行化，具体方法将在下文阐述。III.粗粒度块分配流水线并行算法实现容易想到，将一个具有内在线性的算法，可以按照其线性顺序关系将处理任务分解为一个顺序的流水线，这一流水线是在作业级将算法进行分割，因而属于粗粒度并行算法，这一并行处理机的类型被称为单一程序多数据并行处理机(SPMD).在理想的情况下，一条k段的线性流水线能在k+(n-1)个周期内能处理n个任务，因此总的并行处理时间为:)]1([nkTpar(1-6)其中parT代表并行程序运行时间，定义为单步程序执行时间，k为每一个流水线流程上的程序指令长度，n为流水线上任务的个数。同样地，相同算法的串行执行时间为:knTsen因此，加速比为:1nknkTTSparsen(1-7)通用DSP虽然在处理速度上不如专用DSP，但内存容量却要大得多。[6]各种并行快速算法的提出，其基本思想主要是将以前的较复杂的串行算法分解为可以并行执行的、相对简单的运算，或将标量运算转化为向量运算，或是采用复杂的松弛(Relaxed)算法等，都是依靠较大的内存开销来换取并行处理速度的提高，因此提高并行处理速度就可以采用充分利用多个通用DSP的内存来达到。[7]面向MIMD的并行算法设计的首先是将待处理的输入数据分配给系统中的各个DSP，比较常用的是采用块分配(BlockArranging)的形式，DSP对数据的处理形式也因此被称为块并行处理策略(ParallelBlockProcessingScheme)，其基本思想是将输入数据分割为长度相等的块，再分配给不同的DSP，每个DSP单独计算和输出一块数据，以充分发挥单个DSP的运算性能。[8]同时，块并行处理策略还相应减少了DSP之间通信的额外开销，减小通信对数据处理速度的影响。A.流水线块分配算法的并行化方法图5-2给出了块分配并行处理的调度模型，根据硬件系统结构，我们安排待处理的原始数据由FPGA存入大容量外部SDRAM，各片DSP通过共享64位总线从SDRAM中以DMA方式读取待处理数据。[9]在这里定义idmaT为DSP#i从SDRAM中获取数据的DMA传输时间;当第一批数据到达SDRAM,DSP#1立即开始处理进程，而DSP#2,DSP#3,DSP#4则处在等待数据的空闲状态;定义从第一批数据到达SDRAM和第i批数据到达SDRAM之间的时间差为iidleT，这一时间也就是DSP#2的空闲时间;定义ibfT为DSP#i进行递归蝶形运算时间，itwT为DSP#i进行码元翻转运算时间;定义itotalT为DSP#i对一批数据的总的处理时间。定义串行算法处理i批次数据所采用的总共时间为:isenT图1-2粗粒度FFT算法流水线实现时问一空问拓扑图如上图所示:对于DSP#1处理第一批数据的总时间:11111twbfdmaouttotalTTTTT(1-8)该时间等同于采用串行算法处理一批数据的时间，即:itwibfidmaiouttotalsenTTTTTT11(1-9)而串行算法处理4个批次所需要的总时间为基于DSP的FFT算法并行分析3)(11114twbfdmaoutsenTTTTT(1-10)而利用块分配原则设计的流水线并行处理方法，处理4个批次的数据时间为:414totalidletotalTTT(1-11)当等待4个批次数据传输时间04idleT的时候，加速比为4.另一方面也可以看到，如果仅仅只对一个批次数据进行处理，1total1senTT，采用块分配流水线并行算法并不能提高效率。如果对于单一批次的数据，可以采用数据分割的方法通过分治策略对算法进行优化，使单输入的数据变为多输入的数据进而进行并行处理。其次，影响加速比的总要因素为iidleT，即系统对于数据到来的等待时间决定了加速比。在设计这一类块分配的并行算法的时候，对于单一节点的运算负荷即时间1tw1bf1dma1out1totalTTTTT不需要做太多考虑，因为这一时间不影响系统的加速比，而应该尽量减少节点之间因为等待数据传输而产生空闲时间。块分配流水线并行算法的另外一个优势为各个计算节点的程序可以是一样的，而且可以直接在不改变串行算法的结构的情况下实现并行计算，因此工程实现最为方便。B.针对硬件系统优化的进一步算法重构并行程序的硬件实现和一般并行算法不同点在于，提高系统运行的一个重要手段在于从硬件结构入手，针对硬件平台的实际情况对算法进行优化和重构。由于TS201DSP内部存储器针对顺序读取进行了优化，虽然片内的高速缓存对随机读取进行了优化，然而由于高速缓存容量有限，对于2048点以上的数据就难以存入高速缓存。[10]从图4-2可以知道对于常规的FFT运算来说，每一级蝶形运算步进是双倍的，这一读取就是非顺序的，如果才用常规的FFT运算结构，由于随机读取的低效率会引起运算性能的大大下降。[11]为了进一步优化算法性能，现将频率抽取FFT算法的结构进行重构如图1-3。从图中可以看出，重构后输入和输出顺序和频率抽取FFT算法一致，但是每一级蝶形运算之间的运算步进不变，仍然是两个相邻数据的运算，因此DSP经过优化的存储器顺序访问性能可以使硬件实现的性能大为提高。[12]另外可以得到，对于这一改进的FFT运算，其输出数据序号有如下关系:为奇数221为偶数2)(nNnnnnf(1-12)对于硬件系统来说，当n为偶数是只需将序号n右移一位就可得到序号f(n)，当n为奇数时，将输出数据的序号n最高位置1即可得到f(n)，这一操作等效为1为右旋操作。由于ADSPTigerSHARCDSP支持位反序寻址，因此对于输出数据的码元翻转只需进行位反序寻址即可调整为顺序输出，其效率非常高。图-3重构后FFT蝶形运算结构图从图1-3的重建后蝶形运算结构图可以看出，对蝶形运算输入数据的顺序在不同的阶段是不同的，因此在不同阶段的蝶形运算过程中采用乒乓缓冲来实现数据的交换和移动。C.块分配流水线并行算法的程序实现图1-4为处理节点内部程序框图。由上述分析可以，每一个计算接点内部出题流程相同，因此只画