库利-图基快速傅里叶变换算法(英语:Cooley–Tukey FFT algorithm)[1]是最常见的快速傅里叶变换算法。这一方法以分治法为策略递归地将长度为N = N1N2的DFT分解为长度分别为N1和N2的两个较短序列的DFT,以及与旋转因子的复数乘法。这种方法以及FFT的基本思路在1965年詹姆斯·库利和约翰·图基合作发表《An algorithm for the machine calculation of complex Fourier series》之后开始为人所知。但后来发现,实际上这两位作者只是重新发明了高斯在1805年就已经提出的算法(此算法在历史上数次以各种形式被再次提出)。
库利-图基算法最有名的应用,是将序列长为N的DFT分割为两个长为N/2的子序列的DFT,因此这一应用只适用于序列长度为2的幂的DFT计算,即基2-FFT。实际上,如同高斯和库利与图基都指出的那样,库利-图基算法也可以用于序列长度N为任意因数分解形式的DFT,即混合基FFT,而且还可以应用于其他诸如分裂基FFT等变种。尽管库利-图基算法的基本思路是采用递归的方法进行计算,大多数传统的算法实现都将显示的递归算法改写为非递归的形式。另外,因为库利-图基算法是将DFT分解为较小长度的多个DFT,因此它可以同任一种其他的DFT算法联合使用。
复杂度
离散傅立叶变换的复杂度为 (可参考大O符号)
快速傅立叶变换的复杂度为 ,分析可见下方架构图部分,级数为 而每一级的复杂度为 ,故复杂度为
时域/频域抽取法
在FFT算法中,针对输入做不同方式的分组会造成输出顺序上的不同。如果我们使用时域抽取(Decimation-in-time),那么输入的顺序将会是比特反转排列(bit-reversed order),输出将会依序排列。但若我们采取的是频域抽取(Decimation-in-frequency),那么输出与输出顺序的情况将会完全相反,变为依序排列的输入与比特反转排列的输出。
时域抽取法
我们可以将DFT公式中的项目在时域上重新分组,这样就叫做时域抽取(Decimation-in-time),在这里 将会被代换为旋转因子(twiddle factor) 。
在这边我们要注意的是,我们所替换的G[k]与H[k]具有周期性:
-
还注意到系数具有对称性:
-
上述的推导可以划成下面的图:
划红框处所示的 点DFT架构如下图所示:
划红框处所示的 点DFT架构如下图所示:
下图是一个8点DIT FFT的完整架构图。
频域抽取法
我们可以将DFT公式中的项目在频域上重新分组,这样就叫做频域抽取(Decimation-in-frequency)。
首先先观察在频域上偶数项的部分:
再观察在频域上奇数项的部分:
上述的推导可以画成下面的图:
更进一步的拆解 -point DFT的架构
下图为8点FFT下 -point DFT的架构
总结上述架构,完整的8点DIF FFT架构图为
单一基底
利用库利-图基算法进行离散傅立叶拆解时,能够依需求而以2, 4, 8…等2的幂次方为单位进行拆解,不同的拆解方法会产生不同层数快速傅里叶变换的架构,基底越大则层数越少,复数乘法器也越少,但是每级的蝴蝶形架构则会越复杂,因此常见的架构为2基底、4基底与8基底这三种设计。
2基底
2基底-快速傅立叶算法(Radix-2 FFT)是最广为人知的一种库利-图基快速傅立叶算法分支。此方法不断地将N点的FFT拆解成两个'N/2'点的FFT,利用旋转因子 的对称性借此来降低DFT的计算复杂度,达到加速的功效。
而其实前述有关时域抽取或是频域抽取的方法介绍,即为2基底的快速傅立叶转换法。以下展示其他种2基底快速傅立叶算法的连线方法,此种不规则的连线方法可以让输出与输入都为循序排列,但是连线的复杂度却大大的增加。
4基底
4基底快速傅立叶变换算法则是承接2基底的概念,在此里用时域抽取的技巧,将原本的DFT公式拆解为四个一组的形式:
在这里再利用 的特性来进行与2基数FFT类似的化减方法,以降低算法复杂度。
8基底
在库利-图基算法里,使用的基底(radix)越大,复数的乘法与存储器的访问就越少,所带来的好处不言而喻。但是随之而来的就是实数的乘法与实数的加法也会增加,尤其计算单元的设计也就越复杂,对于可应用FFT之点数限制也就越严格。在表中我们可以看到在不同基底下所需的计算成本。
使用4096点FFT在不同基底下的计算量
动作 |
2基底 |
4基底 |
8基底
|
---|
复数乘法
|
22528 |
15360 |
10752
|
---|
实数乘法
|
0 |
0 |
8192
|
---|
复数加法
|
49152 |
49152 |
49152
|
---|
实数加法
|
0 |
0 |
8192
|
---|
存储器访问
|
49152 |
24576 |
16384
|
---|
在DFT的公式中,我们重新定义x=[x(0),x(1),…,x(N-1)]T, X=[X(0),X(1),…,X(N-1)]T,则DFT可重写为X=FN‧x。FN是一个N×N的DFT矩阵,其元素定义为[FN]ij=WNij(i,j ∈ [0,N-1]),当N=8时,我们可以得到以下的F8矩阵并且进一步将其拆解。
在拆解成三个矩阵相乘之后,我们可以将复数运算的数量从56个降至24个复数的加法。底下是8基底的的SFG。要注意的是所有的输出与输入都是复数的形式,而输出与输入的排序也并非规律排列,此种排列方式是为了要达到接线的优化。
混合基底
2/8基底
在2/8基底的算法中,我们可以看到我们对于偶数项的输出会使用2基底的分解法,对于奇数项的输出采用8基底的分解法。这个做法充分利用了2基底与4基底拥有最少乘法数与加法数的特性。当使用了2基底的分解法后,偶数项的输出如下所示。
奇数项的输出则交由8基底分解来处理,如下四式所述。
以上式子就是2/8基底的FFT快速算法。在架构图上可化为L型的蝴蝶运算架构,如下图所示。
而下图表示的是一个64点的FFT使用2/8基底的架构图。虽然2/8基底的算法缩减了 的乘法量,但是这种算法最大的缺点是跟其他固定基底或是混合基底比较起来,他的架构较为不规则。所以在硬件上比4基底或是2基底更难实现。
2/4/8基底
为了改进Radix 2/8在架构上的不规则性,我们在这里做了一些修改,如下图4.。此修改可让架构更加规则,且所使用的加法器与乘法器数量更加减少,如下表所示。
8n点FFT在不同算法下所需复数乘法量
N=8n |
2基底 |
4基底 |
2/4混合基底 |
2/4/8基底
|
---|
8
|
2 |
- |
2 |
0
|
---|
64
|
98 |
76 |
72 |
48
|
---|
512
|
1538 |
- |
1082 |
824
|
---|
4096
|
18434 |
13996 |
12774 |
10168
|
---|
在这里我们最小的运算单元称为PE(Process Element),PE内部包含了2/8基底、2/4基底、2基底的运算,简化过的信号处理流程与蝴蝶型架构图可见下图
22基底
基底的选择越大会造成蝴蝶形架构更加复杂,控制电路也会复杂化。因此Shousheng He和Mats Torkelson在1996提出了一个2^2基底的FFT算法,利用旋转因子的特性: 。而–j的乘法基本上只需要将被乘数的实部虚部对调,然后将虚部加上负号即可,这样的负数乘法被定义为'简单乘法',因此可以用很简单的硬件架构来实现。利用上面的特性,22基底FFT能用串接的方式将旋转因子以4为单位对DFT公式进行拆解,将蝴蝶形架构层数降到log4N,大幅减少复数乘法器的用量,而同时却维持了2基底蝴蝶形架构的简单性,在性能上获得改进。22基底DIF FFT算法的拆解方法如下列公式所述:
N点DFT公式:
利用线性映射将n与k映射到三个维度上面
然后套用Common Factor Algorithm(CFA)
而蝴蝶型架构会变成以下形式
利用旋转因子 的特性,可以观察出
再将此公式带入原式中可以得到
如上述公式所示,原本的DFT公式会被拆解成多个 ,而 又可分为BF2I与BF2II两个层次结构,分别会对应到之后所介绍的两种硬件架构。
一个16点的DFT公式经过一次上面所述之拆解之后可得下面的FFT架构
可以看出上图的架构保留了2基底的简单架构,然而复数乘法却降到每两级才出现一次,也就是 次。而BF2I以及BF2II所对应的硬件架构下图:
其中BF2II硬件单元中左下角的交叉电路就是用来处理-j的乘法。
一个256点的FFT架构可以由下面的硬件来实现:
其中图下方的为一7比特宽的计数器,而此架构的控制电路相当单纯,只要将计数器的各个比特分别接上BF2I与BF2II单元即可。
下表将2基底、4基底与22基底算法做一比较,可以发现22基底算法所需要的乘法气数量为2基底的一半,加法弃用量是4基底的一半,并维持一样的存储器用量和控制电路的简单性。
乘法器与加法器数量比较
|
乘法数 |
加法数 |
存储器大小 |
控制电路
|
---|
R2SDF
|
2(log4N-1) |
4log4N |
N-1 |
简单
|
---|
R4SDF
|
log4N -1 |
8log4N |
N-1 |
中等
|
---|
R22SDF
|
log4N -1 |
4log4N |
N-1 |
简单
|
---|
23基底
如上所述,22算法是将旋转因子 视为一个简单乘法,进而由公式以及硬件上的化简获得硬件需求上的改进。而借由相同的概念,Shousheng He和Mats Torkelson进一步将旋转因子 的乘法化简成一个简单乘法,而化简的方法将会在下面讲解。
乘法化简
在2基数FFT算法中的基本概念是利用旋转因子 的对称性,4基数算法则是利用
的特性。但是我们会发现在这些旋转因子的对称特性中─
─并没有被利用到。主要是因为 的乘法运算会让整个FFT变得复杂,但是如果借由近似的方法,我们便可以将此一运算化简为12个加法。
我们可以从上式注意到, 可以被近似为五个加法的结果,所以 就可以被简化为只有六个加法,再算入实部与虚部的计算,总共只需12个加法器就可以运用到此一简化特性。
经由与22基底类似的推导,并用串接的方式将旋转因子以8为单位对DFT公式进行拆解,23基底FFT算法进一步将复数乘法器的用量缩减到log8N,并同时维持硬件架构的简单性。
推导的方法与22基底相当类似。借由这样的方法,23基底能将乘法器的用量缩减到2基底的1/3,并同时维持一样的存储器用量以及控制电路的简单性。
除了常在应用中见到与 相关基底的拆解法,对于更加一般性的 点离散傅立叶变换问题,
我们也有办法在理论上进行拆解,将问题化为数个 与 点离散傅立叶变换问题,并可对计算量进行估计。
而特别的是,透过互质在数论上的特性,对于 与 互质的情况,可以进一步节省一些运算,
在下面会特别分开讨论。
非互质
为了避免之后的符号混淆,我们先将 置换为 ,也就是说接着要将 点离散傅立叶变换,
想办法拆解为数个 与 点离散傅立叶变换问题。
接着定义要拆分的问题,要拆分的问题为 点离散傅立叶变换,将 转换至 :
-
直观地说,这个 点离散傅立叶变换,将由 作为参数的函数 ,转换成由 作为参数的函数 ,
并且 都有 个可能的数值。
待定义好要拆分的问题,便可以开始讨论如何进行拆分,基本概念是将有 个可能的数值的 ,
分别化为个使用两个参数进行描述的函数 ,并借此将原问题化为二维度离散傅立叶变换,
便可使用数个较小的离散傅立叶变换问题描述整个过程。
而一种很直觉的转换方式,便是透过将 分别除以 ,
以商数与余数来做为参数描述 的值:
-
-
其中 作为将 除以 的商数,与作为 除以 的余数的 相同,
具有 个可能数值,同理 与 有 个可能数值。
将上述的参数代换及 带入原式,可以得到:
-
将右式的指数部分乘开并分项化简可以得到:
-
最后透过 与 ,可以得到:
-
观察上式,并加上括号辅助厘清运算顺序我们可以得到:
-
最内层的运算可以视为将双参数函数 中的一个参数 ,透过离散傅立叶变换取代为由 描述,
得到一个新函数 (这步因为对每个不同 ,都需要做一次将 取代为 的转换,
共需要 个 点离散傅立叶变换):
-
下一层的运算则可视为单纯的乘法,将 与 相乘,得到
(这步需要的计算量视 的特殊性而会有变化):
-
最后的运算可以视为将函数 中 ,透过离散傅立叶变换取代为由 描述,
得到一个新函数 (这步因为对每个不同 ,都需要做一次将 取代为 的转换,
共需要 个 点离散傅立叶变换):
-
就成功仅使用 与 点离散傅立叶变换,描述了原先的 点离散傅立叶变换。
而在这样的分解下,我们使用了 个 点离散傅立叶变换与 个 点离散傅立叶变换与一些额外的乘法,
并且这些额外使用的复数乘法 ,
在电脑的运算架构中,如果 是 的倍数则不需要使用乘法,
如果 是 的倍数则仅需两个实数乘法,
其他则需三个实数乘法,所以总运算量可以如下方式表示:
-
其中 是 傅立叶所需乘法数, 是 傅立叶所需乘法数,
是需三个实数乘法 组合个数, 是需两个实数乘法 组合个数。
而常见以 为基底的分解则是为了使离散傅立叶变换所需乘法数为零,这样就仅需考虑上面提到的额外乘法,可以提高效率也有较简单的结构。
互质
在 互质的情况下,仍是采取和上面相近的思路来将问题进行拆分,首先,为了避免之后的符号混淆,我们同样将 置换为 。
接着同样定义要拆分的问题:
-
接着就是和上面的算法有最大差异的部分,在将 化为个使用两个参数进行描述的函数 时,
最直觉的作法便是使用商数和余数,但在 互质的情况下,可以有一些其他更具技巧性的选择。
当 互质,对所有 我们可以找到唯一且不重复的一组 使得:
-
其中 ,代表取余数的意思, 是一个整数。
例如假设 ,则 对应到的 就是 ,
有 。
并且对所有 的组合(有 组),都对应到一个特定不重复的 。
同理我们可以把 表示为 的双参数函数:
-
将上述的参数代换及 带入原式,可以得到:
-
接着透过一次的展开化简及应用 可得:
-
再将 带入并再透过一次的展开化简及应用 可得:
-
观察上式,并加上括号辅助厘清运算顺序我们可以得到:
-
内层的运算可以视为将双参数函数 中的一个参数 ,
透过离散傅立叶变换取代为由一个与 有关的变量 描述,
得到一个新函数 (这步因为对每个不同 ,都需要做一次将 取代为 的转换,
共需要 个 点离散傅立叶变换):
-
外层的运算可以视为将函数 中的参数 ,
透过离散傅立叶变换取代为由一个与 有关的变量 描述,
得到一个新函数 (这步因为对每个不同 ,都需要做一次将 取代为 的转换,
共需要 个 点离散傅立叶变换):
-
最后透过 与 在不同 时的点点数值对应关系,
就成功仅使用 与 点离散傅立叶变换,描述了原先的 点离散傅立叶变换。
而这个方法透过聪明的选取表达 的方式,使得拆解的过程中完全不需要多余的乘法运算,
总运算量可以简单表示为:
-
其中 是 傅立叶所需乘法数, 是 傅立叶所需乘法数。
虽然这个方法可以较上面的方法节省运算量,
但这个方法也牵涉较为复杂的 与 转换,较为不直觉且不易理解,
也会遇到许多需要取余数的运算,可能会需要较大的空间建表进行查表法。
最后关于实际上要如何求得 与 的转换关系,
可以先透过辗转相除法获取一对特定的 使得:
-
然后我们可以知道对于任意整数 有:
-
然后就可以得到:
-
满足:
-
参考资料
- Widhe, T., J. Melander, et al. (1997). Design of efficient radix-8 butterfly PEs for VLSI. Circuits and Systems, 1997. ISCAS '97., Proceedings of 1997 IEEE International Symposium on.
- Lihong, J., G. Yonghong, et al. (1998). A new VLSI-oriented FFT algorithm and implementation. ASIC Conference 1998. Proceedings. Eleventh Annual IEEE International.
- Duhamel, P. and H. Hollmann (1984). "Split-radix FFT algorithm." Electronics Letters 20(1): 14-16.
- Vetterli, M. and P. Duhamel (1989). "Split-radix algorithms for length-pm DFT's." Acoustics, Speech and Signal Processing, IEEE Transactions on 37(1): 57-64.
- Richards, M. A. (1988). "On hardware implementation of the split-radix FFT." Acoustics, Speech and Signal Processing, IEEE Transactions on 36(10): 1575-1581.
- Shousheng, H. and M. Torkelson (1996). A new approach to pipeline FFT processor. Parallel Processing Symposium, 1996., Proceedings of IPPS '96, The 10th International.
- Shousheng, H. and M. Torkelson (1998). Designing pipeline FFT processor for OFDM (de)modulation. Signals, Systems, and Electronics, 1998. ISSSE 98. 1998 URSI International Symposium on.