库利-图基快速傅里叶变换算法

库利－图基快速傅里叶变换算法（英语：Cooley–Tukey FFT algorithm）是最常见的快速傅里叶变换算法。这一方法以分治法为策略递归地将长度为N = N₁N₂的DFT分解为长度分别为N₁和N₂的两个较短序列的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算法联合使用。

复杂度

离散傅立叶变换的复杂度为${\mathcal {O}}(n^{2})$（可参考大O符号）

快速傅立叶变换的复杂度为${\mathcal {O}}(N\log N)$，分析可见下方架构图部分，级数为$\log N$而每一级的复杂度为$N$，故复杂度为${\mathcal {O}}(N\log N)$

时域/频域抽取法

在FFT算法中，针对输入做不同方式的分组会造成输出顺序上的不同。如果我们使用时域抽取（Decimation-in-time），那么输入的顺序将会是位元反转排列（bit-reversed order），输出将会依序排列。但若我们采取的是频域抽取（Decimation-in-frequency），那么输出与输出顺序的情况将会完全相反，变为依序排列的输入与位元反转排列的输出。

时域抽取法

我们可以将DFT公式中的项目在时域上重新分组，这样就叫做时域抽取（Decimation-in-time），在这里$e^{-j(2\pi {\frac {nk}{N}})}$将会被代换为旋转因子（twiddle factor）$W_{N}^{nk}$。

$X[k]=\sum _{n=0}^{N-1}x[n]e^{-j(2\pi {\frac {nk}{N}})}\qquad k=0,1,\ldots ,N-1$

$=\sum _{n\ even}x[n]W_{N}^{nk}+\sum _{n\ odd}x[n]W_{N}^{nk}$

$=\sum _{r=0}^{(N/2)-1}x[2r]W_{N}^{2rk}+\sum _{r=0}^{(N/2)-1}x[2r+1]W_{N}^{(2r+1)k}$

$=\sum _{r=0}^{(N/2)-1}x[2r]W_{N}^{2rk}+W_{N}^{k}\sum _{r=0}^{(N/2)-1}x[2r+1]W_{N}^{2rk}$

$=\sum _{r=0}^{(N/2)-1}x[2r]W_{N/2}^{rk}+W_{N}^{k}\sum _{r=0}^{(N/2)-1}x[2r+1]W_{N/2}^{rk}$

$=G[k]+W_{N}^{k}H[k]$

在这边我们要注意的是，我们所替换的G[k]与H[k]具有周期性：

• ${\begin{cases}G[k+{\frac {N}{2}}]=G[k]\\H[k+{\frac {N}{2}}]=H[k]\end{cases}}$

还注意到系数具有对称性：

• $W_{N}^{k+N/2}=-W_{N}^{k}$

上述的推导可以划成下面的图：

划红框处所示的${\frac {N}{2}}$点DFT架构如下图所示：

划红框处所示的${\frac {N}{4}}$点DFT架构如下图所示：

下图是一个8点DIT FFT的完整架构图。

频域抽取法

我们可以将DFT公式中的项目在频域上重新分组，这样就叫做频域抽取（Decimation-in-frequency）。

首先先观察在频域上偶数项的部分：

$X[2r]=\sum _{n=0}^{N-1}x[n]W_{N}^{n(2r)}\ r=0,1,\cdots ,{\frac {N}{2}}-1$

$=\sum _{n=0}^{{\frac {N}{2}}-1}x[n]W_{N}^{2nr}+\sum _{n={\frac {N}{2}}}^{N-1}x[n]W_{N}^{2nr}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}x[n]W_{N}^{2nr}+\sum _{n=0}^{{\frac {N}{2}}-1}x[n+{\frac {N}{2}}]W_{N}^{2r[n+{\frac {N}{2}}]}$

${\color {Gray}\because W_{N}^{2r[n+{\frac {N}{2}}]}=W_{N}^{2rn}W_{N}^{rN}=W_{N}^{2rn}}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}x[n]W_{N}^{2rn}+\sum _{n=0}^{{\frac {N}{2}}-1}x[n+{\frac {N}{2}}]W_{N}^{2rn}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}(x[n]+x[n+{\frac {N}{2}}])W_{\frac {N}{2}}^{rn}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}g[n]W_{\frac {N}{2}}^{rn}$

再观察在频域上奇数项的部分：

$X[2r+1]=\sum _{n=0}^{N-1}x[n]W_{N}^{n(2r+1)}\ r=0,1,\cdots ,{\frac {N}{2}}-1$

$=\sum _{n=0}^{{\frac {N}{2}}-1}x[n]W_{N}^{n(2r+1)}+\sum _{n={\frac {N}{2}}}^{N-1}x[n]W_{N}^{n(2r+1)}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}x[n]W_{N}^{n(2r+1)}+\sum _{n=0}^{{\frac {N}{2}}-1}x[n+{\frac {N}{2}}]W_{N}^{(2r+1)[n+{\frac {N}{2}}]}$

${\color {Gray}\because W_{N}^{(2r+1)[n+{\frac {N}{2}}]}=W_{N}^{(2r+1)n}W_{N}^{(2r+1){\frac {N}{2}}}=-W_{N}^{(2r+1)n}}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}x[n]W_{N}^{(2r+1)n}-\sum _{n=0}^{{\frac {N}{2}}-1}x[n+{\frac {N}{2}}]W_{N}^{(2r+1)n}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}(x[n]-x[n+{\frac {N}{2}}])W_{N}^{n(2r+1)}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}(x[n]-x[n+{\frac {N}{2}}])W_{N}^{n}W_{\frac {N}{2}}^{nr}$

$=\sum _{n=0}^{{\frac {N}{2}}-1}(h[n]W_{N}^{n})W_{\frac {N}{2}}^{rn}$

上述的推导可以画成下面的图：

更进一步的拆解${\frac {N}{2}}$-point DFT的架构

下图为8点FFT下${\frac {N}{4}}$-point DFT的架构

总结上述架构，完整的8点DIF FFT架构图为

单一基底

利用库利－图基算法进行离散傅立叶拆解时，能够依需求而以2, 4, 8…等2的幂次方为单位进行拆解，不同的拆解方法会产生不同层数快速傅里叶变换的架构，基底越大则层数越少，复数乘法器也越少，但是每级的蝴蝶形架构则会越复杂，因此常见的架构为2基底、4基底与8基底这三种设计。

2基底

2基底-快速傅立叶算法（Radix-2 FFT）是最广为人知的一种库利－图基快速傅立叶算法分支。此方法不断地将N点的FFT拆解成两个'N/2'点的FFT，利用旋转因子$W^{nk},W^{\frac {N}{2}}$的对称性借此来降低DFT的计算复杂度，达到加速的功效。

而其实前述有关时域抽取或是频域抽取的方法介绍，即为2基底的快速傅立叶转换法。以下展示其他种2基底快速傅立叶算法的连线方法，此种不规则的连线方法可以让输出与输入都为循序排列，但是连线的复杂度却大大的增加。

4基底

4基底快速傅立叶变换算法则是承接2基底的概念，在此里用时域抽取的技巧，将原本的DFT公式拆解为四个一组的形式：

$X[k]=\sum _{n=0}^{N-1}x[n]e^{-j(2\pi {\frac {nk}{N}})}\qquad k=0,1,\ldots ,N-1$

$=\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+0]W_{N}^{(4n+0)k}+\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+1]W_{N}^{(4n+1)k}$

$+\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+2]W_{N}^{(4n+2)k}+\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+3]W_{N}^{(4n+3)k}$

$=W_{N}^{0}\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+0]W_{\frac {N}{4}}^{nk}+W_{N}^{1k}\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+1]W_{\frac {N}{4}}^{nk}$

$+W_{N}^{2k}\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+2]W_{\frac {N}{4}}^{nk}+W_{N}^{3k}\sum _{n=0}^{{\frac {N}{4}}-1}x[4n+3]W_{\frac {N}{4}}^{nk}$

$W_{N}^{0}F_{0}(k)+W_{N}^{k}F_{1}(k)+W_{N}^{2k}F_{2}(k)+W_{N}^{3k}F_{3}(k)$

在这里再利用$W^{nk+{\frac {N}{4}}}=-W^{nk+{\frac {3N}{4}}}=-jW^{nk}$的特性来进行与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=F_N‧x。F_N是一个N×N的DFT矩阵，其元素定义为[F_N]_ij=W_Nij(i,j ∈ [0,N-1])，当N=8时，我们可以得到以下的F₈矩阵并且进一步将其拆解。

在拆解成三个矩阵相乘之后，我们可以将复数运算的数量从56个降至24个复数的加法。底下是8基底的SFG。要注意的是所有的输出与输入都是复数的形式，而输出与输入的排序也并非规律排列，此种排列方式是为了要达到接线的优化。

混合基底

2/8基底

在2/8基底的算法中，我们可以看到我们对于偶数项的输出会使用2基底的分解法，对于奇数项的输出采用8基底的分解法。这个做法充分利用了2基底与4基底拥有最少乘法数与加法数的特性。当使用了2基底的分解法后，偶数项的输出如下所示。

$C_{2k}=\sum _{n=0}^{{\frac {N}{2}}-1}(x_{2n}+x_{{\frac {N}{2}}+n})W^{2nk}$

奇数项的输出则交由8基底分解来处理，如下四式所述。

$C_{8k+1}=\sum _{n=0}^{{\frac {N}{8}}-1}{\begin{Bmatrix}[(x_{n}-x_{n+{\frac {N}{2}}})-j(x_{n+{\frac {N}{4}}}-x_{n+{\frac {3N}{4}}})]+W^{\frac {N}{8}}[(x_{n+{\frac {N}{8}}}-x_{n+{\frac {5N}{8}}})-j(x_{n+{\frac {3N}{8}}}-x_{n+{\frac {7N}{8}}})]\end{Bmatrix}}W^{8nk}W^{n}$

$C_{8k+5}=\sum _{n=0}^{{\frac {N}{8}}-1}{\begin{Bmatrix}[(x_{n}-x_{n+{\frac {N}{2}}})-j(x_{n+{\frac {N}{4}}}-x_{n+{\frac {3N}{4}}})]-W^{\frac {N}{8}}[(x_{n+{\frac {N}{8}}}-x_{n+{\frac {5N}{8}}})-j(x_{n+{\frac {3N}{8}}}-x_{n+{\frac {7N}{8}}})]\end{Bmatrix}}W^{8nk}W^{5n}$

$C_{8k+3}=\sum _{n=0}^{{\frac {N}{8}}-1}{\begin{Bmatrix}[(x_{n}-x_{n+{\frac {N}{2}}})+j(x_{n+{\frac {N}{4}}}-x_{n+{\frac {3N}{4}}})]+W^{\frac {3N}{8}}[(x_{n+{\frac {N}{8}}}-x_{n+{\frac {5N}{8}}})+j(x_{n+{\frac {3N}{8}}}-x_{n+{\frac {7N}{8}}})]\end{Bmatrix}}W^{8nk}W^{3n}$

$C_{8k+7}=\sum _{n=0}^{{\frac {N}{8}}-1}{\begin{Bmatrix}[(x_{n}-x_{n+{\frac {N}{2}}})+j(x_{n+{\frac {N}{4}}}-x_{n+{\frac {3N}{4}}})]-W^{\frac {3N}{8}}[(x_{n+{\frac {N}{8}}}-x_{n+{\frac {5N}{8}}})+j(x_{n+{\frac {3N}{8}}}-x_{n+{\frac {7N}{8}}})]\end{Bmatrix}}W^{8nk}W^{7n}$

以上式子就是2/8基底的FFT快速算法。在架构图上可化为L型的蝴蝶运算架构，如下图所示。

而下图表示的是一个64点的FFT使用2/8基底的架构图。虽然2/8基底的算法缩减了$W^{\frac {N}{8}},W^{\frac {3N}{8}}$的乘法量，但是这种算法最大的缺点是跟其他固定基底或是混合基底比较起来，他的架构较为不规则。所以在硬件上比4基底或是2基底更难实现。

2/4/8基底

为了改进Radix 2/8在架构上的不规则性，我们在这里做了一些修改，如下图4.。此修改可让架构更加规则，且所使用的加法器与乘法器数量更加减少，如下表所示。

8ⁿ点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基底的运算，简化过的信号处理流程与蝴蝶型架构图可见下图

2²基底

基底的选择越大会造成蝴蝶形架构更加复杂，控制电路也会复杂化。因此Shousheng He和Mats Torkelson在1996提出了一个2^2基底的FFT算法，利用旋转因子的特性：$W_{\frac {N}{4}}=-j$。而–j的乘法基本上只需要将被乘数的实部虚部对调，然后将虚部加上负号即可，这样的负数乘法被定义为'简单乘法'，因此可以用很简单的硬件架构来实现。利用上面的特性，2²基底FFT能用串接的方式将旋转因子以4为单位对DFT公式进行拆解，将蝴蝶形架构层数降到log4N，大幅减少复数乘法器的用量，而同时却维持了2基底蝴蝶形架构的简单性，在性能上获得改进。2²基底DIF FFT算法的拆解方法如下列公式所述：

N点DFT公式： $X(k)=\sum _{n=0}^{N-1}x(n)W_{N}^{nk},0\leqq k<N$

利用线性映射将n与k映射到三个维度上面

${\begin{cases}n=<{\frac {N}{2}}n_{1}+{\frac {N}{4}}n_{2}+n_{3}>_{N}\\k=<k_{1}+2k_{2}+4k_{3}>_{N}\end{cases}}$

然后套用Common Factor Algorithm（CFA）

$X(k_{1}+2k_{2}+4k_{3})=\sum _{n_{3}=0}^{{\frac {N}{4}}-1}\sum _{n_{2}=0}^{1}\sum _{n_{1}=0}^{1}x({\frac {N}{2}}n_{1}+{\frac {N}{4}}n_{2}+n_{3})W_{N}^{({\frac {N}{2}}n_{1}+{\frac {N}{4}}n_{2}+n_{3})(k_{1}+2k_{2}+4k_{3})}$

$=\sum _{n_{3}=0}^{{\frac {N}{4}}-1}\sum _{n_{2}=0}^{1}{\begin{Bmatrix}B_{\frac {N}{2}}^{k_{1}}W_{N}^{({\frac {N}{4}}n_{2}+n_{3})k_{1}}\end{Bmatrix}}W_{N}^{({\frac {N}{4}}n_{2}+n_{3})(2k_{2}+4k_{3})}$

而蝴蝶型架构会变成以下形式

$B_{\frac {N}{2}}^{k_{1}}=x({\frac {N}{4}}n_{2}+n_{3})+(-1)^{k_{1}}x({\frac {N}{4}}n_{2}+n_{3}+{\frac {N}{2}})$

利用旋转因子$W_{\frac {N}{4}}=-j$的特性，可以观察出

$W_{N}^{({\frac {N}{4}}n_{2}+n_{3})(k_{1}+2k_{2}+4k_{3})}=W_{N}^{Nn_{2}n_{3}}W_{N}^{{\frac {N}{4}}n_{2}(k_{1}+2k_{2})}W_{N}^{n_{3}(k_{1}+2k_{2})}W_{N}^{4n_{3}k_{3}}$

$=(-j)^{n_{2}(k_{1}+2k_{2})}W_{N}^{n_{3}(k_{1}+2k_{2})}W_{N}^{4n_{3}k_{3}}$

再将此公式带入原式中可以得到

$X(k_{1}+2k_{2}+4k_{3})=\sum _{n_{3}=0}^{{\frac {N}{4}}-1}[H(k_{1},k_{2},n_{3})W_{N}^{n_{3}(k_{1}+2k_{2})}]W_{N}^{4n_{3}k_{3}}$

$H(k_{1},k_{2},n_{3})={\begin{matrix}\underbrace {{\begin{matrix}BF2I\\\overbrace {[x(n_{3})+(-1)^{k_{1}}(n_{3}+{\frac {N}{2}})]} \end{matrix}}{\begin{matrix}\\\\+(-j)^{k_{1}+2k_{2}}\end{matrix}}{\begin{matrix}BF2I\\\overbrace {[x(n_{3}+{\frac {N}{4}})+(-1)^{k_{1}}(n_{3}+{\frac {3N}{4}})]} \end{matrix}}} \\BF2II\end{matrix}}$

如上述公式所示，原本的DFT公式会被拆解成多个$H(k_{1},k_{2},n_{3})$，而$H(k_{1},k_{2},n_{3})$又可分为BF2I与BF2II两个层次结构，分别会对应到之后所介绍的两种硬件架构。

一个16点的DFT公式经过一次上面所述之拆解之后可得下面的FFT架构

可以看出上图的架构保留了2基底的简单架构，然而复数乘法却降到每两级才出现一次，也就是$log_{4}N$次。而BF2I以及BF2II所对应的硬件架构下图：

其中BF2II硬件单元中左下角的交叉电路就是用来处理-j的乘法。

一个256点的FFT架构可以由下面的硬件来实现：

其中图下方的为一7位元宽的计数器，而此架构的控制电路相当单纯，只要将计数器的各个位元分别接上BF2I与BF2II单元即可。

下表将2基底、4基底与2²基底算法做一比较，可以发现2²基底算法所需要的乘法气数量为2基底的一半，加法弃用量是4基底的一半，并维持一样的存储器用量和控制电路的简单性。

乘法器与加法器数量比较

	乘法数	加法数	存储器大小	控制电路
R2SDF	2(log4N-1)	4log4N	N-1	简单
R4SDF	log4N -1	8log4N	N-1	中等
R22SDF	log4N -1	4log4N	N-1	简单

2³基底

如上所述，2²算法是将旋转因子$W^{\frac {N}{4}}=-j$视为一个简单乘法，进而由公式以及硬件上的化简获得硬件需求上的改进。而借由相同的概念，Shousheng He和Mats Torkelson进一步将旋转因子$W^{\frac {N}{8}}={\frac {\sqrt {2}}{2}}(1-j)$的乘法化简成一个简单乘法，而化简的方法将会在下面讲解。

${\frac {\sqrt {2}}{2}}$乘法化简

在2基数FFT算法中的基本概念是利用旋转因子$W^{nk},W^{\frac {N}{2}}$的对称性，4基数算法则是利用 $W^{nk+{\frac {N}{4}}}=-W^{nk+{\frac {3N}{4}}}=-jW^{nk}$ 的特性。但是我们会发现在这些旋转因子的对称特性中─

$W^{\frac {N}{8}}=-W^{\frac {5N}{8}}={\frac {\sqrt {2}}{2}}(1-j),W^{\frac {3N}{8}}=-W^{\frac {7N}{8}}=-{\frac {\sqrt {2}}{2}}(1+j)$

─并没有被利用到。主要是因为${\frac {\sqrt {2}}{2}}(1-j)$的乘法运算会让整个FFT变得复杂，但是如果借由近似的方法，我们便可以将此一运算化简为12个加法。

${\frac {\sqrt {2}}{2}}=0.70710678=2^{-1}+2^{-3}+2^{-4}+2^{-6}+2^{-8}+2^{-9}$

我们可以从上式注意到，${\frac {\sqrt {2}}{2}}$可以被近似为五个加法的结果，所以${\frac {\sqrt {2}}{2}}(1+j)$就可以被简化为只有六个加法，再算入实部与虚部的计算，总共只需12个加法器就可以运用到此一简化特性。

经由与2²基底类似的推导，并用串接的方式将旋转因子以8为单位对DFT公式进行拆解，2³基底FFT算法进一步将复数乘法器的用量缩减到log₈N，并同时维持硬件架构的简单性。 　　 推导的方法与2²基底相当类似。借由这样的方法，2³基底能将乘法器的用量缩减到2基底的1/3，并同时维持一样的存储器用量以及控制电路的简单性。

一般性分解

除了常在应用中见到与$2$相关基底的拆解法，对于更加一般性的$N=N_{1}N_{2}$点离散傅立叶变换问题， 我们也有办法在理论上进行拆解，将问题化为数个$N_{1}$与$N_{2}$点离散傅立叶变换问题，并可对计算量进行估计。

而特别的是，透过互质在数论上的特性，对于$N_{1}$与$N_{2}$互质的情况，可以进一步节省一些运算， 在下面会特别分开讨论。

$N_{1}N_{2}$非互质

为了避免之后的符号混淆，我们先将$N_{1}N_{2}$置换为$P_{1}P_{2}$，也就是说接着要将$N=P_{1}P_{2}$点离散傅立叶变换， 想办法拆解为数个$P_{1}$与$P_{2}$点离散傅立叶变换问题。

接着定义要拆分的问题，要拆分的问题为$N$点离散傅立叶变换，将$f[n]$转换至$F[m]$：

$F[m]=\sum _{n=0}^{N-1}e^{-i{\frac {2\pi }{N}}mn}f[n]\qquad m,n=0,1,\ldots ,N-1$

直观地说，这个$N$点离散傅立叶变换，将由$n$作为参数的函数$f[n]$，转换成由$m$作为参数的函数$F[m]$， 并且$m,n$都有$N$个可能的数值。

待定义好要拆分的问题，便可以开始讨论如何进行拆分，基本概念是将有$N$个可能的数值的$m,n$， 分别化为个使用两个参数进行描述的函数$m=[m_{1},m_{2}],n=[n_{1},n_{2}]$，并借此将原问题化为二维度离散傅立叶变换， 便可使用数个较小的离散傅立叶变换问题描述整个过程。

而一种很直觉的转换方式，便是透过将$m,n$分别除以$P_{2},P_{1}$， 以商数与余数来做为参数描述$m,n$的值：

$n=n_{1}P_{1}+n_{2}\qquad m=m_{1}+m_{2}P_{2}$

$n_{1},m_{1}=0,1,\ldots ,P_{2}-1\qquad n_{2},m_{2}=0,1,\ldots ,P_{1}-1$

其中$n_{1}$作为将$n$除以$P_{1}$的商数，与作为$m$除以$P_{2}$的余数的$m_{1}$相同， 具有$P_{2}$个可能数值，同理$n_{2}$与$m_{2}$有$P_{1}$个可能数值。

将上述的参数代换及$N=P_{1}P_{2}$带入原式，可以得到：

$F[m_{1}+m_{2}P_{2}]=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}e^{-i{\frac {2\pi }{P_{1}P_{2}}}(m_{1}+m_{2}P_{2})(n_{1}P_{1}+n_{2})}f[n_{1}P_{1}+n_{2}]$

将右式的指数部分乘开并分项化简可以得到：

$F[m_{1}+m_{2}P_{2}]=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}f[n_{1}P_{1}+n_{2}]e^{-i{\frac {2\pi }{P_{2}}}m_{1}n_{1}}e^{-i{\frac {2\pi }{P_{1}}}m_{2}n_{2}}e^{-i{\frac {2\pi }{P_{1}P_{2}}}m_{1}n_{2}}e^{-i2\pi m_{2}n_{1}}$

最后透过$e^{-i2\pi }=1$与$P_{1}P_{2}=N$，可以得到：

$F[m_{1}+m_{2}P_{2}]=\sum _{n_{2}=0}^{P_{1}-1}\sum _{n_{1}=0}^{P_{2}-1}f[n_{1}P_{1}+n_{2}]e^{-i{\frac {2\pi }{P_{2}}}m_{1}n_{1}}e^{-i{\frac {2\pi }{N}}m_{1}n_{2}}e^{-i{\frac {2\pi }{P_{1}}}m_{2}n_{2}}$

观察上式，并加上括号辅助厘清运算顺序我们可以得到：

$F[m_{1}+m_{2}P_{2}]=\sum _{n_{2}=0}^{P_{1}-1}\left\{\left\{\sum _{n_{1}=0}^{P_{2}-1}f[n_{1}P_{1}+n_{2}]e^{-i{\frac {2\pi }{P_{2}}}m_{1}n_{1}}\right\}e^{-i{\frac {2\pi }{N}}m_{1}n_{2}}\right\}e^{-i{\frac {2\pi }{P_{1}}}m_{2}n_{2}}$

最内层的运算可以视为将双参数函数$f[n_{1},n_{2}]$中的一个参数$n_{1}$，透过离散傅立叶变换取代为由$m_{1}$描述， 得到一个新函数$G_{1}[m_{1},n_{2}]$(这步因为对每个不同$n_{2}$，都需要做一次将$n_{1}$取代为$m_{1}$的转换， 共需要$P_{1}$个$P_{2}$点离散傅立叶变换)：

$F[m_{1}+m_{2}P_{2}]=\sum _{n_{2}=0}^{P_{1}-1}\left\{G_{1}[m_{1},n_{2}]e^{-i{\frac {2\pi }{N}}m_{1}n_{2}}\right\}e^{-i{\frac {2\pi }{P_{1}}}m_{2}n_{2}}$

下一层的运算则可视为单纯的乘法，将$G_{1}[m_{1},n_{2}]$与$e^{-i{\frac {2\pi }{N}}m_{1}n_{2}}$相乘，得到$G_{2}[m_{1},n_{2}]$ (这步需要的计算量视${\frac {m_{1}n_{2}}{N}}$的特殊性而会有变化)：

$F[m_{1}+m_{2}P_{2}]=\sum _{n_{2}=0}^{P_{1}-1}G_{2}[m_{1},n_{2}]e^{-i{\frac {2\pi }{P_{1}}}m_{2}n_{2}}$

最后的运算可以视为将函数$G_{2}[m_{1},n_{2}]$中$n_{2}$，透过离散傅立叶变换取代为由$m_{2}$描述， 得到一个新函数$G_{3}[m_{1},m_{2}]$(这步因为对每个不同$m_{1}$，都需要做一次将$n_{2}$取代为$m_{2}$的转换， 共需要$P_{2}$个$P_{1}$点离散傅立叶变换)：

$F[m(=m_{1}+m_{2}P_{2})]=G_{3}[m_{1},m_{2}]$

就成功仅使用$P_{1}$与$P_{2}$点离散傅立叶变换，描述了原先的$N$点离散傅立叶变换。

而在这样的分解下，我们使用了$P_{1}$个$P_{2}$点离散傅立叶变换与$P_{2}$个$P_{1}$点离散傅立叶变换与一些额外的乘法， 并且这些额外使用的复数乘法$G_{1}[m_{1},n_{2}]\times e^{-i{\frac {2\pi }{N}}m_{1}n_{2}}$， 在电脑的运算架构中，如果${\frac {m_{1}n_{2}}{N}}$是${\frac {1}{4}}$的倍数则不需要使用乘法， 如果${\frac {m_{1}n_{2}}{N}}$是${\frac {1}{8}},{\frac {1}{12}}$的倍数则仅需两个实数乘法， 其他则需三个实数乘法，所以总运算量可以如下方式表示：

$P_{2}B_{1}+P_{1}B_{2}+3D_{3}+2D_{2}$

其中$B_{1}$是$P_{1}$傅立叶所需乘法数，$B_{2}$是$P_{2}$傅立叶所需乘法数， $D_{3}$是需三个实数乘法$m_{1}n_{2}$组合个数，$D_{2}$是需两个实数乘法$m_{1}n_{2}$组合个数。

而常见以$2$为基底的分解则是为了使离散傅立叶变换所需乘法数为零，这样就仅需考虑上面提到的额外乘法，可以提高效率也有较简单的结构。

$N_{1}N_{2}$互质

在$N_{1}N_{2}$互质的情况下，仍是采取和上面相近的思路来将问题进行拆分，首先，为了避免之后的符号混淆，我们同样将$N_{1}N_{2}$置换为$P_{1}P_{2}$。

接着同样定义要拆分的问题：

$F[m]=\sum _{n=0}^{N-1}e^{-i{\frac {2\pi }{N}}mn}f[n]\qquad m,n=0,1,\ldots ,N-1$

接着就是和上面的算法有最大差异的部分，在将$m,n$化为个使用两个参数进行描述的函数$m=[m_{1},m_{2}],n=[n_{1},n_{2}]$时， 最直觉的作法便是使用商数和余数，但在$P_{1},P_{2}$互质的情况下，可以有一些其他更具技巧性的选择。

当$P_{1},P_{2}$互质，对所有$n=0,1,\ldots ,N-1$我们可以找到唯一且不重复的一组$n_{1},n_{2}$使得：

$n=((n_{1}P_{1}+n_{2}P_{2}))_{N}=n_{1}P_{1}+n_{2}P_{2}+c_{1}N\qquad 0\leq n_{1}<P_{2}\qquad 0\leq n_{2}<P_{1}$

其中$((a))_{N}=a\mod N$，代表取余数的意思，$c_{1}$是一个整数。

例如假设$N=15,P_{1}=3,P_{2}=5$，则$n=1$对应到的$(n_{1},n_{2})$就是$(2,2)$， 有$2\times 3+2\times 5\mod 15=16\mod 15=1$。

并且对所有$n_{1},n_{2}$的组合(有$P_{1}\times P_{2}=N$组)，都对应到一个特定不重复的$n$。

同理我们可以把$m=0,1,\ldots ,N-1$表示为$m_{1},m_{2}$的双参数函数：

$m=((m_{1}P_{1}+m_{2}P_{2}))_{N}=m_{1}P_{1}+m_{2}P_{2}+c_{2}N\qquad 0\leq m_{1}<P_{2}\qquad 0\leq m_{2}<P_{1}$

将上述的参数代换及$N=P_{1}P_{2}$带入原式，可以得到：

$F[((m_{1}P_{1}+m_{2}P_{2}))_{N}]=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}e^{-i{\frac {2\pi }{N}}(m_{1}P_{1}+m_{2}P_{2}+c_{2}N)(n_{1}P_{1}+n_{2}P_{2}+c_{1}N)}f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]$

接着透过一次的展开化简及应用$e^{-i2\pi }=1$可得：

${\begin{aligned}F[((m_{1}P_{1}+m_{2}P_{2}))_{N}]&=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}e^{-i{\frac {2\pi }{N}}(m_{1}P_{1}+m_{2}P_{2})(n_{1}P_{1}+n_{2}P_{2})}e^{-i2\pi c_{2}(n_{1}P_{1}+n_{2}P_{2}+c_{1}N)}e^{-i2\pi c_{1}(m_{1}P_{1}+m_{2}P_{2}+c_{2}N)}e^{-i2\pi c_{1}c_{2}N}f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]\\&=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}e^{-i{\frac {2\pi }{N}}(m_{1}P_{1}+m_{2}P_{2})(n_{1}P_{1}+n_{2}P_{2})}f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]\end{aligned}}$

再将$N=P_{1}P_{2}$带入并再透过一次的展开化简及应用$e^{-i2\pi }=1$可得：

${\begin{aligned}F[((m_{1}P_{1}+m_{2}P_{2}))_{N}]&=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]e^{-i{\frac {2\pi }{P_{2}}}m_{1}P_{1}n_{1}}e^{-i{\frac {2\pi }{P_{1}}}m_{2}P_{2}n_{2}}e^{-i2\pi m_{1}n_{2}}e^{-i2\pi m_{2}n_{1}}\\&=\sum _{n_{1}=0}^{P_{2}-1}\sum _{n_{2}=0}^{P_{1}-1}f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]e^{-i{\frac {2\pi }{P_{2}}}m_{1}P_{1}n_{1}}e^{-i{\frac {2\pi }{P_{1}}}m_{2}P_{2}n_{2}}\end{aligned}}$

观察上式，并加上括号辅助厘清运算顺序我们可以得到：

$F[((m_{1}P_{1}+m_{2}P_{2}))_{N}]=\sum _{n_{2}=0}^{P_{1}-1}\left\{\sum _{n_{1}=0}^{P_{2}-1}f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]e^{-i{\frac {2\pi }{P_{2}}}m_{1}P_{1}n_{1}}\right\}e^{-i{\frac {2\pi }{P_{1}}}m_{2}P_{2}n_{2}}$

内层的运算可以视为将双参数函数$f[((n_{1}P_{1}+n_{2}P_{2}))_{N}]$中的一个参数$n_{1}$， 透过离散傅立叶变换取代为由一个与$m_{1}$有关的变量$m_{3}=((m_{1}P_{1}))_{P_{2}}$描述， 得到一个新函数$G_{1}[m_{3},n_{2}]$(这步因为对每个不同$n_{2}$，都需要做一次将$n_{1}$取代为$m_{3}$的转换， 共需要$P_{1}$个$P_{2}$点离散傅立叶变换)：

$F[((m_{1}P_{1}+m_{2}P_{2}))_{N}]=\sum _{n_{2}=0}^{P_{1}-1}G_{1}[m_{3},n_{2}]e^{-i{\frac {2\pi }{P_{1}}}m_{2}P_{2}n_{2}}$

外层的运算可以视为将函数$G_{1}[m_{3},n_{2}]$中的参数$n_{2}$， 透过离散傅立叶变换取代为由一个与$m_{2}$有关的变量$m_{4}=((m_{2}P_{2}))_{P_{1}}$描述， 得到一个新函数$G_{2}[m_{3},m_{4}]$(这步因为对每个不同$m_{3}$，都需要做一次将$n_{2}$取代为$m_{4}$的转换， 共需要$P_{2}$个$P_{1}$点离散傅立叶变换)：

$F[m]=F[((m_{1}P_{1}+m_{2}P_{2}))_{N}]=G_{2}[m_{3},m_{4}]=G_{2}[((m_{1}P_{1}))_{P_{2}},((m_{2}P_{2}))_{P_{1}}]$

最后透过$F$与$G_{2}$在不同$m_{1},m_{2}$时的点点数值对应关系， 就成功仅使用$P_{1}$与$P_{2}$点离散傅立叶变换，描述了原先的$N$点离散傅立叶变换。

而这个方法透过聪明的选取表达$m,n$的方式，使得拆解的过程中完全不需要多余的乘法运算， 总运算量可以简单表示为：

$P_{2}B_{1}+P_{1}B_{2}+3D_{3}+2D_{2}$

其中$B_{1}$是$P_{1}$傅立叶所需乘法数，$B_{2}$是$P_{2}$傅立叶所需乘法数。

虽然这个方法可以较上面的方法节省运算量， 但这个方法也牵涉较为复杂的$m,n$与$m_{1},m_{2},n_{1},n_{2}$转换，较为不直觉且不易理解， 也会遇到许多需要取余数的运算，可能会需要较大的空间建表进行查表法。

最后关于实际上要如何求得$n$与$n_{1},n_{2}$的转换关系， 可以先透过辗转相除法获取一对特定的$n_{11},n_{21}$使得：

$((n_{11}P_{1}+n_{21}P_{2}))_{N}=1$

然后我们可以知道对于任意整数$0\leq k<N$有：

$((kn_{11}P_{1}+kn_{21}P_{2}))_{N}=k=((kn_{11}P_{1}-c_{1}N+kn_{21}P_{2}-c_{2}N))_{N}=(((kn_{11}-c_{1}P_{2})P_{1}+(kn_{21}-c_{2}P_{1})P_{2}))_{N}$

然后就可以得到：

$n_{1k}=((kn_{11}))_{P_{2}}\qquad n_{2k}=((kn_{21}))_{P_{1}}$

满足：

$((n_{1k}P_{1}+n_{2k}P_{2}))_{N}=k\qquad 0\leq n_{1k}<P_{2}\qquad 0\leq n_{2k}<P_{1}$

参考资料

1. 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.
2. Lihong, J., G. Yonghong, et al. (1998). A new VLSI-oriented FFT algorithm and implementation. ASIC Conference 1998. Proceedings. Eleventh Annual IEEE International.
3. Duhamel, P. and H. Hollmann (1984). "Split-radix FFT algorithm." Electronics Letters 20(1): 14-16.
4. 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.
5. Richards, M. A. (1988). "On hardware implementation of the split-radix FFT." Acoustics, Speech and Signal Processing, IEEE Transactions on 36(10): 1575-1581.
6. Shousheng, H. and M. Torkelson (1996). A new approach to pipeline FFT processor. Parallel Processing Symposium, 1996., Proceedings of IPPS '96, The 10th International.
7. 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.