快速傅里叶变换

时间：2020-05-06 23:37:15

一个五项余弦级数的时域信号。频率为

{displaystyle 10{sqrt {2}}}

缩放后的五项余弦级数的DFT。注意到显式积分更细的步长大小比FFT更精确地再现了峰值（4）和频率（56.569 Hz），代价是速度慢了上千倍。

快速傅里叶变换（英语：Fast Fourier Transform, FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法。傅里叶分析将信号从原始域（通常是时间或空间）转换到频域的表示或者逆过来转换。FFT会通过把DFT矩阵分解为稀疏（大多为零）因子之积来快速计算此类变换。因此，它能够将计算DFT的复杂度从只用DFT定义计算需要的 ${displaystyle O(n^{2})}$ ，其中 ${displaystyle n}$

直接按这个定义求值需要 O(N²) 次运算：X_k 共有 N 个输出，每个输出需要 N 项求和。直接使用DFT运算需使用N个复数乘法(4N 个实数乘法)与N-1个复数加法(4N-4个实数加法)，因此，计算使用DFT所有N点的值需要N²复数乘法与N²-N 个复数加法。FFT则是能够在 O(N log N) 次操作计算出相同结果的任何方法。更准确的说，所有已知的FFT算法都需要 O(N log N) 次运算（技术上O只标记上界），虽然还没有已知的证据证明更低的复杂度是不可能的。（Johnson and Frigo, 2007）

要说明FFT节省时间的方式，就得考虑复数相乘和相加的次数。直接计算DFT的值涉及到 N² 次复数相乘和 N(N−1) 次复数相加（可以通过削去琐碎运算（如乘以1）来节省 O(N) 次运算）。众所周知的基2库利-图基算法，N 为2的幂，可以只用 (N/2)log₂(N) 次复数乘法（再次忽略乘以1的简化）和 Nlog₂(N) 次加法就可以得到相同结果。在实际中，现代计算机通常的实际性能通常不受算术运算的速度和对复杂主体的分析主导（参见Frigo & Johnson, 2005），但是从 O(N²) 到 O(N log N) 的总体改进仍然能够体现出来。

一般的简化理论

假设一个M*N型矩阵S可分解成列向量以及行向量相乘：

${displaystyle mathbf {S} ={begin{bmatrix}a_{1}a_{2}vdots a_{m}end{bmatrix}}{begin{bmatrix}b_{1}&b_{2}&cdots &b_{n}end{bmatrix}}}$ 有 ${displaystyle M_{0}}$ where ${displaystyle mneq n}$ 有 ${displaystyle N_{0}}$ 个乘法。

${displaystyle {begin{bmatrix}Z[1]Z[2]vdots Z[N]end{bmatrix}}=mathbf {S} {begin{bmatrix}X[1]X[2]vdots X[N]end{bmatrix}}={begin{bmatrix}a_{1}a_{2}vdots a_{m}end{bmatrix}}{begin{bmatrix}b_{1}&b_{2}&cdots &b_{n}end{bmatrix}}{begin{bmatrix}X[1]X[2]vdots X[N]end{bmatrix}}}$

Step 2： ${displaystyle Z[1]=a_{1}Z_{a},Z[2]=a_{2}Z_{a},cdots ,Z[N]=a_{m}Z_{a}}$

${displaystyle S_{1}}$ 有 ${displaystyle P_{1}}$ 的乘法量上限为 ${displaystyle M_{0}+N_{0}+P_{1}}$ ，其中 ${displaystyle P_{1},P_{2},cdots ,P_{k}}$ 点DFT的乘法量为 ${displaystyle mathbf {B_{k}} }$ 点DFT的乘法量为：

{displaystyle {frac {N}{P_{1}}}B_{1}+{frac {N}{P_{2}}}B_{2}+cdots cdots +{frac {N}{P_{k}}}B_{k}}

，P是一个素数。

若 ${displaystyle mathbf {N_{1}=P^{a}} }$

且 ${displaystyle n_{1}times n_{2}}$ ）

有 ${displaystyle D_{1}}$ 及 ${displaystyle {frac {N}{8}}}$ 个值为 ${displaystyle {frac {N}{12}}}$ 的倍数，但不为 ${displaystyle {frac {N}{4}}}$

库利-图基算法

库利-图基算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为 ${displaystyle N=N_{1}N_{2}}$ 的 ${displaystyle N_{2}}$ 个旋转因子的复数乘法。

这种方法以及FFT的基本思路在1965年J. W. Cooley和J. W. Tukey合作发表An algorithm for the machine calculation of complex Fourier series之后开始为人所知。但后来发现，实际上这两位作者只是重新发明了高斯在1805年就已经提出的算法（此算法在历史上数次以各种形式被再次提出）。

库利-图基算法最有名的应用，是将序列长为N 的DFT分割为两个长为N/2 的子序列的DFT，因此这一应用只适用于序列长度为2的幂的DFT计算，即基2-FFT。实际上，如同高斯和Cooley与Tukey都指出的那样，Cooley-Tukey算法也可以用于序列长度N 为任意因数分解形式的DFT，即混合基FFT，而且还可以应用于其他诸如分裂基FFT等变种。尽管Cooley-Tukey算法的基本思路是采用递归的方法进行计算，大多数传统的算法实现都将显式的递归算法改写为非递归的形式。另外，因为Cooley-Tukey算法是将DFT分解为较小长度的多个DFT，因此它可以同任一种其他的DFT算法联合使用。

设计思想

下面，我们用N次单位根 ${displaystyle W_{N}}$ 。

${displaystyle W_{N}}$ 具有周期 ${displaystyle N}$

对称性：

{displaystyle W_{N}^{k+{frac {N}{2}}}=-W_{N}^{k}}

是

{displaystyle N}

为了简单起见，我们下面设待变换序列长度 ${displaystyle n=2^{r}}$ 时，可以将求和区间分为两部分：

$\begin{matrix} y_{k} = \sum_{n = 2 t} W_{N}^{k n} x_{n} + \sum_{n = 2 t + 1} W_{N}^{k n} x_{n} \\ = \sum_{t} W_{\frac{N}{2}}^{k t} x_{2 t} + W_{N}^{k} \sum_{t} W_{\frac{N}{2}}^{k t} x_{2 t + 1} \\ = F_{e v e n} (k) + W_{N}^{k} F_{o d d} (k) & (i \in Z) \end{matrix} {displaystyle {begin{matrix}y_{k}=sum _{n=2t}W_{N}^{kn}x_{n}+sum _{n=2t+1}W_{N}^{kn}x_{n}=sum _{t}W_{frac {N}{2}}^{kt}x_{2t}+W_{N}^{k}sum _{t}W_{frac {N}{2}}^{kt}x_{2t+1}=F_{even}(k)+W_{N}^{k}F_{odd}(k)&&&&&&(iin mathbb {Z} )end{matrix}}}$ 和 ${displaystyle F_{even}(k)}$ 奇数号和偶数号序列N/2点变换。由此式只能计算出 ${displaystyle y_{k}}$ 和 ${displaystyle F_{even}(k)}$

{displaystyle y_{k}=F_{even}(k)+W_{N}^{k}F_{odd}(k)}

算法实现

蝶形结网络和位反转（Bit Reversal）：

首先将

{displaystyle n=2^{N}}

逐层加入到蝶形网络中。

算法复杂度

由于按蝶形结网络计算n点变换要进行log n层计算，每层计算n个点的变换，故算法的时间复杂度为 ${displaystyle mathrm {O} (nlog n)}$ 且 ${displaystyle N_{1}}$ 互质的序列，可以采用基于剩余定理的互素因子算法将N长序列的DFT分解为两个子序列的DFT。与Cooley-Tukey算法不同的是，互素因子算法不需要旋转因子。

Rader-Brenner算法是类似于Cooley-Tukey算法，但是采用的旋转因子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radix variant of Cooley-Tukey所取代，与Rader-Brenner算法相比，有一样多的乘法量，却有较少的加法量，且不牺牲数值的准确性。

Bruun以及QFT算法是不断的把DFT分成许多较小的DFT运算。（Rader-Brenner以及QFT算法是为了2的指数所设计的算法，但依然可以适用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字）。尤其是Bruun算法，把FFT看成是 ${displaystyle z^{N}-1}$ 与 ${displaystyle z^{2M}+az^{M}+1}$ 分解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。这样只需要很少的乘法量（如果有需要的话），所以winograd是可以得到最少乘法量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd算法让DFT可以用 ${displaystyle 2^{k}}$ = ${displaystyle mathbf {X} _{k}^{*}}$ 点的DFT，也还没能够严谨的证明出FFT至少需要 ${displaystyle Omega (NlogN)}$ 。当 ${displaystyle N=2^{k}}$ 次无理实数的乘法即可以计算出来。更详尽，且也能趋近此下限的算法也一一被提出（Heideman & Burrus, 1986; Duhamel, 1990）。很可惜的是，这些算法，都需要很大量的加法计算，目前的硬件无法克服这个问题。

对于所需加法量的数目，虽然我们可以在某些受限制的假设下，推得其下限，但目前并没有一个精确的下限被推导出来。1973年，Morgenstern在乘法常量趋近巨大的情形下（对大部分的FFT算法为真，但不是全部）推导出加法量的下限： ${displaystyle Omega left(Nlog Nright)}$ ，但一般来说，此假设相当的不明确。长度为 ${displaystyle N=2^{k}}$ 是最少的。到目前为止，在长度为 ${displaystyle N=2^{k}}$ 还少。

还有一个问题是如何把乘法量与加法量的总和最小化，有时候称作"演算复杂度"（在这里考虑的是实际的运算量，而不是渐近复杂度）。同样的，没有一个严谨下限被证明出来。从1968年开始， ${displaystyle N=2^{k}}$ 的情形下，其需要 ${displaystyle 4Nlog _{2}N-6N+8}$ 。（Johnson and Frigo, 2007; Lundy and Van Buskirk, 2007）

大多数尝试要降低或者证明FFT复杂度下限的人都把焦点放在复数数据输入的情况，因其为最简单的情形。但是，复数数据输入的FFT算法，与实数数据输入的FFT算法，离散余弦变换（DCT），离散哈特列变换（DHT），以及其他的算法，均有很大的关连性。故任何一个算法，在复杂度上有任何的改善的话，其他的算法复杂度也会马上获得改善（Duhamel & Vetterli, 1990）。

快速算法查表:

当输入信号长度为N时:

N = 1~60

N	乘法数	加法数	N	乘法数	N	乘法数	N	乘法数
1	0	0	11	46	24	28	39	182
2	0	4	12	8	25	148	40	100
3	2	12	13	52	26	104	42	124
4	0	16	14	32	27	114	44	160
5	10	34	15	40	28	64	45	170
6	4	36	16	20	30	80	48	92
7	16	72	18	32	32	72	52	208
8	4	52	20	40	33	160	54	228
9	16	72	21	62	35	150	56	156
10	20	88	22	80	36	64	60	160

N < 1000

N	乘法数	N	乘法数	N	乘法数	N	乘法数
63	256	96	280	192	752	360	1540
64	204	104	468	204	976	420	2080
66	284	108	456	216	1020	480	2360
70	300	112	396	224	1016	504	2300
72	164	120	380	240	940	512	3180
80	260	128	560	252	1024	560	3100
81	480	144	436	256	1308	672	3496
84	248	160	680	288	1160	720	3620
88	412	168	580	312	1324	784	4412
90	340	180	680	336	1412	840	4580

N > 1000

N	乘法数	N	乘法数	N	乘法数	N	乘法数
1008	5356	1440	8680	2520	16540	4032	29488
1024	7436	1680	10420	2688	19108	4096	37516
1152	7088	2016	12728	2880	20060	4368	35828
1260	7640	2048	16836	3369	24200	4608	36812
1344	8252	2304	15868	3920	29900	5040	36860

参阅

离散傅里叶变换
并行快速傅里叶变换