21、离散傅里叶变换（DFT）与快速傅里叶变换（FFT）详解

离散傅里叶变换（DFT）与快速傅里叶变换（FFT）详解

1. 离散傅里叶变换（DFT）的应用与相关问题

在信号处理领域，离散傅里叶变换（DFT）是一个非常重要的工具，它可以将时域信号转换为频域信号，从而方便我们对信号的频率成分进行分析。下面我们来看一些DFT的应用实例和相关问题。

1.1 频率样本间距计算

假设有一个连续时间信号 $x_c(t)$，以 $f_s = 2$ kHz 的采样频率进行采样。如果对 1000 个采样点进行 1000 点的 DFT 计算，那么频率样本 $X(k)$ 在模拟频率上的间距 $\Delta f$ 为 2 kHz。

1.2 加窗序列的 DFT 计算

已知 $x(n)$ 的 $N$ 点 DFT 为 $X(k)$，对于加窗序列 $y(n) = w(n)x(n)$，其中 $w(n)$ 是布莱克曼窗（Blackman window），其表达式为：
[w(n) = 0.42 - 0.5\cos\left(\frac{2\pi n}{N - 1}\right) + 0.08\cos\left(\frac{4\pi n}{N - 1}\right),\quad 0\leq n\leq N - 1]
则 $y(n)$ 的 $N$ 点 DFT $Y(k)$ 可通过以下公式计算：
[Y(k) = 0.42X(k) - 0.25[X((k - 1))_N + X((k + 1))_N] + 0.04[X((k - 2))_N + X((k + 2))_N]]

2. 快速傅里叶变换（FFT）概述

虽然 DFT 非常有用，但直接计算 DFT 的计算复杂度较高，尤其是当序列长度 $N$ 较大时。为了提高计算效率，人们开发了快速傅里叶变换（FFT）算法。FFT 算法的基本策略是“分而治之”，即将一个 $N$ 点的 DFT 分解为多个更小的 DFT，从而减少计算量。

3. 基 - 2 FFT 算法

基 - 2 FFT 算法是最常见的 FFT 算法之一，它要求序列长度 $N$ 是 2 的幂次方，即 $N = 2^v$。下面我们详细介绍基 - 2 FFT 算法的两种实现方式：按时间抽取（Decimation - in - Time）FFT 和按频率抽取（Decimation - in - Frequency）FFT。

3.1 按时间抽取（Decimation - in - Time）FFT

按时间抽取 FFT 算法的核心思想是将输入序列 $x(n)$ 分解为更小的子序列，然后通过这些子序列的 DFT 来计算原序列的 DFT。具体步骤如下：
- 序列分解 ：设 $x(n)$ 是长度为 $N = 2^v$ 的序列，将其分解为两个长度为 $N/2$ 的子序列。第一个子序列 $g(n)$ 由 $x(n)$ 的偶数索引项组成，第二个子序列 $h(n)$ 由 $x(n)$ 的奇数索引项组成，即：
[g(n) = x(2n),\quad n = 0, 1, \cdots, \frac{N}{2} - 1]
[h(n) = x(2n + 1),\quad n = 0, 1, \cdots, \frac{N}{2} - 1]
- DFT 计算 ：原序列 $x(n)$ 的 $N$ 点 DFT $X(k)$ 可以表示为：
[X(k) = G(k) + W_N^kH(k),\quad k = 0, 1, \cdots, N - 1]
其中 $G(k)$ 是 $g(n)$ 的 $N/2$ 点 DFT，$H(k)$ 是 $h(n)$ 的 $N/2$ 点 DFT，$W_N^k = e^{-j\frac{2\pi k}{N}}$ 称为旋转因子（twiddle factor）。

下面是一个 8 点按时间抽取 FFT 的流程图示例：

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px

    A([输入序列 x(n)]):::startend --> B(分解为 g(n) 和 h(n)):::process
    B --> C(计算 G(k) 和 H(k)):::process
    C --> D(计算 X(k) = G(k) + W_N^kH(k)):::process
    D --> E([输出 X(k)]):::startend

按时间抽取 FFT 算法的优点是可以实现原位计算，即只需要一个长度为 $N$ 的数组来存储中间结果。但输入序列需要按位反转的顺序存储。

3.2 按频率抽取（Decimation - in - Frequency）FFT

按频率抽取 FFT 算法则是将输出序列 $X(k)$ 分解为更小的子序列，通过对输入序列进行适当的组合和变换来计算这些子序列。具体步骤如下：
- 样本分离 ：分别计算 $X(k)$ 的偶数索引和奇数索引样本。偶数样本 $X(2k)$ 可以通过将输入序列的前 $N/2$ 点和后 $N/2$ 点相加后进行 $N/2$ 点 DFT 得到；奇数样本 $X(2k + 1)$ 也可以通过类似的方式计算。
- 递归分解 ：重复上述过程，直到最终得到长度为 2 的 DFT。

下面是一个 8 点按频率抽取 FFT 的流程图示例：

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px

    A([输入序列 x(n)]):::startend --> B(分离 X(2k) 和 X(2k + 1)):::process
    B --> C(计算相关序列的 N/2 点 DFT):::process
    C --> D(递归分解):::process
    D --> E([输出 X(k)]):::startend

按频率抽取 FFT 算法的复杂度与按时间抽取 FFT 相同，也可以实现原位计算，但输出序列需要按位反转的顺序读取。

4. 复合长度序列的 FFT 算法

在实际应用中，序列的长度 $N$ 不一定是 2 的幂次方。但如果 $N$ 可以分解为多个因子的乘积，即 $N = N_1 \cdot N_2$，我们仍然可以使用高效的 FFT 算法。具体步骤如下：
- 序列分解 ：将序列 $x(n)$ 分解为 $N_2$ 个长度为 $N_1$ 的子序列，并将这些子序列排列成一个二维数组。
- 行 DFT 计算 ：对二维数组的每一行进行 $N_1$ 点 DFT 计算，得到一个新的二维数组。
- 旋转因子相乘 ：将新数组的每个元素乘以相应的旋转因子。
- 列 DFT 计算 ：对二维数组的每一列进行 $N_2$ 点 DFT 计算，最终得到 $N$ 点 DFT 的结果。

例如，对于长度为 $N = 15$ 的序列，我们可以将其分解为 5 个长度为 3 的子序列，或者 3 个长度为 5 的子序列，然后按照上述步骤进行计算。

5. 素因子 FFT 算法

对于某些特定的 $N$ 值，通过适当的索引映射，可以完全消除旋转因子。如果 $N = N_1 \cdot N_2$，且 $N_1$ 和 $N_2$ 互质（即它们没有共同的因子），则可以找到一组合适的整数 $A$、$B$、$C$ 和 $D$，使得旋转因子可以被消除。

例如，对于 $N = 12$，$N_1 = 3$，$N_2 = 4$，可以通过特定的索引映射将 12 点 DFT 分解为 3 点 DFT 和 4 点 DFT 的组合，从而减少计算量。

6. 快速傅里叶变换（FFT）的计算复杂度分析

为了更直观地了解 FFT 算法的优势，我们来分析一下不同长度 DFT 计算的时间复杂度。

DFT 长度	直接计算乘法次数	FFT 计算乘法次数	时间节省比例
1024	$1024^2 = 1048576$	$512 \cdot \log_2{1024} = 5120$	约 99.5%
4096	$4096^2 = 16777216$	$2048 \cdot \log_2{4096} = 24576$	约 99.9%

从表格中可以看出，随着 DFT 长度的增加，FFT 算法相对于直接计算 DFT 的优势越来越明显。

7. 快速傅里叶变换（FFT）的实际应用案例

下面我们通过几个实际应用案例来进一步说明 FFT 算法的重要性。

7.1 信号卷积

在信号处理中，卷积是一个常见的操作。对于长度为 $N$ 的序列 $x(n)$ 和长度为 $L$ 的序列 $h(n)$，直接进行卷积需要 $N \cdot L$ 次复数乘法。而使用重叠相加法（overlap - add method）结合 1024 点 FFT 进行卷积，则可以大大减少计算量。

例如，对于长度为 $N = 8192$ 的序列 $x(n)$ 和长度为 $L = 512$ 的序列 $h(n)$，直接卷积需要 $8192 \cdot 512 = 4194304$ 次复数乘法；而使用重叠相加法结合 1024 点 FFT 只需要约 $4.5\%$ 的计算量。

7.2 实时语音处理

在实时语音处理中，需要对采样后的语音信号进行快速处理。假设语音信号以 10 kHz 的采样率进行采样，每次处理 1024 个样本。使用 FFT 算法计算 1024 点 DFT 和 1024 点逆 DFT 后，还可以剩余 61.44 ms 的时间进行其他处理，从而保证了实时性。

8. 快速傅里叶变换（FFT）的优化技巧

为了进一步提高 FFT 算法的效率，我们可以采用以下优化技巧：

8.1 排除乘以 $\pm 1$ 的乘法

在按时间抽取和按频率抽取 FFT 算法中，有些旋转因子的值为 $\pm 1$，可以通过编写专门的程序来排除这些乘法，从而减少计算量。

例如，对于 8 点按时间抽取 FFT，排除乘以 $\pm 1$ 的乘法后，实际只需要 5 次复数乘法；对于 16 点按时间抽取 FFT，需要 17 次复数乘法。

8.2 优化复数乘法

复数乘法通常需要 4 次实数乘法和 3 次实数加法。但可以通过一些技巧将其优化为 3 次实数乘法和 4 次实数加法。

设两个复数 $a + jb$ 和 $c + jd$，它们的乘积 $(a + jb)(c + jd) = (ac - bd) + j(ad + bc)$。可以通过以下方式优化计算：
[p_1 = a(c + d)]
[p_2 = d(a + b)]
[p_3 = c(b - a)]
则乘积的实部为 $p_1 - p_2$，虚部为 $p_1 + p_3$。

9. 利用 FFT 计算不同类型 DFT 的方法

在实际应用中，我们可能会遇到各种不同类型的 DFT 计算需求，下面介绍几种利用 FFT 计算不同类型 DFT 的方法。

9.1 计算两个实值序列的 N 点 DFT

如果有两个实值序列 $x_1(n)$ 和 $x_2(n)$，可以将它们组合成一个复数序列 $x(n) = x_1(n) + jx_2(n)$，然后计算 $x(n)$ 的 $N$ 点 DFT $X(k)$。最后，通过 $X(k)$ 的共轭对称和共轭反对称部分分别得到 $x_1(n)$ 和 $x_2(n)$ 的 DFT。具体公式如下：
[X_1(k) = \frac{1}{2}[X(k) + X^ ((N - k))_N]]
[X_2(k) = \frac{1}{2j}[X(k) - X^ ((N - k))_N]]

9.2 计算 2N 点实值序列的 DFT

对于长度为 $2N$ 的实值序列 $g(n)$，可以将其分解为两个长度为 $N$ 的实值序列 $x_1(n)$ 和 $x_2(n)$，然后组合成复数序列 $x(n) = x_1(n) + jx_2(n)$。计算 $x(n)$ 的 $N$ 点 DFT $X(k)$ 后，再通过以下公式得到 $g(n)$ 的 $2N$ 点 DFT $G(k)$：
[G(k) = X_1(k) + W_{2N}^kX_2(k),\quad k = 0, 1, \cdots, 2N - 1]
其中 $X_1(k)$ 和 $X_2(k)$ 可以通过上述方法从 $X(k)$ 中提取。

9.3 利用 FFT 计算逆 DFT

已知一个 FFT 程序可以计算序列的 $N$ 点 DFT，我们可以利用它来计算逆 DFT。方法有两种：
- 方法一 ：计算 $X^ (k)$ 的 DFT，然后将结果取共轭并除以 $N$ 得到 $x(n)$。
- 方法二 *：计算 $X(k)$ 的 DFT，然后将结果进行适当的变换得到 $x(n)$。

10. 特殊序列的 DFT 特性及计算方法

在某些情况下，序列具有特殊的性质，我们可以利用这些性质来简化 DFT 的计算。

10.1 只含奇次谐波的序列

设序列 $x(n)$ 满足 $x(n) = -x(n + N/2)$，其中 $N$ 为偶数。则该序列的 $N$ 点 DFT $X(k)$ 只包含奇次谐波，即 $X(k) = 0$ （$k$ 为偶数）。

为了计算该序列的 $N$ 点 DFT，我们可以将前 $N/2$ 点的 $x(n)$ 乘以 $2W_N^k$ 得到序列 $y(n)$，然后计算 $y(n)$ 的 $N/2$ 点 DFT 即可。

总结

快速傅里叶变换（FFT）是信号处理领域中一个非常重要的算法，它通过“分而治之”的策略大大减少了离散傅里叶变换（DFT）的计算量。本文详细介绍了基 - 2 FFT 算法（包括按时间抽取和按频率抽取）、复合长度序列的 FFT 算法、素因子 FFT 算法等，并通过实际应用案例和计算复杂度分析展示了 FFT 算法的优势。同时，还介绍了一些 FFT 算法的优化技巧和利用 FFT 计算不同类型 DFT 的方法，以及特殊序列 DFT 的特性和计算方法。希望通过本文的介绍，读者能够对 DFT 和 FFT 有更深入的理解，并在实际应用中灵活运用这些算法。

离散傅里叶变换（DFT）与快速傅里叶变换（FFT）详解

11. FFT 算法的硬件实现考虑

在实际工程应用中，除了软件实现 FFT 算法，还常常会涉及到硬件实现，以满足实时性和高性能的要求。以下是一些硬件实现 FFT 算法时的重要考虑因素。

11.1 资源占用

硬件实现 FFT 算法需要消耗一定的硬件资源，如逻辑门、寄存器和乘法器等。不同的 FFT 算法和实现方式对资源的需求差异较大。例如，按时间抽取和按频率抽取的基 - 2 FFT 算法在硬件实现时，其资源占用情况会受到序列长度、数据位宽等因素的影响。为了降低资源占用，可以采用一些优化技术，如复用乘法器、共享寄存器等。

11.2 并行性

FFT 算法具有一定的并行性，可以通过并行计算来提高计算速度。在硬件实现中，可以采用多处理器并行计算、流水线技术等方法来充分利用并行性。例如，在大规模 FFT 计算中，可以将序列分成多个子序列，由多个处理器同时进行计算，最后将结果合并。

11.3 数据存储

FFT 算法在计算过程中需要大量的数据存储，包括输入序列、中间结果和旋转因子等。在硬件实现时，需要合理设计数据存储结构，以提高数据读取和写入的效率。例如，可以采用双端口 RAM 来实现数据的并行读写，避免数据访问冲突。

12. FFT 算法在不同领域的应用拓展

FFT 算法不仅在信号处理领域有着广泛的应用，还在其他许多领域发挥着重要作用。

12.1 图像处理

在图像处理中，FFT 算法可以用于图像的频域滤波、图像压缩等。通过将图像从空间域转换到频域，可以方便地对图像的频率成分进行分析和处理。例如，通过低通滤波可以去除图像中的高频噪声，通过高通滤波可以增强图像的边缘信息。

12.2 通信系统

在通信系统中，FFT 算法是正交频分复用（OFDM）技术的核心。OFDM 技术通过将高速数据信号分成多个低速子信号，在多个子载波上并行传输，从而提高了频谱利用率和抗干扰能力。FFT 算法用于将时域的 OFDM 信号转换为频域信号，以便进行调制和解调。

12.3 雷达系统

在雷达系统中，FFT 算法可以用于目标检测和测距。通过对雷达回波信号进行 FFT 变换，可以得到信号的频谱信息，从而分析目标的速度和距离。例如，在脉冲多普勒雷达中，FFT 算法用于计算多普勒频移，从而测量目标的径向速度。

13. FFT 算法的发展趋势

随着科技的不断发展，FFT 算法也在不断演进和完善。以下是一些 FFT 算法的发展趋势。

13.1 更高的计算效率

为了满足日益增长的大数据处理需求，FFT 算法需要不断提高计算效率。未来的 FFT 算法可能会采用更加先进的并行计算技术、优化的硬件架构和算法设计，以实现更快的计算速度和更低的能耗。

13.2 适应不同的序列长度

目前的 FFT 算法主要适用于序列长度为 2 的幂次方或可分解为多个因子乘积的情况。未来的 FFT 算法可能会更加灵活，能够直接处理任意长度的序列，而无需进行复杂的序列扩展或分解。

13.3 与其他算法的融合

FFT 算法可能会与其他信号处理算法和机器学习算法进行融合，以实现更加复杂和智能的信号处理任务。例如，将 FFT 算法与深度学习算法相结合，可以用于语音识别、图像分类等领域，提高算法的性能和准确性。

14. 常见问题解答

在使用 FFT 算法时，可能会遇到一些常见的问题，以下是对这些问题的解答。

14.1 为什么 FFT 算法要求序列长度为 2 的幂次方？

基 - 2 FFT 算法要求序列长度为 2 的幂次方，是因为该算法基于“分而治之”的思想，通过不断将序列分解为更小的子序列来减少计算量。当序列长度为 2 的幂次方时，可以方便地进行递归分解，从而实现高效的计算。如果序列长度不是 2 的幂次方，可以通过补零的方式将序列扩展为 2 的幂次方长度。

14.2 如何选择合适的 FFT 算法？

选择合适的 FFT 算法需要考虑多个因素，如序列长度、计算资源、实时性要求等。如果序列长度是 2 的幂次方，基 - 2 FFT 算法是一个不错的选择，因为它具有较高的计算效率和简单的实现方式。如果序列长度可以分解为多个因子的乘积，可以选择复合长度序列的 FFT 算法。对于某些特定的序列长度，素因子 FFT 算法可能会更加高效。

14.3 FFT 算法的计算结果是否准确？

FFT 算法的计算结果在理论上是准确的，但在实际应用中，由于有限字长效应和数值计算误差等因素的影响，可能会存在一定的误差。为了减少误差，可以采用更高的数据位宽、优化的算法实现和数值计算方法等。

15. 总结与展望

离散傅里叶变换（DFT）是信号处理领域的基础工具，但直接计算 DFT 的计算复杂度较高。快速傅里叶变换（FFT）算法通过“分而治之”的策略，将一个 $N$ 点的 DFT 分解为多个更小的 DFT，大大减少了计算量，提高了计算效率。

本文详细介绍了多种 FFT 算法，包括基 - 2 FFT 算法（按时间抽取和按频率抽取）、复合长度序列的 FFT 算法、素因子 FFT 算法等，还探讨了 FFT 算法的计算复杂度、实际应用案例、优化技巧以及在不同领域的应用拓展。同时，对 FFT 算法的硬件实现考虑、发展趋势和常见问题进行了解答。

随着科技的不断进步，FFT 算法将在更多领域发挥重要作用，并且会不断发展和完善。未来，我们可以期待更加高效、灵活和智能的 FFT 算法的出现，为信号处理和其他相关领域带来新的突破和发展。

希望本文能够帮助读者深入理解 DFT 和 FFT 算法，并在实际应用中合理选择和使用这些算法，解决各种复杂的信号处理问题。

以下是一个简单的 FFT 算法在 Python 中的实现示例：

import numpy as np

def fft(x):
    N = len(x)
    if N <= 1:
        return x
    even = fft(x[0::2])
    odd = fft(x[1::2])
    T = [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)]
    return [even[k] + T[k] for k in range(N // 2)] + [even[k] - T[k] for k in range(N // 2)]

# 示例使用
x = [1, 2, 3, 4]
X = fft(x)
print("FFT 结果:", X)

这个示例代码实现了一个简单的按时间抽取的基 - 2 FFT 算法。通过递归调用 fft 函数，将序列不断分解为更小的子序列，最终计算出 FFT 的结果。

通过学习和掌握 FFT 算法，我们可以更好地处理和分析各种信号，为解决实际问题提供有力的支持。

以下是一个总结 FFT 算法不同类型及其特点的表格：
| FFT 算法类型 | 序列长度要求 | 计算复杂度 | 原位计算 | 输入/输出顺序 |
| ---- | ---- | ---- | ---- | ---- |
| 按时间抽取基 - 2 FFT | 2 的幂次方 | $O(N\log_2N)$ | 可以 | 输入按位反转 |
| 按频率抽取基 - 2 FFT | 2 的幂次方 | $O(N\log_2N)$ | 可以 | 输出按位反转 |
| 复合长度序列 FFT | 可分解为因子乘积 | 取决于因子 | 部分可以 | 无特殊要求 |
| 素因子 FFT | 特定 $N$ 值（$N = N_1 \cdot N_2$ 且 $N_1$ 和 $N_2$ 互质） | 较低 | 可以 | 无特殊要求 |

通过这个表格，可以更直观地比较不同 FFT 算法的特点，以便在实际应用中做出合适的选择。

最后，为了更清晰地展示 FFT 算法在整个信号处理流程中的位置，以下是一个 mermaid 流程图：

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px

    A([输入信号]):::startend --> B(采样):::process
    B --> C(加窗):::process
    C --> D(FFT 计算):::process
    D --> E(频域处理):::process
    E --> F(逆 FFT 计算):::process
    F --> G([输出信号]):::startend

这个流程图展示了一个典型的信号处理流程，其中 FFT 计算和逆 FFT 计算是关键步骤。通过 FFT 计算将信号从时域转换到频域，在频域进行处理后，再通过逆 FFT 计算将信号转换回时域，最终得到处理后的输出信号。