MindSpore AI科学计算系列 | “没有最快，只有更快”，快速傅里叶变换详解

背景

离散傅里叶变换（Discrete Fourier Transform，缩写为DFT）将信号的时域采样变换为频域的采样，是傅里叶变换在时域和频域上的离散形式。DFT是非常重要的离散变换，在许多实际应用中用于执行傅里叶分析。例如在数字信号处理中，声压、无线电信号或每日温度读数，都可以看作随时间变化的函数在有限时间间隔（窗口）内的采样。在图像处理中，样本可以是沿光栅图像的行或列的像素值。在数学领域中，DFT还用于有效求解偏微分方程、卷积或大整数相乘等运算。在AI4Science领域，DFT也是很多算法(如Fourier Neural Operator，FNO)的基础。

图1. 傅里叶变换示意图

对于给定长度N的复数序列，DFT如式1所示，其中，x为N个复数的输入序列，X为对应的输出序列，指数幂项称为常被称为“旋转因子”(twiddle factor)。

式1

式2

从式2可以看出，DFT本质上是向量矩阵操作，计算N个点的DFT需要N^2次复数乘法，N(N-1)次复数加法，总的时间复杂度为O(N^2×mul+N(N-1)×add)。

快速傅里叶变换（Fast Fourier Transform, 缩写为FFT）是高效、快速计算离散傅里叶变换（DFT）或其逆变换的统称。FFT利用DFT的特性，例如对称性和周期性，去减少DFT中的重复计算，降低复杂度。FFT算法的提出，使DFT在工程、科学和数学领域等领域得到了更加广泛的应用，被IEEE科学与工程计算期刊评选为20 世纪十大算法之一。

1、方法

1.1 库利-图基算法

库利-图基算法（Cooley-Tukey Algorithm）是目前应用最广泛、最流行的FFT算法。这一方法以分治法为策略递归地将长度为N的DFT分解为长度N1和N2的子序列的DFT，以及与旋转因子的复数乘法。此处以基为2的Cooley-Tukey算法为例，即N的长度为2的幂次(不满足时补零)，将序列长为N的DFT分割为两个长为N/2的子序列的DFT。

Cooley-Tukey算法利用了DFT的对称性和周期性，减少DFT中的重复计算，降低复杂度。推导如下：

可以看出，蝶形网络中的每次蝶形运算需要一次复数乘法以及两次复数加法。对于n=8的采样序列x(n)，可推出以下蝶式网络

图3. 8点蝶形网络

按照Cooley-Tukey算法流程，在以2为基底的条件下，假设给定长度为N=2^L长度的复数序列，共需要2^(L-1)次复数乘法和2^L次复数加法，依次递推后，总的计算复杂度为O(N/2*logN×mul+NlogN×add)，相较于DFT，显著的降低了计算复杂度。

进一步的，Cooley-Tukey算法还可以选择不同的基底为单位对序列进行拆解，不同的拆解方法会产生不同层数快速傅里叶变换的架构，基底越大则层数越少，复数乘法器也越少，但是每级的蝴蝶形架构则会越复杂，因此常见的架构为2基底、4基底与8基底这三种设计。

1.2 分裂基算法(Split-radix)

分裂基算法（Split-radix FFT Algorithm, SRA）由Duhamel和Hollman于1984提出。SRA算法是目前众多 FFT 算法中乘法和加法次数最少的算法，而且它具有良好的运算结构，以及较短的运算程序。SRA算法基本思路就是对偶序列使用radix-2FFT算法，对奇序列使用radix-4FFT 算法。

图4. 分裂基算法示意图

1.3 互质因子算法(Prime Factor Algorithm, PFA)

互质因子算法（Prime-factor FFT algorithm, PFA），又称为Good-Thomas算法，也是一种快速傅里叶变换。它是把输入规模为N的DFT，分解为N_1 × N_2大小的二维DFT，其中N_1与N_2互质。变成大小为N1和N2的DFT之后，可以继续递归使用 PFA，或选择其他 FFT 算法(例如 Cooley-Tukey)来计算。

互质因子算法将一维的 DFT 问题转变为多维问题来进行计算。总的时间复杂度：O(N(N_1+N_2)(mul+add))。

1.4 Rader's FFT算法

Rader算法是针对输入大小为质数DFT的快速算法，该算法将长度为N的DFT重新表示为长度为N-1的循环卷积

式5

循环卷积又是可分解的，再根据卷积定理，将N-1的循环卷积再分解为更小规模的DFT问题。

2、结论

通过学习上述的几种FFT方法，可总结出以下几点：

2.1 蝶形网络的组织方式将影响整个优化

蝶形网络决定了数据访问模式和蝶形计算执行顺序，相同的蝶形网络不同的实现和优化可能导致不同的性能。例如相对于Cooley-Tukey算法来说，PFA算法虽然在浮点计算量上有所减少，但是因为其复杂的映射关系，这将导致更多的数据存取，因此PFA算法适用于数据存取开销小于浮点计算开销的平台架构。

2.2 蝶形计算的性能将直接影响FFT算法的最终性能

在FFT计算过程中蝶形计算被反复调用，例如，根据DFT本身的周期性与对称性，可以挖掘蝶形计算过程中的共同项，供各级重复利用，这将减少大量的浮点计算开销。

2.3 实现高性能FFT需要与硬件架构相结合

例如，蝶形计算时互相独立的，采用SIMD可并行处理多个蝶形计算，并行的数量又依赖其余部件的情况，如浮点寄存器的个数。所以，硬件本身能否充分利用，直接关系到FFT程序性能的提升。

3、总结

随着异构计算框架的发展，使用不同类型指令集和体系架构的计算单元组成系统的计算方式已经非常普遍。部分研究人员对于FFT算法的研究，从降低FFT本身的计算开销转移至在特定的硬件架构下实现高性能的FFT。在这方面，已经有一些成熟的高性能快速傅里叶变换算法库，例如FFTW、ARMPL、Intel MKL、鲲鹏的KML_FFT、AOCL（AMD optimizing CPU libraries）、CUFFT（CUDA fast Fourier transform library）等。

昇腾异构计算架构也在结合自身优势，从提升计算效率、降低功耗等方面，实现高性能的FFT库，预计不久将与大家见面。

参考文献

[1]Cooley J W, Tukey J W. An algorithm for the machine calculation of complex Fourier series[J]. Mathematics of computation, 1965, 19(90): 297-301.

[2]Good I J. The interaction algorithm and practical Fourier analysis[J]. Journal of the Royal Statistical Society Series B: Statistical Methodology, 1958, 20(2): 361-372. https://doi.org/10.1111/j.2517-6161.1958.tb00300.x

[3]Thomas L H. Using a computer to solve problems in physics[J]. Applications of digital computers, 1963: 44-45.

[4]Duhamel P, Hollmann H. ‘Split radix’FFT algorithm[J]. Electronics letters, 1984, 20(1): 14-16.

[5]Rader C M. Discrete Fourier transforms when the number of data samples is prime[J]. Proceedings of the IEEE, 1968, 56(6): 1107-1108. DOI: 10.1109/PROC.1968.6477

[6]陈暾, 李志豪, 贾海鹏, 等. 基于 ARMv8 平台的多维 FFT 实现与优化研究[J]. 计算机学报, 2019, 42(11): 2384-2402.

[7]郭金鑫, 张广婷, 张云泉, 等. Cooley-Tukey FFT 算法高性能实现与优化研究[J]. 计算机科学与探索, 2022, 16(6): 1304. DOI: 10.3778/j.issn.1673-9418.2011092

往期回顾

MindSpore AI科学计算系列 | 有效改善预报模糊问题，生成式模型在短临降水预报中大放异彩

MindSpore AI科学计算系列 | 基于深度学习模型来替代传统DFT模型以及DeephE3模型的分析综述

MindSpore AI科学计算系列 | 以MindSpore Elec为例的智能电磁计算若干进展综述

MindSpore AI科学计算系列 | 最新综述文章梳理量子到宏观尺度AI4S共性介绍