FFT 求特征值

快速傅里叶变换（FFT）不仅是信号处理的利器，也是数值线性代数的优雅工具。 本文讨论 FFT 与特征值之间的内在联系。

离散傅里叶变换

DFT 定义为：

$$ X_k = \sum_{n=0}^{N-1} x_n \, e^{-2\pi i k n / N}, \quad k = 0, 1, \dots, N-1 $$

行内公式示例：频率为 $k$ 的分量幅值 $|X_k|$ 由上式给出。 勾股定理是最基础的例子：$a^2 + b^2 = c^2$。

循环矩阵与特征值

DFT 矩阵 $F_N$ 的列向量是循环矩阵的特征向量。设 $C$ 为循环矩阵，则：

$$ C \, \mathbf{v}_k = \lambda_k \, \mathbf{v}_k $$

其中特征值 $\lambda_k$ 是 $C$ 第一行元素离散傅里叶变换的结果。更一般地，对厄米矩阵 $A = A^\dagger$，存在酉矩阵 $U$ 使得：

$$ A = U \Lambda U^\dagger, \quad \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n) $$

其中 $\lambda_i \in \mathbb{R}$。牛顿第二定律的矩阵形式也遵循类似结构：

$$ \mathbf{F} = m\mathbf{a} $$

复杂度对比

方法时间复杂度空间复杂度
直接 DFT$O(N^2)$$O(N)$
FFT (Cooley–Tukey)$O(N \log N)$$O(N)$
特征值分解$O(N^3)$$O(N^2)$

代码实现

#include <bits/stdc++.h>
using namespace std;

using cd = complex<double>;
const double PI = acos(-1);

void fft(vector<cd>& a, bool invert) {
    int n = a.size();
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    for (int len = 2; len <= n; len <<= 1) {
        double ang = 2 * PI / len * (invert ? -1 : 1);
        cd wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            cd w(1);
            for (int j = 0; j < len / 2; j++) {
                cd u = a[i + j], v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert) {
        for (cd& x : a) x /= n;
    }
}

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])
    factor = np.exp(-2j * np.pi * np.arange(n // 2) / n)
    return np.concatenate([even + factor * odd, even - factor * odd])

Bash 示例：

echo "hello" | awk '{print toupper($0)}'

频谱图

FFT 频谱示意

注意事项

使用 FFT 时注意数值精度：$N > 10^6$ 时相对误差约为 $10^{-12}$ 量级。

  1. 输入长度应为 2 的幂，否则先补零1
  2. 复数乘法顺序会影响精度
  3. 逆变换记得除以 $N$

参考

  • Cooley, J. W.; Tukey, J. W. (1965), An algorithm for the machine calculation of complex Fourier series
  • FFT 笔记

  1. Bluestein 算法可以处理任意长度，但常数较大。 ↩︎