基于超标量处理器的高效FFT映射方法

高立宁，朱亮，刘腾飞，刘峰

(北京理工大学 信息与电子学院，北京 100081)



基于超标量处理器的高效FFT映射方法

高立宁，朱亮，刘腾飞，刘峰

(北京理工大学 信息与电子学院，北京 100081)

针对超标量处理器的结构特点，研究新的映射方法，实现高效FFT运算. 对现代超标量结构处理器进行建模，分析FFT算法在其上执行情况，得出内存访问是FFT算法执行的关键点. 并进一步对FFT的内访问过程进行建模分析，最终实现了一种基于cache优化的高效FFT映射方法，该方法将FFT进行拆分实现，充分发挥了cache的作用，进而提高了处理性能. 最后在ADI公司的TS201数字信号处理器上，以该映射方法为指导实现了基2FFT算法，实验结果显示在处理点数超出cache容量时，本映射方法可以大幅度提高处理性能.

快速傅里叶变化(FFT)；高速缓存(cache)；超标量处理器

快速傅里叶变化(fast Fourier transform， FFT)是现代化雷达信号处理中的关键技术之一，由于雷达系统是一种强实时处理系统，FFT作为系统中的重要组成部分，必须在限定的时间内完成处理，因此高效的FFT处理是必要的[1].

目前高效FFT算法实现主要有两种：一种是通过ASIC或FPGA实现，这种硬件实现方法适应算法特点，功耗低且实时性高，但是这种方法缺乏灵活性且开发周期长. 常用的另一种方法是基于处理器的编程实现，该方法通过程序代码将FFT实现，开发周期短且灵活，并且处理器的主频在不断提高，FFT的执行时间也在不断缩短，因此基于处理器程序编写的方法有着更多的优势. 然而处理器结构固定，并非完全按照FFT算法结构特点设计，这就要求在编程实现时，找到一种新的映射方法将FFT映射到处理器中，使处理器的处理效率最大化发挥，进而满足实时性需求.

针对上述分析，本文对超标量处理器结构进行研究，将与处理相关的因素进行参量化，建立了超标量处理器的程序执行模型. 以基2FFT算法为例，分析了常规FFT算法映射方法的缺点，并针对这些问题提出解决方法，提出了可拆分的FFT映射算法，该映射算法充分利用了处理器上的处理资源，同时利用cache特点大大缩短了内存访问时间，进而使处理器性能最大发挥. 最后以该映射算法为指导，在ADI公司的TS201数字信号处理器(DSP)上实现了基2FFT算法，并取得了满意的结果.

1 超标量处理器的结构特点

现代化处理器中均采用了超标量流水线技术，超标量流水线是并行流水线，其将指令执行分为若干部分(通常包括取址、译码、执行、访存、写回)，然后通过流水线技术，使同一时刻多条指令处于不同的处理阶段. 其利用时间并行方法，提高了处理器的吞吐率. 并且处理器中采用不同的执行单元，执行不同功能指令，可以在每个机器周期中启动对多个指令的处理，提高了指令的执行效率[2-3]，图1显示了超标量处理器的工作原理.

图1所示为同时有4条指令执行的超标量流水原理图，从图中可以看出指令执行被分为6个阶段：① 读取阶段，主要负责从内存中获取执行代码；② 解码阶段，主要任务是将指令翻译为机器可识别的操作码；③ 调度阶段，负责将指令按照功能分配到各个功能执行部件上执行，同时在该阶段负责查找指令中的相关指令，并进行调度；④ 执行阶段，主要完成指令的功能；⑤、⑥完成和退出阶段，这两个阶段主要负责修改指令完成后的机器状态，同时保证指令的顺序完成. 在超标量架构的处理器中，如上处理中阶段①，②，③，⑤及阶段⑥针对所有指令工作机理一样，因此在设计中常常通过集中处理方式完成这⑤个阶段，针对阶段④，不同功能指令执行过程会不同，因此这一阶段设计采用了分布式处理方式[2-3].

以上介绍均是假设流水处理中每个阶段执行时间相同，然而实际应用中由于内存访问速度远远低于内核速度，导致数据访问时间很长，这种情况会打破流水线运行，降低处理效率. 为解决这一问题，处理器中引入了同内核速度相当的高速缓存(cache). cache利用程序的局部性原理，将最近访问的数据存储在cache中，当访问内存时首先会查询cache中是否缓存，如果在缓存命中(hit)，那么直接获取数据，减少时间开销；如果没有命中(miss)，将会启动内存访问，并将访问结果缓存[4-6]. 通过合理设置cache容量，可以有效缓和内存速度与处理器内核速度不匹配问题，充分发挥处理器性能.

综上所述，基于超标量结构的现代化处理器，通过并行流水增加处理器吞吐量，并提高处理效率；为了缓和内存与内核的不均衡，现代处理器中引入cache，使内核计算与数据访问相匹配. 通过如上技术改进，可以不断提高处理器的处理能力. 在实际应用中程序执行的实时性，一方面与处理器的处理能力相关；另一方面也与应用算法在处理器上的映射相关，需要采用合适的映射算法，使程序执行过程适合超标量架构特点，最大化发挥处理能力，才可实现实时处理.

2 基于超标量处理器的FFT映射方法

如上讨论可知，现代处理器中不同类别指令分别由独立的执行单元运行，因此处理器可以抽象为若干相互关联的执行单元的集合，设处理器中执行单元数为s，其中任意一个单元i的执行速度为vi，那么所有执行单元的处理速度(对于不同执行单元处理速度定义也不同，内存控制执行单元以访问带宽表示；内核计算单元以单位时间内执行的操作次数表示)可构成集合V={v1，v2，…，vs}. 运行于处理器上的程序，也可按照执行单元划分为s个处理任务(任务包括内存读写、乘法操作等)，设其中任意一个处理任务i的任务量(对于不同任务的任务量定义不同，内存访问以访问的字节数表示，数据计算以运算次数表示任务量)为ci，那么程序的任务量可由集合C={c1，c2，…，cs}表示. 当程序在各个部件执行时，每个单元的执行时间可组成集合

T={t1，t2，…，ts}={c1/v1，c2/v2，…，cs/vs}.

集合T中每个元素表示程序中各个任务的执行时间，这些执行相互交叠，构成了程序最终处理时间. 因此为研究实时性，应对集合T中的元素间关系进行研究，为方便讨论符号定义如下：

①αi|j表示任务j的执行时间与任务i执行时间重叠部分占任务j执行时间的百分比，这里称该百分比为重叠度(使用重叠度，主要反映各个执行单元间的并行度，该比值与映射算法及处理器资源相关)；

②t(i，j)表示任务i，j的总执行时间.

由如上定义可得处理过程中任务i，j的总执行时间为

(1)

如果将任务i，j合并为一个大的任务I，那么该任务与任务集C中的任务k并行执行时间可表示为

(2)

对如上计算过程进行递归计算，那么可得程序中s个任务并行执行总时间为

(3)

式(3)表明了执行过程中，程序每部分任务处理时间与总时间关系. 由于集合T中的各个元素具有对等性，不失一般性，这里令t1为所有执行时间中数值最大的，即满足如下:

(4)

这里称t1为最大执行时间，那么可得如下关系:

(5)

由式(4)和(5)可知当所有重叠度均为1时，执行时间最小，此时的物理意义为：程序在处理器上各执行单元的处理过程均隐藏在执行时间最长的处理过程中. 因此应用算法映射到处理器上执行时，应结合处理器特点合理安排处理顺序，尽量使所有处理均能隐藏在处理时间最长的任务中，即要求程序在实现时应充分发掘指令间的并行性，合理安排指令流水，使程序执行时间趋于最大处理时间.

基于如上分析，对程序在现代处理器上的执行时间建立了模型. 当FFT算法映射到处理器上执行时，按照如上所述可划分为4个处理任务：乘法任务、加法任务、减法任务和读写任务. 针对现代化处理器，如上4个任务可分别在乘法器、加法器、减法器和内存访问控制器中执行. 对指令流水进行合理安排，FFT执行时序如图2所示(该图所示时序图为处理数据一部分位于cache，一部分位于内存).

流水建立分为填充、流水、排空3个阶段[6]. 当进入流水阶段后，所有的处理(乘法运算、加减法运算)均隐藏在内存访问中，只有在填充和排空阶段部分隐藏. 因此可以把内存执行时间作为最大执行时间，即任务1，令乘法处理为任务2，加法处理为任务3，减法为任务4，设处理中每级FFT的蝶形计算个数为n，那么可得每级处理中的重叠度为

(6)

式(6)表明FFT处理中各个任务的重叠度与处理中每级蝶形运算个数n成反比，当n取值较大时，可认为所有任务执行的重叠度为1，即所有内核计算隐藏在内存访问中进行. 对于雷达系统而言，FFT处理通常是在几千点以上，每级的蝶形计算次数也在几千个以上，因此内存访问便成为FFT映射中的关键性问题. 通过第1节介绍可知，内存访问速度低于内核处理速度，故引入了缓存解决该问题. 通过图2中所示时序图可以看出，当数据在cache中时，内存读或写的时间与内核计算时间相匹配，但是访问数据没有命中时，会在流水线中引入额外的时间消耗[7-9]. 可见cache的利用对FFT程序执行很重要. 如下主要针对FFT在处理器中的内存访问过程进行建模并分析，首先对一些相关因素符号化：

① cache命中率设为P，该值表示数据访问过程中，数据在cache中的概率；

② 丢失时引入的额外时间设为tmiss，当访问数据不在cache中，需从内存获取数据而产生的额外消耗时间；

③ 高速缓存的存储容量设为cacheL，为方便讨论以复数点为单位，其中每个复数点包括IQ两个数据，每个数据由b个字节组成;

④ 高速缓存一行可以缓存的数据容量为w，为了讨论方便这里w表示的是复数点，其中每个复数点包括IQ两个数据，每个数据由b个字节组成；

⑤ 数据的访问时间Tio，该值表示数据的读写时间；

⑥ 处理的数据长度为L，该值以复数点为单位且为2的整数次幂，其中每个复数点包括IQ两个数据，每个数据由b个字节组成；

⑦ cache访问带宽为Bc，单位为byte/s.

目前常用的FFT算法大部分是基2的快速变换算法，因此本文中以基2算法为例对内存访问过程建模，所得结果对其他快速变换算法同样适用[10-11]. 基2FFT算法的蝶形数据流图如图3所示.

图4为蝶形运算的传统实现方法流程图.

依据图4可知，传统的FFT映射算法需要做n级(n=lbL)蝶形运算，每一级蝶形运算访问的数据量为L，因此，如果每一级访问的cache命中率为Pi，时间消耗为tmiss，那么FFT的执行时间中因cache丢失所引入的时间消耗为

(7)

有前面论证可知FFT执行中主要时间为数据访问(包括数据读、写)，由于每级处理数据访问量2L(包括读一次，写一次数据访问)，根据前面的参数假设，可得处理长度为L的数据时内存访问时间(传统FFT中有效执行时间)为

(8)

综上所述，传统的FFT映射方法的执行时间近似为

(9)

由于基2FFT数据流为原位处理，所以数据写回内存时不会发生丢失，故每级处理中丢失次数为L(1-Pi). 可用如下公式表示正常映射方式的内存访问时间:

(10)

在式(10)的表示含义为：处理过程中的第一级处理必然存在丢失现象，由于内存控制器以cache的行长w为单位访问数据，那么访问w个数据后便会丢失一次，所以第一级丢失概率为1/w；如果处理长度小于cache容量，以后每级处理均可在cache中命中数据，丢失率为0；如果处理长度大于cache容量，以后的每级都与第一级处理过程一样，每级的丢失率为1/w.

为了充分发挥cache的作用，缩短大点数处理时间，应提高访问过程中的命中率，减少丢失次数. 一种解决方法便是考虑映射实现时将大点数划分为小点数分段处理，这样每次小点数处理时可以充分发挥cache作用，缩短了处理时间[12]. 为实现这个目的可以从离散傅里叶变换公式入手.

(11)

对式(11)中的变量l，k分别进行分解：

(12)

将式(12)带入式(11)中可得

(13)

式(13)中第二项求和表示按列方向进行离散傅里叶变换，那么通过式(13)可以获得拆分的方法：

① 将一维的数据序列分为一个F*Y矩阵；

② 对矩阵的每一列进行傅里叶变换；

④ 再逐行进行离散傅里叶变换.

基于如上思想FFT算法映射到处理器执行时，可以先拆分，然后对每维采用FFT算法进行处理，为了缩短处理时间，可将第3步融合到每列FFT计算的最后一级. 通过该映射方式映射后的程序执行流程图如图5所示.

采用本文介绍的映射方法实现快速傅里叶变化，可根据式(12)计算每步执行时间.

对二维矩阵所有列FFT变换的时间消耗：

(14)

因子相乘，可与FFT中最后一级的运算放在一起进行，那么乘以因子的时间消耗主要是内存读时间，可表示如下：

(15)

对矩阵所有行FFT变换的时间消耗为

(16)

按照如上3个步骤执行完后，总时间消耗为

(17)

按照将大点数数据拆分成小段去运算的思路，使得每一小段能够放到cache中，将cache充分利用起来，提高了cache命中率，进而总运算时间得以下降. 如式(17)表示，采用本映射算法映射后的FFT算法，在处理数据长度小于cache容量时，并没有优势，但是当处理长度大于cache容量时，要优于正常映射方法. 为了验证本映射方法的有效性，如下以该映射算法为指导，在ADI公司的TS201处理器上实现基2FFT算法.

3 基于TS201处理器的高效FFT映射实现

TS201处理器是ADI公司推出的一款高性能静态超标量信号处理器，主频为600 MHz，采用超标量流水线技术每个时钟周期可执行4条指令，在内核计算方面采用SIMD(单指令多数据流)方式，每个时钟周期可以执行2个浮点乘法，2个浮点加法和2个浮点减法；存储性能方面，处理器内部集成了24 Mbit内存，共分为6个分区(bank)；TS201内部有4条128 bit宽度的内部总线，通过J和K两个整型ALU可以实现一个周期内4个32 bit数据的读取和4个32 bit数据的存储；为解决内核速度与内存速度不均衡问题，每个bank带有一个容量为128 kbit(2 048个浮点型复数)的高速缓存，每个cache的行长度为256 bit(4个浮点型复数)，当数据位于高速缓存中，一个时钟周期便可完成数据访问，如果数据不在cache中，将会花费7～8个时钟周期完成数据访问任务.

按照映射算法将FFT映射到TS201上执行，首先应确定拆分后的二维矩阵大小. 这里以使两维处理均衡为目的，按照如下准则划分：

(18)

在实际应用中由于处理的数据量很大，常常将数据存储在DSP片外存储空间中(通常为SDRAM)，因此在从SDRAM中访问数据时可以使用DSP处理器中的DMA按照式(18)的划分方式将数据搬移到片内存储空间. 由于第一步处理是按列进行，所以搬移可按照列连续，行跳跃的方式读取数据，目前所有DMA均支持该方式. 数据获取完毕，便可访问数据进行处理，在处理过程中为了充分发挥处理器的存储资源，本程序将旋转因子和数据放在不同bank中，这样便可在读取处理数据的同时，访问旋转因子. 另外对于程序中旋转因子长度设置，只需根据Y值设置，根据参考文献[13]可得长度应为

(19)

这里用符号step表示跳变步长，n表示处理的级数，那么在处理长度为Y的数据时，跳变步长与级数关系为

(20)

当处理长度为F的数据时，如果F与Y值相同那么可按照式(20)计算步进长度，如果F

(21)

通过如上介绍方法便可完成旋转因子的快速访问，在进行数据访问时每级处理需进行反序计算，这里可直接使用TS201处理其中的反序寻址指令. 获取数据后便利用内核的乘法器等资源完成蝶形运算. 为了充分发挥处理器的性能，程序中同时访问两个蝶形，分别放入x，y运算单元中，这样可在同一刻进行两个蝶形计算. 前面第2节介绍可知，合理安排流水后内核计算应隐藏在内存访问中，可把其作为优化目标，对蝶形计算进行软件流水建立，参考文献[14]中详细介绍了DSP排流水的方法，这里不再详细介绍. 排完指令流水后结果如图6所示.

如上便完成了每级蝶形流水线的建立，为了缩短处理时间，程序实现中将列处理的最后一级蝶形计算单独编写，以便可以同乘补偿因子操作融合在一起. 具体操作方法如前面所述.

程序优化后对其进行测试结果如表1所示.

通过表1可以看出采用本文介绍的映射方法实现基2FFT，在小点数时(小于2 048点)处理时间比正常映射方式要长，原因是小点数处理可存在cache中(TS201cache可缓存2 048点浮点型复数点数据)，在这种情况下采用本文介绍的映射方法会多出乘以补偿因子的时间. 当处理的数据长度大于2 048点以后，由于本文的映射方法充分利用cache，处理时间要明显优于正常的映射方法，速度约提高1倍. 另外通过该表可以看出理论计算时间同实际测试时间相近，可验证本文模型建立的正确性. 虽然表1所示正常映射方法中，大点数(2 048点)执行时间与理论时间存在较大偏差，其原因是DSP内部存储器采用的是DRAM，访问过程中在换页情况下会额外引入页激活时间，因此在换页和数据丢失同时存在时，额外开销要比手册给出的时间长，而本文理论计算时额外开销选取的只是手册提供数值，故大点数处理情况下会存在一定偏差.

表1 FFT执行时间对比表

4 结 论

针对雷达信号处理中常用的FFT算法到处理器上的映射进行研究. 首先分析了现代处理器的结构特点，基于分析结果将与处理性能相关的因素参量化，建立了基于现代化处理器的程序执行模型. 进而基于该模型对FFT算法在处理器上的执行过程分析，得知处理中的内存访问是影响实时性的关键因素. 进一步研究FFT的内访问特点，并建立模型进行分析，发现cache命中率直接影响算法执行的实时性. 对于采用正常映射的FFT方法，由于没有充分利用cache，导致大点数处理时效率低. 为了提高命中率，文中实现了一种可拆分的FFT映射方法，其将大点数划分成小点数处理，对cache进行充分利用，进而大大提高处理性能. 最后在TS201数字信号处理器上，以该映射方法为指导，实现了基2FFT算法，验证了本文所提供映射方法的有效性，同时验证了模型建立的正确性.

[1] Yang J， Xiong T， Peng Y. A fuzzy approach to filtering interferometric SAR data [J]. International Journal of Remote Sensing， 2007，28(6):1375-1382.

[2] Geer D. Chip makers turn to multicore processors[J]. Computer， 2005，38(5):11-13.

[3] Shen J P， Lipasti M H. Modern processor design: fundamentals of superscalar processors[M]. [S.l.]： Waveland Press， 2013.

[4] Suh G E， Rudolph L， Devadas S. Dynamic partitioning of shared cache memory[J]. The Journal of Supercomputing， 2004，28(1):7-26.

[5] Mazreah A A， Shalmani M T M. Low-leakage soft error tolerant dual-port SRAM cells for cache memory applications [J]. Microelectronics Journal， 2012，43(11):766-792.

[6] 范长永.32位RISC微处理器模块设计[D].北京：北京工业大学，2003.

Fan Changyong. 32 bit microprocessor module design[D]. Beijing: Beijing University of Technology， 2003. (in Chinese)

[7] Austin T， Larson E， Ernst D. Simple Scalar: an infrastructure for computer system modeling[J]. Computer， 2002，35(2):59-67.

[8] Benedetto S， Montorsi G. Unveiling turbo codes: some results on parallel concatenated coding schemes[J]. IEEE Transactions on Information Theory， 1996，42(2):409-428.

[9] 刘莉，高梅国，周闰,等.大点数FFT的多DSPs并行处理算法及实现[J].系统工程与电子技术，2003，25(10):1193-1196.

Liu Li， Gao Meiguo， Zhou Run, et al. Algorithm and implementation of a large-point FFT under the master-slave parallel multi-processor architecture[J]. Systems Engineering and Electronics, 2003，25(10):1193-1196. (in Chinese)

[10] 乔树山，黑勇，吴斌,等.块浮点FFT处理器的有限字长效应分析[J].电子科技大学学报，2008，37(1):58-60.

Qiao Shushan， Hei Yong, Wu Bin, et al. Finite word-length effects in implementation of a block floating-point FFT processor[J]. Journal of UEST of China, 2008，37(1):58-60. (in Chinese)

[11] Yi J J， Eeckhout L， Lilja D J， et al. The future of simulation: a field of dreams[J]. Computer， 2006，39(11):22-29.

[12] Magnusson P S， Christensson M， Eskilson J， et al. Simics: a full system simulation platform[J]. Computer， 2002，35(2):50-58.

[13] 王世一.数字信号处理[M].北京：北京理工大学出版社，1997.

Wang Shiyi. Digital signal processing[M]. Beijing: Publishing House of Electronics Industry， 1997. (in Chinese)

[14] 李方慧，王飞，何佩琨.TMS320C6000系列DSPs原理与应用[M]. 北京：电子工业出版社，2003.

Li Fanghui, Wang Fei, He Peikun. Principle and application of TMS320C6000 series DSPs[M]. Beijing: Publishing House of Electronics Industry, 2003. (in Chinese)

(责任编辑：刘芳)

An Efficient FFT-Mapping Method Based on Superscalar Processor

GAO Li-ning，ZHU Liang，LIU Teng-fei，LIU Feng

(School of Information and Electronics， Beijing Institute of Technology， Beijing 100081， China)

Considering the technology background and the structure characteristic of superscalar processor, novel mapping method was studied to implement the FFT effectively. Firstly, the structure of modern superscalar processor was modeled, and then the effect of FFT implementation on the processor was analyzed. The analysis results show that the EMS memory accessing is the key point of FFT algorithm implementation. Further more, the accessing process of FFT processing was modeled to implement FFT mapping effectively based on cache optimization. The new mapping method splits sequence of FFT into short sequence to exert the function of cache sufficiently. Finally, a radix-2 FFT was implemented on the ADI’s TS201 digital signal processor based on the new mapping method. The result shows that the implementation time of FFT is improved greatly.

fast Fourier transform (FFT); cache; superscalar processor

2015-04-02

国家自然科学基金资助项目(61370017)

高立宁(19XX—)，男，讲师，硕士生导师，E-mail:gaolining@bit.edu.cn.

TB 114.3

A

1001-0645(2016)09-0940-08

10.15918/j.tbit1001-0645.2016.09.011