【FFT算法深度解析】：从原理到应用，一文读懂快速傅里叶变换

# 1. FFT算法基础理论

快速傅里叶变换（FFT）是一种算法，用于计算离散傅里叶变换（DFT）。DFT将时域信号转换为频域表示，从而揭示信号的频率分量。

FFT算法利用了DFT的周期性和对称性，通过分治策略将DFT的计算复杂度从O(N^2)降低到O(N log N)，其中N是信号的长度。FFT算法的核心思想是将信号分解为较小的子序列，对每个子序列进行DFT，然后将结果合并以得到原始信号的DFT。

傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。

# 2.1 FFT算法的数学原理

### 傅里叶变换

FFT算法是基于傅里叶变换的一种快速算法。傅里叶变换是一种数学变换，它将一个时域信号（如时间序列）转换为一个频域信号（如频率谱）。

时域信号表示信号在时间轴上的变化，而频域信号表示信号中不同频率分量的幅度和相位。傅里叶变换将时域信号分解为一系列正弦波和余弦波，每个正弦波和余弦波对应于信号中的一个特定频率分量。

根据原信号的不同类型，傅里叶变换可以分为四种类别：

(1)非周期性连续信号傅里叶变换

(2)周期性连续信号傅里叶级数

(3)非周期性离散信号离散时域傅里叶变换

(4)周期性离散信号离散傅里叶变换

### 离散傅里叶变换（DFT）

傅里叶变换是连续函数的变换，而计算机只能处理离散数据。因此，为了在计算机上实现傅里叶变换，需要使用离散傅里叶变换（DFT）。

DFT将时域信号离散化为一组等间隔的采样点，并计算每个采样点的傅里叶变换。DFT的公式如下：

对于一个长度为 $N$ 的离散信号 $x[n]$，其DFT $X[k]$ 定义为：

$$
X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} kn}
$$

其中：
- X[k]$是信号在频域中的第 $k$ 个频率分量。
- x[n]是时间域中的第 $n$ 个信号样本。
- N是信号的总样本数。
- j 是虚数单位，满足 j^2 = -1。
- k是频率索引，取值范围为 0, 1, 2, ..., N-1。
- n是时间索引，取值范围为 0, 1, 2, ..., N-1。

#### 公式推导

1. 离散时间信号的表示

假设我们有一个长度为 $N$ 的离散时间信号 $x[n]$，其中 $n$ 的取值范围为 $0, 1, 2, \ldots, N-1$。

2. 傅里叶级数的基本概念

傅里叶级数表示一个周期信号为一系列正弦波和余弦波的叠加。对于连续信号，傅里叶级数的公式为：

$$
f(t) = \sum_{k=-\infty}^{\infty} C_k e^{j 2\pi k t / T}
$$

其中 $C_k$ 是傅里叶系数，$T$ 是信号的周期。

3. 离散傅里叶变换的定义

对于离散信号，我们希望找到类似的表示方法。我们将离散信号 $x[n]$ 看作是一个周期为 $N$ 的周期信号的一部分。离散傅里叶变换的目标是将这个离散信号表示为一系列复指数函数的叠加。

4. 复指数函数的性质

复指数函数 $e^{j \omega n}$ 可以表示为：

$$
e^{j \omega n} = \cos(\omega n) + j \sin(\omega n)
$$

其中 $\omega$ 是角频率。

5. DFT的推导

我们希望找到一组复指数函数，使得它们的线性组合可以表示原始信号 $x[n]$。假设我们有 $N$ 个这样的复指数函数，其频率分别为 $\frac{2\pi k}{N}$，其中 $k = 0, 1, 2, \ldots, N-1$。

我们假设信号 $x[n]$ 可以表示为这些复指数函数的线性组合：

$$
x[n] = \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} kn}
$$

其中 $X[k]$ 是我们需要找到的系数。

6. 求解系数 $X[k]$

为了求解 $X[k]$，我们利用正交性条件。将上式两边同时乘以 $e^{-j \frac{2\pi}{N} mn}$ 并对 $n$ 从 $0$ 到 $N-1$ 求和：

$$
\sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} mn} = \sum_{n=0}^{N-1} \left( \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} kn} \right) e^{-j \frac{2\pi}{N} mn}
$$

交换求和次序：

$$
\sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} mn} = \sum_{k=0}^{N-1} X[k] \sum_{n=0}^{N-1} e^{j \frac{2\pi}{N} (k-m)n}
$$

利用复指数函数的正交性：

$$
\sum_{n=0}^{N-1} e^{j \frac{2\pi}{N} (k-m)n} =
\begin{cases}
N, & \text{if } k = m \\
0, & \text{if } k \neq m
\end{cases}
$$

因此，当 $k = m$ 时：

$$
\sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn} = N X[k]
$$

解得：

$$
X[k] = \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn}
$$

为了简化表示，我们通常省略 $1/N$ 因子，得到标准形式的DFT公式：

$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn}
$$

7. 反向离散傅里叶变换（IDFT）

为了从频域信号 $X[k]$ 恢复时间域信号 $x[n]$，我们使用反向离散傅里叶变换，其公式为：

$$
x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} kn}
$$

这样，我们就完成了DFT公式的推导。

#### 反向傅里叶变换（IDFT）

反向离散傅里叶变换（Inverse Discrete Fourier Transform, IDFT）用于将频域信号转换回时间域信号，其公式为：

$$
x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \cdot e^{j \frac{2\pi}{N} kn}
$$

通过这些公式，可以在时间域和频域之间进行转换，分析信号的频率成分。

### 快速傅里叶变换（FFT）

DFT的计算复杂度为O(N^2)，对于大规模数据集来说非常耗时。FFT算法是一种快速计算DFT的算法，其复杂度为O(NlogN)。

FFT算法通过将DFT分解为一系列较小的DFT来实现加速。具体来说，FFT算法将N个采样点的DFT分解为两个N/2个采样点的DFT，然后递归地将N/2个采样点的DFT分解为更小的DFT，直到每个DFT只有两个采样点。

##### FFT算法的数学原理

FFT算法基于傅里叶变换，通过将DFT分解为一系列较小的DFT来实现快速计算。FFT算法的数学原理可以总结如下：

* 傅里叶变换将时域信号转换为频域信号，分解信号为一系列正弦波和余弦波。
* DFT是傅里叶变换的离散形式，将时域信号离散化为采样点并计算每个采样点的傅里叶变换。
* FFT算法通过将DFT分解为一系列较小的DFT来实现加速，其复杂度为O(NlogN)。

##### FFT算法的推导

为了更详细地推导从DFT到FFT的过程，我们可以将整个推导过程分为以下几个部分：

1. **DFT公式的回顾**
2. **信号的分解**
3. **DFT公式的分解**
4. **递归关系的建立**
5. **FFT算法的实现**

我们将逐步进行每个部分的推导和解释。首先，我们从第一部分开始。

###### 1. DFT公式的回顾

离散傅里叶变换（DFT）的公式如下：

$$
X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} kn}
$$

其中：
- $X[k]$ 是频域中的第 $k$ 个频率分量。
- $x[n]$ 是时间域中的第 $n$ 个信号样本。
- $N$ 是信号的总样本数。
- $j$ 是虚数单位，满足 $j^2 = -1$。
- $k$ 是频率索引，取值范围为 $0, 1, 2, \ldots, N-1$。
- $n$ 是时间索引，取值范围为 $0, 1, 2, \ldots, N-1$。

这个公式表示的是如何将一个长度为 $N$ 的离散时间信号 $x[n]$ 转换为频域信号 $X[k]$。

接下来，我们进入第二部分：信号的分解。

###### 2.信号的分解

为了高效计算DFT，我们可以利用信号的对称性和周期性，将信号分解为奇数和偶数两部分。假设信号 $x[n]$ 的长度为 $N$，我们将其分解为两部分：

- 偶数索引部分：$x[2m]$，其中 $m = 0, 1, 2, \ldots, \frac{N}{2}-1$
- 奇数索引部分：$x[2m+1]$，其中 $m = 0, 1, 2, \ldots, \frac{N}{2}-1$

这样，我们可以将原始信号 $x[n]$ 表示为：

$$
x[n] = \begin{cases}
x[2m], & \text{if } n = 2m \\
x[2m+1], & \text{if } n = 2m+1
\end{cases}
$$

通过这种分解，我们将原始的长度为 $N$ 的信号分解为两个长度为 $N/2$ 的信号。接下来，我们将利用这种分解来简化DFT的计算。

###### 3. DFT公式的分解

将DFT公式中的 $x[n]$ 替换为上述两部分：

$$
X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} kn}
$$

可以分解为：

$$
X[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot e^{-j \frac{2\pi}{N} k(2m)} + \sum_{m=0}^{N/2-1} x[2m+1] \cdot e^{-j \frac{2\pi}{N} k(2m+1)}
$$

进一步简化：

$$
X[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot e^{-j \frac{2\pi}{N/2} km} + e^{-j \frac{2\pi}{N} k} \sum_{m=0}^{N/2-1} x[2m+1] \cdot e^{-j \frac{2\pi}{N/2} km}
$$

我们可以定义两个新的DFT：

- $E[k]$ 表示偶数索引部分的DFT：
$$
E[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot e^{-j \frac{2\pi}{N/2} km}
$$

- $O[k]$ 表示奇数索引部分的DFT：
$$
O[k] = \sum_{m=0}^{N/2-1} x[2m+1] \cdot e^{-j \frac{2\pi}{N/2} km}
$$

于是，原始DFT公式可以表示为：

$$
X[k] = E[k] + e^{-j \frac{2\pi}{N} k} O[k]
$$

接下来，我们将进入第四部分：递归关系的建立。

###### 4. 递归关系的建立

在上一部分中，我们将DFT公式分解为两个长度为 $N/2$ 的DFT。现在，我们需要进一步分析这些分解后的DFT，以建立递归关系。

**分解后的DFT公式**

我们已经得到了以下公式：

$$
X[k] = E[k] + e^{-j \frac{2\pi}{N} k} O[k]
$$

其中：
- $E[k]$ 是偶数索引部分的DFT：
$$
E[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot e^{-j \frac{2\pi}{N/2} km}
$$

- $O[k]$ 是奇数索引部分的DFT：
$$
O[k] = \sum_{m=0}^{N/2-1} x[2m+1] \cdot e^{-j \frac{2\pi}{N/2} km}
$$

**频率索引的范围**

注意到 $k$ 的取值范围是 $0, 1, 2, \ldots, N-1$。为了进一步简化计算，我们将 $k$ 分为两部分：

- 当 $k = 0, 1, 2, \ldots, N/2-1$ 时：
$$
X[k] = E[k] + W_N^k O[k]
$$

- 当 $k = N/2, N/2+1, \ldots, N-1$ 时：
$$
X[k] = E[k - N/2] - W_N^{k - N/2} O[k - N/2