快速傅里叶变换

用FFT计算DFT会得到与直接用DFT定义计算相同的结果；最重要的区别是FFT更快。（由于舍入误差的存在，许多FFT算法还会比直接运用定义求值精确很多，后面会讨论到这一点。）

令 x₀, ...., x_N-1 为复数。DFT由下式定义

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{i2\pi k{\frac {n}{N}}}}\qquad k=0,\dots ,N-1.}$ 

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

要说明FFT节省时间的方式，就得考虑复数相乘和相加的次数。直接计算DFT的值涉及到 N² 次复数相乘和 N(N−1) 次复数相加（可以通过削去琐碎运算（如乘以1）来节省 O(N) 次运算）。众所周知的基2库利-图基算法，N 为2的幂，可以只用 (N/2)log₂(N) 次复数乘法（再次忽略乘以1的简化）和 Nlog₂(N) 次加法就可以得到相同结果。在实际中，现代计算机通常的实际性能通常不受算术运算的速度和对复杂主体的分析主导^[7]，但是从 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 {\begin{bmatrix}a_{1}&a_{2}&\cdots &a_{m}\end{bmatrix}}^{T}}$ 有${\displaystyle M_{0}}$ 个相异的非平凡值（${\displaystyle a_{m}\neq \pm 2^{k},a_{m}\neq \pm 2^{k}a_{n}}$  where ${\displaystyle m\neq n}$ ）　

${\displaystyle {\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}}$ 有${\displaystyle N_{0}}$ 个相异的非平凡值

则S共需要${\displaystyle M_{0}*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 1：${\displaystyle Z_{a}=b_{1}X[1]+b_{2}X[2]+\cdots +b_{n}X[N]}$ 

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

简化理论的变型：

${\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}}+\mathbf {S} _{1}}$ 

${\displaystyle S_{1}}$ 也是一个M*N的矩阵。

若${\displaystyle S_{1}}$ 有${\displaystyle P_{1}}$ 个值不等于0，则${\displaystyle \mathbf {S} }$ 的乘法量上限为${\displaystyle M_{0}+N_{0}+P_{1}}$ 。

假设${\displaystyle N=P_{1}\times P_{2}\times \cdots \times P_{k}}$ ，其中${\displaystyle P_{1},P_{2},\cdots ,P_{k}}$ 彼此互素

${\displaystyle \mathbf {P_{k}} }$ 点DFT的乘法量为${\displaystyle \mathbf {B_{k}} }$ ，则${\displaystyle \mathbf {N} }$ 点DFT的乘法量为：

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

假设${\displaystyle \mathbf {N=P^{c}} }$ ，P是一个素数。

若${\displaystyle \mathbf {N_{1}=P^{a}} }$ 点的DFT需要的乘法量为${\displaystyle \mathbf {B_{1}} }$ 

且${\displaystyle n_{1}\times n_{2}}$ 当中（${\displaystyle n_{1}=0,1,\cdots ,N_{1}-1,\quad n_{2}=0,1,\cdots ,N_{2}-1}$ ）

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

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

则N点DFT的乘法量为：

${\displaystyle \mathbf {N_{2}B_{1}+N_{1}B_{2}+3D_{1}+2D_{2}} }$ 

库利-图基算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为${\displaystyle N=N_{1}N_{2}}$ 的离散傅里叶变换分解为长度为${\displaystyle N_{1}}$ 的${\displaystyle N_{2}}$ 个较短序列的离散傅里叶变换，以及与${\displaystyle \mathrm {O} (N)}$ 个旋转因子的复数乘法。

这种方法以及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 e^{-j{\frac {2\pi }{N}}}}$ 。

${\displaystyle W_{N}}$ 的性质：

1. 周期性，${\displaystyle W_{N}}$ 具有周期${\displaystyle N}$ ，即${\displaystyle W_{N}^{k+N}=W_{N}^{k}}$ 
2. 对称性：${\displaystyle W_{N}^{k+{\frac {N}{2}}}=-W_{N}^{k}}$ 。
3. 若${\displaystyle m}$ 是${\displaystyle N}$ 的约数，${\displaystyle W_{N}^{mkn}=W_{\frac {N}{m}}^{kn}}$ 

为了简单起见，我们下面设待变换序列长度${\displaystyle n=2^{r}}$ 。根据上面单位根的对称性，求级数${\displaystyle y_{k}=\sum _{n=0}^{N-1}W_{N}^{kn}x_{n}}$ 时，可以将求和区间分为两部分：

${\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)&&&&&&(i\in \mathbb {Z} )\end{matrix}}}$ 

${\displaystyle F_{odd}(k)}$ 和${\displaystyle F_{even}(k)}$ 是两个分别关于序列${\displaystyle \left\{x_{n}\right\}_{0}^{N-1}}$ 奇数号和偶数号序列N/2点变换。由此式只能计算出${\displaystyle y_{k}}$ 的前N/2个点，对于后N/2个点，注意${\displaystyle F_{odd}(k)}$ 和${\displaystyle F_{even}(k)}$ 都是周期为N/2的函数，由单位根的对称性，于是有以下变换公式：

• ${\displaystyle y_{k+{\frac {N}{2}}}=F_{even}(k)-W_{N}^{k}F_{odd}(k)}$ 
• ${\displaystyle y_{k}=F_{even}(k)+W_{N}^{k}F_{odd}(k)}$ 。

这样，一个N点变换就分解成了两个N/2点变换。照这样可继续分解下去。这就是库利-图基快速傅里叶变换算法的基本原理。根据主定理不难分析出此时算法的时间复杂度为${\displaystyle \mathrm {O} (N\log N)}$ 

算法实现

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

首先将${\displaystyle n=2^{N}}$ 个输入点列按二进制进行编号，然后对各个编号按位倒置并按此重新排序。例如，对于一个8点变换，

001倒置以后变成 100

000 → 000

001 → 100

010 → 010

011 → 110

100 → 001

101 → 101

110 → 011

111 → 111

倒置后的编号为{0,4,2,6,1,5,3,7}。

然后将这n个点列作为输入传送到蝶形结网络中，注意将因子${\displaystyle W_{N}^{k}}$ 逐层加入到蝶形网络中。

算法复杂度

由于按蝶形结网络计算n点变换要进行log n层计算，每层计算n个点的变换，故算法的时间复杂度为${\displaystyle \mathrm {O} (n\log n)}$ 。

在Cooley-Tukey算法之外还有其他DFT的快速算法。对于长度${\displaystyle N=N_{1}N_{2}}$ 且${\displaystyle N_{1}}$ 与${\displaystyle N_{2}}$ 互质的序列，可以采用基于中国剩余定理的互素因子算法将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^{M-1}}$ 与${\displaystyle z^{2M}+az^{M}+1}$ 的形式。

另一个从多项式观点的快速傅立叶变换法是Winograd算法。此算法把${\displaystyle z^{N}-1}$ 分解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。这样只需要很少的乘法量（如果有需要的话），所以winograd是可以得到最少乘法量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd算法让DFT可以用${\displaystyle 2^{k}}$ 点的DFT来简化，但减少乘法量的同时，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。

Rader算法提出了利用点数为N（N为素数）的DFT进行长度为N-1的回旋折积来表示原本的DFT，如此就可利用折积用一对基本的FFT来计算DFT。另一个prime-size的FFT算法为chirp-Z算法。此法也是将DFT用折积来表示，此法与Rader算法相比，能运用在更一般的转换上，其转换的基础为Z转换（Rabiner et al., 1969）。

在许多的运用当中，要进行DFT的资料是纯实数，如此一来经过DFT的结果会满足对称性：

${\displaystyle \mathbf {X} _{N-k}}$ =${\displaystyle \mathbf {X} _{k}^{*}}$ 

而有一些算法是专门为这种情形设计的（e.g. Sorensen, 1987）。另一些则是利用旧有的算法（e.g. Cooley-Tukey），删去一些不必要的演算步骤，如此省下了记忆体的使用，也省下了时间。另一方面，也可以把一个偶数长度且纯实数的DFT，用长度为原本一半的复数型态DFT来表示（实数项为原本纯实数资料的偶数项，虚数项则为奇数项）。

在Matlab上用一次N点FFT计算两个N点实数信号x,y的FFT:

function [X,Y] = Real2FFT( x,y )

% N-point FFT is used for 2 real signals both with length N

N = length(x);

z = x+y*i ;

Z = fft(z);

X = 0.5*(Z+conj(Z(1+mod(0:-1:1-N,N))));

Y = 0.5*(Z-conj(Z(1+mod(0:-1:1-N,N))))/i;

end

一度人们认为，用离散哈特利转换（Discrete Hartley Transform）来处理纯实数的DFT会更有效率，但接着人们发现，对于同样点数的纯实数DFT，经过设计的FFT，可以比DHT省下更多的运算。Bruun算法是第一个试着从减少实数输入的DFT运算量的算法，但此法并没有成为人们普遍使用的方法。

对于实数输入，且输入为偶对称或奇对称的情形，可以更进一步的省下时间以及记忆体，此时DFT可以用离散余弦转换或离散正弦转换来代替（Discrete cosine/sine transforms）。由于DCT/DST也可以设计出FFT的算法，故在此种情形下，此方法取代了对DFT设计的FFT算法。

DFT可以应用在频谱分析以及做折积的运算，而在此处，不同应用可以用不同的算法来取代，列表如下：

用来做频谱分析的情况下，DFT可用下列的算法代替：

• DCT
• DST
• DHT
• 正交基底的扩展（orthogonal basis expantion）包括正交多项式（orthogonal polynomials）以及CDMA.
• Walsh（Hadamard）转换
• Haar转换
• 小波（wavelet）转换
• 时频分布（time-frequency distribution）

用来做折积的情况下，DFT可用下列的算法代替：

• DCT
• DST
• DHT
• 直接做折积（direct computing）
• 分段式DFT折积（sectioned DFT convolution）
• 威诺格拉德快速傅里叶变换算法
• 沃尔什（Walsh、Hadamard）转换
• 数论转换

长久以来，人们对于求出快速傅里叶变换的复杂度下限以及需要多少的运算量感到很有兴趣，而实际上也还有许多问题需要解决。即使是用较简单的情形，即${\displaystyle 2^{k}}$ 点的DFT，也还没能够严谨的证明出FFT至少需要${\displaystyle \Omega (N\log N)}$ （不比${\displaystyle N\log N}$ 小）的运算量，目前也没有发现复杂度更低的算法。通常量学运算量的多寡会是运算效率好坏最主要的因素，但在现实中，有许多因素也会有很大的影响，如缓冲存储器以及CPU均有很大的影响。

在1978年，Winograd率先导出一个较严谨的FFT所需乘法量的下限：${\displaystyle \Theta (N)}$ 。当${\displaystyle N=2^{k}}$ 时，DFT只需要${\displaystyle 4N-2\log _{2}^{2}N-2\log _{2}N-4}$ 次无理实数的乘法即可以计算出来。更详尽，且也能趋近此下限的算法也一一被提出（Heideman & Burrus, 1986; Duhamel, 1990）。很可惜的是，这些算法，都需要很大量的加法计算，目前的硬件无法克服这个问题。

对于所需加法量的数目，虽然我们可以在某些受限制的假设下，推得其下限，但目前并没有一个精确的下限被推导出来。1973年，Morgenstern在乘法常量趋近巨大的情形下（对大部分的FFT算法为真，但不是全部）推导出加法量的下限：${\displaystyle \Omega \left(N\log N\right)}$ 。Pan（1986）在假设FFT算法的不同步的情形有其极限下证明出加法量的下限${\displaystyle \Omega (NlogN)}$ ，但一般来说，此假设相当的不明确。长度为${\displaystyle N=2^{k}}$ 的情形下，在某些假设下，Papadimitriou（1979）提出使用Cooley-Tukey算法所需的复数加法量${\displaystyle N\log _{2}N}$ 是最少的。到目前为止，在长度为${\displaystyle N=2^{k}}$ 情况，还没有任何FFT的算法可以让复数的加法量比${\displaystyle N\log _{2}N}$ 还少。

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

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

当输入信号长度为N时:

N = 1~60

N < 1000

N > 1000

• 离散傅里叶变换
• 并行快速傅里叶变换
• 快速数论变换