基于超大点数FFT优化算法的研究与实现

高立宁 马 潇刘腾飞 吴 金

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

1 引言

快速傅里叶变换(FFT)在雷达、数字通信和图像处理等科学技术领域有着广泛的应用，这使得FFT的工程实现具有十分重要的意义。随着现代数字信号处理技术的发展，应用系统中对实现超大点数FFT的需求越来越高，这也意味着对数字信号处理器(DSP)的处理性能和内存资源都提出了更高的要求和挑战[1]。

在雷达系统中，高分辨大测绘带宽的合成孔径雷达的飞速发展，对信号处理系统中实现超大点数FFT提出了更高的要求。然而，超大点数FFT的实现对资源的消耗随着点数的增加而增大，且对超大点数 FFT的优化效率往往和资源开销的程度直接相关，这都说明处理器资源已经成为制约大点数FFT实现的关键因素。因此，如何在节约现有平台资源开销并保证执行效率的条件下实现超大点数FFT，成为当前研究的重要方向。

很多文献对大点数 FFT的实现进行了深入的研究，对其性能的优化主要采用了两种方式：一种是通过算法的优化，使FFT处理更加适合处理器架构；另一种是增加系统资源，达到并行处理。例如，文献[2,3]中给出了 SingLeton结构实现定点/浮点FFT的方法，采用该结构对蝶形进行重排，使除了第1级外的其它级数据都是顺序读取，具有较高运算效率，但是这种方法要采用乒乓缓存，浪费了一倍的存储器；文献[4]采用多片DSP分级并行处理一个大点数FFT，将蝶形运算分配到多个DSP分级并行计算，进而提高FFT的执行效率，但是该处理方法需要多片DSP并行计算，造成了资源浪费，在实际工程应用中成本较高，且增加了开发难度；文献[5]提出了一种基于FPGA的高性能并行FFT处理方案，其采用4个蝶形运算单元进行并行处理，优化了处理效率，其缺点是增加了4倍的资源消耗；文献[6]采用Winograd算法将大点数FFT拆成小点数进行处理，虽然在处理中引入了额外的铰链因子相乘运算和3次显性转置操作，但是该算法能够有效地降低资源的消耗，能够在现有平台上实现更大点数的FFT运算。

本文针对处理器硬件资源对实现超大点数FFT的制约，提出一种优化的 Winograd算法实现超大点数FFT的方法。相比传统的Winograd处理方法，该实现方法通过优化铰链因子存储，改变矩阵访问方式和限制行列划分等手段使 FFT处理能够更好的适应处理器的架构，使其在现有平台上可以实现更大点数的FFT运算。

2 资源优化分析

文献[6]中Winograd算法实现FFT的主要流程为：(1)将1维序列拆分成2维矩阵并转置；(2)行方向FFT处理；(3)乘以铰链因子；(4)对(3)的结果转置存储；(5)列方向FFT处理；(6)将处理结果转置存储，得到序列的FFT结果。该处理算法将大点数FFT拆成小点数运算，但是该算法在实现过程中存在两个问题：一是引入的铰链因子需要额外的存储空间，同时增加了乘法运算次数；二是3次显性转置需要较多额外的内存空间，且存储转置处理引起的Cache丢失对时间影响很大。针对上述问题，本节给出具体的优化方案。

2.1 铰链因子存储优化

为了便于说明，这里重设2维序列N=ML×,M为列数，L为行数，并设1n和0n分别表示时域行列序号，0k和1k分别表示频域的行列序号。由文献[6]中给出的处理步骤可知，在第(3)步时对行方向FFT的处理结果乘以铰链因子，其值如式(1)：

由式(1)可推导：

由式(2)可以看出，铰链因子与0n和k两个变量有关，且可以将第k行的铰链因子分解为k个相乘，因此，在存储铰链因子时可以只存储第1行的铰链因子，而其它行的铰链因子可由第1行的铰链因子计算得出。虽然采用计算得出铰链因子会引入新的额外乘法运算，但该处理方式可以大大节省内存资源，这对在资源有限的情况下实现超大点数FFT十分重要。为了得出新的额外乘法运算对超大点数 FFT执行效率的影响，下面对该运算时间在FFT总运算时间中所占的比重作进一步分析[7]。

Winograd算法实现FFT运算所需运算量为

采用节省内存空间的方式存储铰链因子会额外引入N次乘法运算，又由于处理器进行乘法运算的时间比加法多得多，因此，如果只考虑FFT计算时间与乘法次数成正比，则新引入的额外乘法运算时间占FFT运算时间的比例为

由柯西不等式可以得

由式(6)和图1可见，新引入的额外乘法运算随着FFT点数的增加逐步减小，其对FFT的执行效率的影响也在逐渐减低。例如，按照理论分析在256k点时，额外乘法运算只占 FFT乘法总运算时间的18.49%；然而，该处理方法可以节约的内存空间，这为在现有处理器资源不变的情况下实现更大点数的FFT提供了可能。

经以上推导分析，为了节省处理器内存资源，本文采取只存储第1行的方法来存储铰链因子，且由此引入的新的额外乘法运算对 FFT执行时间的影响随点数的增加而逐渐减小；所以，对于超大点数FFT的实现来说，该部分的时间开销可以忽略不计。

2.2 矩阵访问优化

图1 随点数增加额外乘法运算时间/FFT乘法总运算时间曲线图

相对于 1维序列的处理，文献[6]的 Winograd算法将序列映射成2维序列处理，对2维矩阵的行列访问采用转置成正常顺序后连续读取的方式，对于FFT处理，该方法引入了3次显性转置，而转置处理需要大量的额外转置内存空间，且矩阵转置会对FFT的执行速度产生严重的影响。目前主流DSP处理器都采用分级存储结构，一级缓存(Cache)会对数据进行预处理来协助对矩阵数据的访问，从而减少直接对DRAM的访问次数；当数据在Cache中命中，可以直接从Cache中读取所需数据，否则需要重新到DRAM中读取新的数据页面，此时会引起数据访问时延[8]。文献[9]采用变换行列号方式实现对2维矩阵的访问，即通过行列号来计算数据首地址和访问步长，然后进行数据访问。但其带来的问题是，读取一列数据时，需要每读一个数据均进行行切换，而行切换需要切换时间；读取一行数据时，行数据不能太长，如果超出一级Cache的容量，在行FFT处理过程中，则会出现因无法从Cache命中数据带来的各种问题。因此，需要通过优化行列的拆分规则来最大发挥Cache协助数据访问的特点，实现对2维矩阵高效访问。

为了尽量减少行切换，可以在列处理时一次读取几列数据，读取列数的多少以及每列含数据量大小均需要根据Cache的容纳量计算，以实现对Cache的充分利用，但要保证读取列的总数据小于 Cache大小，否则在列FFT处理时会出现Cache频繁丢失的情况；图2所示为在进行一次列数据访问读取第i列的数据时，数据在Cache中的布局图。

在行FFT处理时为了充分利用Cache，只有当读一行的数据尽可能填充Cache空间时才能充分利用Cache的优势来协助数据访问。矩阵行处理访问内存后数据在Cache中分布如图3所示。

图2 矩阵列方向访问后的缓存布局

图3 矩阵行方向访问后的缓存布局

由此，采用行列号方式访问2维矩阵，节省了额外的转置存储空间，解决了3次显性转置引起时间开销对FFT执行效率影响的问题，并结合处理器分级存储的特点对 2维矩阵行列划分规则进行优化，最大限度发挥Cache的优点来提高矩阵访问速度。

上述两小节分别从铰链因子和矩阵转置两个方面对资源及效率进行优化：一是依据铰链因子的规律，采取减少存储铰链因子的数据量来节省内存空间，该方法引入了额外的乘法运算，但经过理论分析该额外运算引起的时间开销对超大点数 FFT的执行效率影响较小；二是改变矩阵的访问方式，节省了显性转置增加的额外存储空间，同时通过限制行列划分规则，在一定程度上提高了数据访问的速度。表1给出了本文方法采用Winograd算法的优化方法，实现超大点数FFT时对内存资源的需求。其中表示取其中最大的缓冲区，该缓冲区行列可复用。下面通过工程实例来具体说明采用Winograd算法的优化方法对超大点FFT的实现流程。

3 工程实例

TS201是ADI公司的TigerSHARC系列DSP，其采用超级哈佛结构，最高可支持600 MHz的内核时钟，具有双运算单元和双整数ALU，其内部集成24M bit的DRAM存储器，整个DRAM划分为6个分区(BANK)，每个 DRAM 存储器块带有一个128k bit的 Cache，用来协助数据访问。它的静态超标量结构使其能够在1个周期内执行4条指令，完成24次16 bit定点运算，12次32 bit定点运算6次浮点运算[10]。

表1 本文方法实现FFT对资源的需求

由于在一个BANK里面只能存储128k复数点的数据量，因此考虑将原始数据存储在两个不同的BANK中，采用J和K总线同时访问两个BANK的数据；另外，TS201的双运算单元可以一次读取4列8个复数点的数据，Cache的大小为128k bit，按照第2节所述的行列划分规则对1维矩阵进行划分：M=2048, L=64, i=8；即行处理时，每次读取一行的数据是 2048个复数点可以刚好填充 Cache的空间，在每次列处理时，分8次读取数据，一次读取4列，也刚好可以填充一个Cache的空间。因此，在行列FFT处理时，只在蝶形运算第1级和最后一级存在Cache丢失的情况下，其它级运算因为数据都缓存在Cache中而提高了数据命中概率。图4为128k复数浮点FFT的程序设计流程图。

本文主要从指令优化和汇编优化两个方面对128k复数浮点FFT的程序设计进行优化：在指令优化方面，TS201提供了合适复数运算的指令集，其内部有4条128 bit宽度的内部总线，通过J和K总线可以同时实现一个周期内4个32 bit数据的读和写操作，且其存在两个完全对称的X/Y计算单元和指令对齐缓冲器(IAB)为单个时钟周期内运行 4条指令提供了硬件基础，所以，理论上一个周期内可以完成两个蝶形运算；在汇编优化方面，TS201采用并行的超标量流水线技术，通过合理安排指令流水，使程序进入核循环时内核运算单元(乘法运算、加减法运算)的操作与内存访问的IO操作均并行起来，使程序只在填充和排空时存在一定的串行操作，保证了算法的高效运算[11]。

4 资源开销及速度对比

验证实验采用北京理工大学雷达所研制的8TS201通用信号处理板卡进行实测，即通过在FFT函数运行前后分别读取TS201 内部的时钟计数器，计算二者差值并除以TS201的处理器主频可以来得到运行时间，运行时间见表2。

由表 2可以看出，相对与文献[6]的 Winograd算法实现FFT，采用本文的方法对Winograd算法实现进行优化，其在小点数时因引入了额外的乘法运算致使执行效率大概降低16%；但是，在大点数时FFT的执行效率却得到了提升，这主要是因为该处理方法采用行列号方式优化矩阵访问使其省去了3次显性转置，并通过限制行列划分规则来减少Cache的丢失对执行效率的影响。另外，由表3可以看出，在超大点数时，本文 Winograd算法的实现方法比传统方法节省了近一半的内存资源，较其它算法也有巨大优势，这为同一处理器平台上实现更大点数的FFT提供了一种可行的方法。

5 结束语

本文针对现有处理器硬件平台对实现超大点数FFT的限制因素进行了深入分析，对Winograd算法的实现方法进行了资源优化处理，并依据处理器分级存储结构的特点优化了行列划分规则，从而提高了FFT处理过程中的行列访问效率；同时，给出了在TS201上实现了128k复数浮点FFT的实现流程。实验验证及资源对比表明，在超大点数FFT时，本文的处理方法能够节约一半的处理器资源开销，明显提升了大点数FFT处理的处理效率，该处理方法实现的超大点数FFT具有很高的工程应用价值。

图4 128k复数浮点FFT程序实现流程图

表2 FFT运行时间对比(s)μ

表3 FFT内存开销对比(kB)

[1] 李斌, 田素雷, 孙雪晶. 大点数 FFT 设计中提高资源利用率的方法[J]. 无线电工程, 2011, 41(1): 54-57.Li Bin, Tian Su-lei, and Sun Xue-jing. Method of improving resource utilization for large point FFT design[J]. Radio Engineering, 2011, 41(1): 54-57.

[2] Boris Lerner. Writing efficient floating-point FFTs for ADSP-TS201 TigerSHARC[OL]. http://www.Analog.com/dsp, 2012.

[3] 李欣, 刘峰, 龙腾. 定点FFT在TS201上的高效实现[J]. 北京理工大学学报, 2010, 30(1): 88-91.Li Xin, Liu Feng, and Long Teng. Efficient implementation of fixed-point FFT on TS201[J]. Transactions of Beijing Institute of Technology, 2010, 30(1): 88-91.

[4] 刘莉, 高梅国, 周闰, 等. 大点数FFT的多DSPs并行处理算法及实现[J]. 系统工程与电子技术, 2003, 25(10): 1193-1197.Liu Li, Gao Mei-guo, 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-1197.

[5] 石长振, 杨雪, 王贞松. 高性能并行 FFT 处理器的设计与实现[J]. 计算机工程, 2012, 38(2): 242-245.Shi Chang-zhen, Yang Xue, and Wang Zhen-song. Design and realization of high performance parallel FFT processor[J].Computer Engineering, 2012, 38(2): 242-245.

[6] Boris Lerner. Parallel implementation of fixed-point FFTs on TigerSHARC processors[OL]. http://www.Analog.com/dsp,2012.

[7] 王世一. 数字信号处理[M]. 北京: 北京理工大学出版社, 1997:123-131.

[8] 李浩, 谢伦国. 片上多处理器末级Cache优化技术研究[J]. 计算机研究与发展, 2012, 49(增刊): 172-179.Li Hao and Xie Lun-guo. Research development of optimization technology on last level cache in chip multiprocessors[J]. Journal of Computer Research and Development, 2012, 49(Suppl.): 172-179.

[9] 周永彬, 张军超. 基于软硬件的协同支持在众核上对 1-DFT算法的优化研究[J]. 计算机学报, 2008, 31(11): 2005-2014.Zhou Yong-bin and Zhang Jun-chao. Software ＆amp; hardware co-design for 1-D FFT optimization on many-core architecture[J]. Chinese Journal of Computers, 2008, 31(11):2005-2014.

[10] 刘书明, 罗勇江. ADSP TS20XS系列 DSP原理与应用设计[M]. 北京: 电子工业出版社, 2007: 6-7.

[11] Analog Device Inc. ADSP-TS201 TigerSHARC Processor Programming Reference[M]. Norwood: Mass, US, Analog Device Incorporation, 2004: 172-199.

[12] 刘志哲, 仲顺安. 基于分级存储并行运算的 FFT处理器设计[J]. 北京理工大学学报, 2011, 32(6): 691-684.Liu Zhi-zhe and Zhong Shun-an. Design of FFT processor based on grade-memory and parallel-computation[J].Transactions of Beijing Institute of Technology, 2011, 32(6):691-694.

[13] 苏涛, 庄德靖. 大点数 FFT 算法的改进及其实现[J]. 现代雷达, 2005, 27(7): 23-27.Su Tao and Zhuang De-jing. Improvement and implementation of FFT algorithm for long sequences[J].Modern Radar, 2005, 27(7): 23-27.

[14] Analog Device Inc. TigerSHARC DSP 32 bit REAL/COMPL EX FFT example [EB/OL]. http://www.Analog.com/dsp,2012.

[15] Analog Device Inc. TigerSHARC DSP complex fixed point FFT example for TS201 and TS101[EB/OL]. http://www.Analog.com/dsp, 2012.

[16] Bailey D H. FFTs in external or hierarchical memory[J].Journal of Supercomputing, 1990, 4(1): 23-35.