快速傅里叶变换（FFT）基本原理

引言

快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的加速算法，而DFT则是将连续的傅里叶变换离散化（在时域和频域离散），连续傅里叶变换可由傅里叶展开式推导得出。

正向逻辑为： 傅里叶展开式 \(\rightarrow\) 推导得到 \(\rightarrow\) 连续傅里叶变换对 \(\rightarrow\) 离散后得到 \(\rightarrow\) 离散傅里叶变换 \(\rightarrow\) 开发加速算法 \(\rightarrow\) 快速傅里叶变换 ## 傅里叶变换 ### 三角函数形式 对于一个周期为 \(T\) 的信号 \(x(t)\)，其基频为 \(f_1=1/T\)。对 \(x(t)\) 进行傅里叶展开：

\[ x(t) = \frac{a_0}{2} +\sum^{\infty}_ {n=1} \left[ a_n \cos\left(2\pi n f_1 t \right) + b_n \sin\left(2\pi n f_1 t \right) \right] \]

上式中的 \(a_n\) 和 \(b_n\) 是傅里叶系数：

\[ a_n = \frac{2}{T} \int^T_ 0 x(t) \cos \left( 2\pi n f_1 t \right) dt \]

\[ b_n = \frac{2}{T} \int^T_ 0 x(t) \sin \left( 2\pi n f_1 t \right) dt \]

指数形式

根据欧拉公式，可以将傅里叶变换式写为指数形式，指数形式看起来更为简洁，利于推导。 根据欧拉公式，有以下变换：

\[ \cos \left( 2\pi n f_1 t \right) = \frac{1}{2} \left[ e^{j2\pi n f_1 t} + e^{-j2\pi n f_1 t} \right] \]

\[ \sin \left( 2\pi n f_1 t \right) = \frac{1}{2j} \left[ e^{j2\pi n f_1 t} - e^{-j2\pi n f_1 t} \right] \]

欧拉公式：\(e^{ix} = \cos x + i \sin x\)

带入傅里叶展开式，可得：

\[ x(t) = \frac{a_0}{2} +\sum^{\infty}_ {n=1} \left[ \frac{a_n}{2} \left[ e^{ j 2\pi n f_1 t } + e^{ -j 2\pi n f_1 t } \right] + \frac{b_n}{2j} \left[ e^{ j 2\pi n f_1 t } - e^{ -j 2\pi n f_1 t } \right] \right] \]

合并同类项可得：

\[ x(t) = \frac{a_0}{2} +\sum^{\infty}_ {n=1} \left[ \frac{a_n - j b_n}{2} e^{ j 2\pi n f_1 t } + \frac{a_n + b_n}{2} e^{ -j 2\pi n f_1 t } \right] \]

观察指数项前的系数，令

\[ X(2\pi n f_1) = \frac{a_n - j b_n}{2} \]

则有：

\[ X(-2\pi n f_1) = \frac{a_n + j b_n}{2} \]

发现指数项前的系数具有奇偶性，利用奇偶性并改变求和符号的下限，将展开式改写为：

\[ x(t) = \sum^{\infty}_ {n=-\infty} X(2\pi n f_1) e^{ j 2\pi n f_1 t } \]

注意，此时求和符号的下限由一变为负无穷，\(1 \rightarrow -\infty\)。式中的 $X(2n f_1) $ 展开后为：

\[ \begin{gathered} X(2\pi n f_1) =\frac{1}{2} \left[ \frac{2}{T} \int^T_ 0 x(t) \cos \left( 2\pi n f_1 t \right) dt - \frac{2j}{T} \int^T_ 0 x(t) \sin \left( 2\pi n f_1 t \right) dt \right] \\\\ = \frac{1}{T} \int^T_ 0 x(t) \left[ \cos \left( 2\pi n f_1 t \right) - j \sin \left( 2\pi n f_1 t \right) \right] dt \end{gathered} \]

再利用欧拉公式，可得：

\[ X(2\pi n f_1) = \frac{1}{T} \int^T_ 0 x(t) e^{ -j 2\pi n f_1 t} dt \]

总结上述推导过程，得到指数形式的傅里叶变换对：

\[ x(t) = \sum^{\infty}_ {n=-\infty} X(2\pi n f_1) e^{ j 2\pi n f_1 t } \]

\[ X(2\pi n f_1) = \frac{1}{T} \int^T_ 0 x(t) e^{ -j 2\pi n f_1 t} dt \]

以上两个公式中，第二式才是傅里叶变换，第一式是傅里叶变换反演公式。

当信号没有周期，即周期无穷大时 \(T\rightarrow \infty\)，有\(f_1\rightarrow 0\)，此时第一式的求和符号可以改写为积分号，\(f_1\) 可以被当作是一个连续变化的量，傅里叶变换对就被改写为了连续形式。

傅里叶变换的实际意义

傅里叶变换对第一式：表示被采集信号 \(x(t)\) 可以由具有不同频率（\(2\pi n f_1\)）和幅值（\(X(2\pi n f_1)\)）的周期信号叠加得到。

傅里叶变换对第二式：为傅里叶变换，表示不同频率分量前的系数，可以理解为，该频率分量在被采集信号中的占比。

傅里叶变换的目的：得到被采集信号包含的频率成分和每个频率成分在原信号中的占比，用图形表示则为频谱图。

离散傅里叶变换

离散傅里叶变换对应的是在时域、频域都有限长，且都是离散的一类变换。

现有一个时域模拟信号 \(x(t)\)，对其进行采样，采样长度为 \(T\), 采样点有 \(N\) 个，采样时间间隔为 \(\Delta t\)。

则采样长度与采样点的关系为：

\[ T=N\Delta t \]

采样频率为：

\[ f_s = \frac{1}{\Delta t} \]

频谱的基频是\(f_1\)，即：

\[ f_1=\frac{1}{T} \]

根据傅里叶变换公式，由傅里叶变换求得的频谱谱线都是基频的整数倍，即频谱谱线的间隔：

\[ \Delta f=f_1=\frac{1}{T} = \frac{1}{N\Delta t} \]

记采样得到的离散的时域数字信号为\(x(k\Delta t)\)，把傅里叶变换式中的连续变量替换为离散变量

\[ t=k\Delta t,\quad 2\pi nf_1=2\pi n \Delta f \]

再将\(T=N\Delta t\)带入傅里叶变换对，对时间积分变为有限时间段内求和，即，\(\int \rightarrow \sum\)，\(dt\rightarrow\Delta t\)，得到：

\[ x(k\Delta t) = x(k)= \sum^{N-1}_ {n=0} X(2\pi n \Delta f) e^{ j 2\pi n k / N } \]

\[ X(2\pi n \Delta f) =X(n)= \frac{1}{N} \sum^{N-1}_ {k=0} x(k\Delta t) e^{ -j 2\pi n k /N} \]

习惯上将标定因子 \(N\) 移至反变换式中去，并且用 \(n\) 表示第 \(n\) 个采样点，用 \(k\) 表示第 \(k\) 条谱线（频率分量），即将上式中 \(n\) 和 \(k\) 的位置互换。总结上述结论有：

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

\[ X(k) = \sum^{N-1}_ {k=0} x(n) e^{ -j 2\pi n k /N} \]

上式就是离散傅里叶变换式。式中：

\[ k= 0,1,2\dots,N-1; \quad n= 0,1,2\dots,N-1 \]

根据离散傅里叶变换式可以求出第 \(n\) 条谱线对应的傅里叶系数的值，即该频率分量在整个信号中的占比。

我们注意到离散傅里叶变换式中求和符号的上下限改变了，不再是正负无穷。再改为正负无穷可以吗？改了之后两个公式含义还一样吗？答案是一样的，要用采样定理回答该问题。虽然标注是无穷，但是采样频率一定要大于2倍的被采样信号最高频率：

\[ f_s=\frac{1}{\Delta t}>2f_m \]

当采样频率已经确定为 \(1/ \Delta t\) 时，满足采样定理的 \(x(t)\) 的最高频率分量是：

\[ f_m<\frac{1}{2}f_s=\frac{1}{2\Delta t} = \frac{1}{2} N \Delta f \]

即满足采样定理的最后一根谱线所在的频率为\(1/(2N\Delta f)\)，换而言之 \(k\) 最大为 \(N/2\)。

快速傅里叶变换

FFT提出的背景

使用离散傅里叶变换时，计算一点的 \(X(n)\) 需要 \(N\) 次复数乘法和 \(N-1\) 次复数加法运算；计算全部的频谱需要 \(N^2\) 次复数乘法和 \(N(N-1)\) 次复数加法。而实现1次复数乘法需要4次实数乘法和2次实数加法，实现1次复数加法，需要2次实数加法。\(N=1024\)时，需要 \(1048576\) 次复数乘法运算。

有没有什么方法可以缩短计算时间呢？

1965年J. W. Cooley 和 J. W. Tukey 提出了一种快速求解DFT的算法，将乘法运算量由 \(N^2\) 降低至 \(\frac{N}{2} \log_ 2 N\) 次。以 \(N=1024\) 为例，他们提出的算法，只需要 \(5120\) 次复数乘法运算。

值得注意的是FFT不是特指某一种算法，而是指一类算法，自1965年后也有许多优秀学者提出了新的FFT算法。

FFT基本原理

对于 \(N\) 点序列 \(x(n)\)，其 DFT变换对为：

\[ x(n) = \frac{1}{N} \sum^{N-1}_ {n=0} X(k) W^{-nk}_ N \]

\[ X(k) = \sum^{N-1}_ {k=0} x(n) W^{nk}_ N \]

其中，\(W_ N = e^{-j2\pi /N}\) 。

DFT还可以写成矩阵形式，定义 \(\mathbf{W}_N\) 为：

\[ \mathbf{W}_N = \left[ W^{nk} \right]= \begin{bmatrix} W^0 & W^0 & W^0 & \cdots & W^0 \\\\ W^0 & W^1 & W^2 & \cdots & W^{N-1} \\\\ W^0 & W^2 & W^4 & \cdots & W^{2(N-1)} \\\\ \vdots & \vdots & \vdots & \vdots & \vdots \\\\ W^0 & W^{(N-1)} & W^{2(N-1)} & \cdots & W^{(N-1)(N-1)} \end{bmatrix} \]

\[ \bar{X}_N = \begin{bmatrix} X(0) & X(1) & \cdots & X(N-1) \end{bmatrix}^\mathrm{T} \]

\[ \bar{x}_N = \begin{bmatrix} x(0) & x(1) & \cdots & x(N-1) \end{bmatrix}^\mathrm{T} \]

则DFT的矩阵表达式为：

\[ \bar{X}_N = \mathbf{W}_N \bar{x}_N \]

观察矩阵形式的DFT表达式，可以发现DFT中包含大量重复运算，矩阵 \(\mathbf{W}_N\) 虽然有 \(N^2\) 个元素，但是其中只有 \(N\) 个独立值，并且一部分元素的值极为简单。\(W\) 因子的取值有如下特点：

• \(W^0=1\), \(W^{N/2}=-1\)
• \(W^{N+r}=W^r\), \(W^{N/2+r}=-W^r\)

对于4点的DFT，需要 \(4^2=16\) 次复数乘法，但是利用 \(W\) 因子的周期性和对称性可以大大简化计算，相关的DFT矩阵格式为：

\[ \begin{bmatrix} X(0) \\\\ X(1) \\\\ X(2) \\\\ X(3) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\\\ 1 & W^1 & -1 & -W^1 \\\\ 1 & -1 & 1 & -1 \\\\ 1 & -W^1 & -1 & W^1 \end{bmatrix}\begin{bmatrix} x(0) \\\\ x(1) \\\\ x(2) \\\\ x(3) \end{bmatrix} \]

将矩阵的第二列和第三列交换，可得：

\[ \begin{bmatrix} X(0) \\\\ X(1) \\\\ X(2) \\\\ X(3) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\\\ 1 & -1 & W^1 & -W^1 \\\\ 1 & 1 & -1 & -1 \\\\ 1 & -1 & -W^1 & W^1 \end{bmatrix}\begin{bmatrix} x(0) \\\\ x(2) \\\\ x(1) \\\\ x(3) \end{bmatrix} \]

由此可得：

\[ \begin{gathered} X(0) = [x(0) + x(2)] + [x(1) + x(3)] \\\\ X(1) = [x(0) - x(2)] + [x(1) - x(3)] W^1 \\\\ X(2) = [x(0) + x(2)] - [x(1) + x(3)] W^1 \\\\ X(3) = [x(0) - x(2)] - [x(1) - x(3)] \end{gathered} \]

利用上式计算DFT只需要 \(1\) 次复数乘法运算。处理长序列时，只需要将长序列分成类似于上式或比上式更简短的序列后，进行简单的运算，再按一定方式组合起来即可。

因此，FFT的基本原理是：先将原始的序列分解为一系列的短序列，充分利用 \(W\) 因子的周期性和对称性，进而求出这些短序列相应的DFT并进行适当组合，达到删除重复计算，减少乘法运算和简化结构的目的。

对于大部分FFT算法，有一些通用的概念和规律： * 级：将原信号每次折半，分为更小的单元，每折半一次，多出一级 * 蝶形单元：在FFT计算过程中，计算并不是顺序的，画出信号流图会发现，信号流图中包含着大量8字型（类似于蝴蝶）的计算单元 * 组：每一级的蝶形单元，按照其特性可以分为若干组 * \(W\)因子的分布：每一级的 \(W\) 因子具有一定的分布规律 * 码位倒置：使用FFT时，输出序列 \(\bar{X}_N\) 依照正序排列，但是输入序列的次序不再是自然序列，其排布次序和二进制码翻转、二进制与十进制转换有关

本文不再详细叙述具体的FFT算法。