面向GPU 平台的二维FFT 的加速技术研究

时间：2024-05-04

陈博伦，何卫锋

（上海交通大学微纳电子学系，上海200240）

0 引言

快速傅里叶变换（Fast Fourier Transform，FFT）算法是图像处理的基础算法之一，广泛应用于图像滤波、图像降噪、图像压缩等领域，同时也是基于频域分析的图像处理算法的基础。因此，针对FFT 运算的硬件加速对提升系统性能有重要的意义。

主要的硬件加速平台包括图形处理器（Graphics Processing Unit，GPU）、现场可编程逻辑门阵列（Field Programmable Gate Array，FPGA）等。其中，GPU 因其先进的并行架构体系，能够支持海量的多线程、高内存带宽以及大计算能力，在数千个核心上分配大规模多线程，并行执行独立任务而备受青睐。由于FFT 算法往往具有计算并行度高、数据流规整的特点，可以有效发挥GPU 大规模并行计算的潜力，因而受到了广泛的关注。

张全等人[1]分析了全局内存合并访问事务大小与GPU 占用率的关系，优化了二维FFT 算法在GPU 上的全局内存访问效率，使GPU 上二维FFT 算法的运行速度达到CUFFT 6.5 的1.3 倍。Govindaraju 等人[2]介绍了分层混合基数FFT 算法，有效利用了GPU 上的共享内存作为FFT 逐级运算时的缓存。另外，使用一些基于GPU 平台开发的商用算法库可以方便快捷地实现图像处理加速。例如NVIDIA 发布的CUDNN 算法库，能提供相关的函数接口，而且能自动适应不同的GPU 硬件平台，可移植性强[3]。然而，这些算法库都是闭源的，无法对其进行修改、分解和组合，灵活性差，不能适应所有的应用场景。

1 GPU上一维FFT算法的整体结构设计

1.1 算法原理

本文选取Cooley-Tukey 算法来实现FFT。Cooley-Tukey 算法是最常见的FFT 算法之一，其原理是用多个长度较小的离散傅里叶变换（Discrete Fourier Transform，DFT）来计算长度较大的DFT，经过反复迭代，最终分解成多个最小尺寸DFT。这个最小尺寸DFT 称作蝶形运算单元，而最小尺寸称作基数。在这基础上，根据基数的性质，又有高基和分裂基等多种算法。FFT分解基数的大小决定了分解后的级数。

图1 展示了以2 为基数的8 点Cooley-Tukey FFT算法，蝶形运算结构共3 级。最后一级的输出结果是顺序的，而第一级的输入序列的顺序需要重新排列，称为倒序排列。每个蝶形运算单元有自己的旋转因子，蝶形单元的输入需要先乘上对应的旋转因子，再进行蝶形运算。旋转因子的底数取决于FFT 点数，而指数与蝶形单元所在的级数及其本身的序数有关。

图1 基-2 Cooley-Tukey FFT算法

1.2 算法整体结构设计

基于一维FFT 算法的原理，我们首先分析一维FFT 算法的可并行性。以单一基数的一维FFT 为例，假设基数为R，长度为N=RM的FFT 算法运算结构共M级，每级包含N/R 个蝶形单元。同一级的蝶形单元互不相关，可并行执行。但后一级的蝶形单元的输入数据来源于前一级中R 个相关蝶形单元的输出数据，相邻两级存在依赖关系，不可并行，必须按顺序串行执行。

根据上述分析结果，本文设计了一维FFT 算法在GPU 上的并行实现过程，如图2 所示。图中的主要模块包括分解基数选取、倒序排列、旋转因子计算以及蝶形运算。分解基数的选取决定了每一级的蝶形运算结构。每一级的基数不必相同，基数不同时，蝶形单元的结构也不同。倒序排列发生在输入序列的读取阶段。每一级中的蝶形单元结构相同，但是指向不同的内存地址。两级蝶形结构之间通过共享内存进行数据通信。

图2 GPU上一维FFT算法并行实现过程

（1）旋转因子计算：使用常量内存查找表或CUDA内置函数库中的快速三角函数__sinf（）、__cosf（）来得到旋转因子θ的三角函数值。比较设备的运算能力和常量内存的带宽大小，从而决定到底使用哪种旋转因子θ的计算方式。计算密集型内核采取查找表方法，内存密集型则采取实时计算方法。

（2）倒序排列：倒序排列会导致全局内存的不连续访问。线程束中的线程访问全局内存不连续时，GPU会发出超过原本次数的访存事务，称为全局负载重播。全局负载重播会消耗过多的内存带宽，导致严重的性能问题。不产生负载重播的全局内存访问称为合并访问。设计中引入一级额外的共享内存作为缓存，在此缓存空间中进行倒序排列，之后倒序的序列传向后递到第一级蝶形运算结构。这种方法实际上将不连续的内存访问风险从全局内存转嫁到共享内存。幸而共享内存在不连续访问的情况下仍然可以实现较高的带宽利用率，从而解决全局内存访问不连续的问题。

（3）蝶形运算：蝶形运算包括两个部分：共享内存访问和运算。运算是简单的乘累加运算，而共享内存的访问则较为复杂。因为每一级蝶形结构的访存跨度在不断变化，所以共享内存访问过程中极容易发生bank conflict。以图1 为例，每一级蝶形单元的地址访问跨度Ns 依次为1、2、4，均小于32，因此线程有可能同时占用两个或以上的存储体，导致共享内存发生bank conflict。为了避免上述问题，本文设计了以蝶形单元访存跨度为依据的共享内存访问方式，写入数据时，每隔R×Ns 个数据填充一个无效数据。

（4）分解基数选取：FFT 算法分解过程中，基数越大，蝶形单元的级数越少，则同步操作次数越少，而且内存的访问次数降低。但是大基数也存在相应的劣势，基数越大，每个蝶形单元的规模更大，计算量更高，占用寄存器数量越大。由此可以推断，随着基数的增大，程序的性能瓶颈从内存密集型过渡到计算密集型。

比较GPU 设备的共享内存带宽和浮点运算吞吐量，可以定量分析FFT 的并行效率，从而选择合适的基数进行分解。另外，GPU 寄存器数量的限制意味着大基数分解的FFT 往往是不可取的，减少线程数可以减少寄存器占用的总数和使用的共享内存的数量，但是线程太少，线程束数量不足以隐藏访存延迟。

一个大小为R 的蝶形单元的浮点运算次数为5·r·R，其中r=log2R，每级共N/R 个蝶形单元，其理论运算时间如下所示，其中ComputeCapability 是GPU 的浮点运算能力。

假设数据缓存在共享内存，而且共享内存的访问效率为100%，那么共享内存访问的理论时间如下所示，其中DataWidth 是数据位宽，Bandwidth 是共享内存带宽。

图3 展示了共享内存访问和蝶形运算理论性能的比较。假设访存时间与运算时间相等时的基数等于Rx，此时内核同时达到共享内存带宽瓶颈以及计算单元运算速度瓶颈；当基数小于Rx 时，内核性能受限于共享内存有效带宽；当基数大于Rx 时，内核性能受限于GPU 运算能力。

图3 共享内存访问和蝶形运算的理论性能比较

2 GPU上二维FFT算法的并行设计

二维FFT 可以分解成行方向一维FFT 和列方向一维FFT。行变换和列变换不能同时进行，必须先完成其中之一，再进行另一变换。从整幅图像变换的可并行性来看，每一行或列的变换可以并行执行，即将图像的行或列变换分解成与图像行或列数相同的一系列独立的一维FFT 任务，然后将这些任务分配给线程块或线程。这些任务可以有限地合并在一起，因为图像任意两行或两列的一维FFT 的算法结构完全相同，可以共享旋转因子，从而减少旋转因子的计算次数。

一维FFT 的并行实现已在前文给出，在此基础上实现二维FFT 的关键是对全局内存的访问。全局内存中地址连续的数据被同一个线程束中的线程访问，而且满足内存对齐条件时，这些内存访问请求会被合并。合并的内存区间越大，全局内存的访问速度越快。GPU支持32byte、64byte、128byte 的合并内存模式。

行变换时的全局内存访问是连续的，如图4 所示。基数为16 的1024 点序列分割成16 个字段，每个字段长度为64。64 个线程依次顺序读取每一个字段到共享内存，然后进行一维FFT 运算。因此，行变换时的全局内存访问满足内存合并的要求，其带宽利用率理论上为100%。但是当基数增大时，每个字段的长度也相应地减短，导致线程数量减少，不利于程序实现充分的线程级并行。

图4 行变换时的全局内存访问

列变换时，全局内存访问不再连续，这会导致访存效率显著降低。为了提高访存效率，本文重新设计了列变换时的全局内存访问模式，使其满足内存合并条件，如图5 所示。线程块同时读取多列数据，当每个数据宽度为8Byte 时，连续4 个线程读取的数据合并为32Byte，满足最低限度的内存合并要求。

图5 列变换时的全局内存访问

3 二维FFT算法的性能测试结果及分析

实验使用的硬件平台为JETSON TX2，其搭载的Tegra X2 SoC 配置及性能如表1 所示。

表1 Tegra X2 配置及性能

首先进行二维FFT 的行变换和列变换的性能测试。二维图像的尺寸为1024×1024，每个数据的位宽为8Byte。其中行变换采取不同的数据访问模式，分别是32bit、64bit 以及128bit，研究数据访问模式对性能的影响。列变换采取不同的批量处理列数，研究批量处理的数据列数对性能的影响。列数越大，意味着同一个线程束中合并访问全局内存的线程越多。由于FFT 的运算量较低，理论运算时间与内核实际执行时间相比几乎可以忽略不计，因此这里使用全局内存的有效带宽来表征内核的性能。

基于上述测试方法得到的行变换性能如图6 所示。线程数量从32 增长到256 时，内核的性能有了10.5%至43.5%不等程度的提升。这是因为占用率逐渐增大，线程级并行度逐渐升高。随后线程数量继续增长，内核性能又降低了9.5%到11.4%。这是因为占用率不再增大，而线程的指令级并行度却逐渐降低。从不同数据访问模式的性能来看，64bit 和128bit 数据访问模式的性能相近，而二者的性能相比32bit 访问模式提升了15%左右，说明内核开发时应当尽可能选择64bit 和128bit 数据访问模式。如果数据位宽小于64bit，则应该将多个数据合并读取，从而提高全局内存访问速度。

图6 行变换性能测试结果

同理，我们进行了列变换的性能测试，结果如图7所示。随着线程束内合并访问全局内存的线程数量上升，内核性能有着300%左右的显著提升。当合并访问的全局内存宽度大于64Byte 后，性能趋于稳定。列变换中，全局内存的峰值有效带宽达到41GByte/s，与行变换时的带宽相同。因此本文的二维FFT 算法已经充分发挥了GPU 存储器的有效带宽，实现了峰值性能。

图7 列变换性能测试结果

基于行变换和列变换的性能测试及分析结果，我们接下来测试二维FFT 算法整体的性能，并将其与CUFFT 库函数的性能作比较，结果如图8 所示。随着图像尺寸的增大，FFT 算法的运算吞吐率逐渐升高，1024×1024 尺寸下的运算吞吐率相比128×128 尺寸提升了将近400%。当尺寸继续增大，运算吞吐率的增长幅度开始变缓，4096×4096 尺寸下的运算吞吐率相比1024×1024 尺寸只提升了42%。

图8 使用库函数与本文方法实现二维FFT的性能比较

随着二维FFT 的尺寸继续增大，运算吞吐率很难再继续提升，小范围内有可能反而会降低。因为随着尺寸增大，共享内存空间需求在不断增大，直至达到临界点。此时最大共享内存空间无法满足一次行变换所需，需要引入全局内存作为缓存。这意味着算法访问全局内存中的次数增多，会导致程序运行时间显著上升。