频谱分析与 FFT 实战

引言

FFT（快速傅里叶变换）是音频工程里最常用的工具：频谱可视化、音高检测、滤波器设计验证、时间伸缩与变调、降噪、卷积混响，全部建立在它之上。它也是最容易被误用的工具——窗函数选错、分辨率不足、补零当精度用，都会得到看似合理但完全错误的结论。

工程上的核心认知是：FFT 的输出是一个"时频折中"的结果。窗越长，频率分辨率越高但时间分辨率越低；窗越短，能看清瞬态但看不清频率。这个不确定性原理（Heisenberg–Gabor 极限）无法绕过，只能根据任务选择合适的时间窗。

本文从 DFT 的定义出发，逐层拆解窗函数、分辨率、STFT、相位处理，然后进入频域处理（相位声码器）与特征提取（Mel、倒谱），最后给出实时频谱仪与常见算法的实现要点。

目录

1. DFT 与 FFT：复杂度与实现
2. 窗函数与谱泄漏
3. 频率分辨率与栅栏效应
4. STFT 与短时分析
5. 相位、群延迟与相位展开
6. 倒谱与 Mel 频谱
7. 实时频谱分析与可视化
8. 频域处理：相位声码器
9. 实现选型与性能优化

1. DFT 与 FFT：复杂度与实现

离散傅里叶变换（DFT）把 N 点时域序列映射到 N 点频域序列：

X[k] = Σ(n=0..N-1) x[n] · e^(-j2πkn/N)

直接计算需要 O(N²) 次复数乘加。FFT 利用旋转因子的对称性与周期性，把复杂度降到 O(N log N)。

1.1 基 2 与混合基

最常用的是 radix-2（N 为 2 的幂）。工程上的实现选择：

• KissFFT：轻量、任意 N（混合基），代码可读，适合嵌入式。
• FFTW：性能最强，支持任意 N 与 SIMD，GPL/商业双许可。
• pffft：仅支持 2 的幂与 3 的幂，SIMD 优化，BSD 许可。
• Ooura FFT：经典 C 实现，轻量，性能不错。
• Intel MKL / vDSP：平台专属，性能最好但不可移植。

// KissFFT 用法示例
#include "kiss_fft.h"
kiss_fft_cfg cfg = kiss_fft_alloc(N, 0, nullptr, nullptr);
kiss_fft_cpx in[N], out[N];
for (int i = 0; i < N; i++) { in[i].r = x[i]; in[i].i = 0; }
kiss_fft(cfg, in, out);          // 输出为复数，索引 0..N-1
kiss_fft_free(cfg);

1.2 实数输入的优化

音频是实数信号，其频谱共轭对称：X[N-k] = conj(X[k])。所以只需计算前 N/2+1 个点。实数 FFT（rfft） 通过把 N 点实数打包成 N/2 点复数，再用 N/2 点复数 FFT 实现，算力约为复数 FFT 的一半。

1.3 输出格式：幅度与 dB

工程上很少直接用复数输出，通常转成：

幅度：mag[k] = sqrt(re[k]² + im[k]²)
dBFS：db[k] = 20·log10(mag[k] / (N/2))     // 归一化到满量程

归一化因子 N/2 的由来：单频正弦的 FFT 峰值约为 N/2（其余能量在镜像频率）。若忘记归一化，得到的 dB 值会随 N 变化，无法跨窗长比较。

2. 窗函数与谱泄漏

FFT 假设输入是周期信号的一个整周期。现实信号不满足这个假设，截断处的不连续会产生谱泄漏（Spectral Leakage）：一个纯正弦的能量扩散到整个频谱。

加窗的作用是让截断处的幅度平滑过渡到 0，从而减少不连续。

2.1 主瓣与旁瓣的权衡

• 主瓣宽度决定频率分辨率：主瓣越宽，两个邻近频率越难分开。
• 旁瓣高度决定动态范围：强信号的旁瓣会掩盖邻近的弱信号。

没有同时最优的窗。测失真（看 -100 dB 的谐波）用 Blackman-Harris；测瞬态用矩形；通用分析用 Hann。

2.2 Hann 窗的公式

import numpy as np

def hann(N):
    return 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1))

def blackman_harris(N):
    a = [0.35875, 0.48829, 0.14128, 0.01168]
    n = np.arange(N)
    return (a[0] - a[1]*np.cos(2*np.pi*n/(N-1))
                 + a[2]*np.cos(4*np.pi*n/(N-1))
                 - a[3]*np.cos(6*np.pi*n/(N-1)))

2.3 加窗的幅度补偿

加窗会降低信号总能量，必须补偿。相干增益（coherent gain）是窗的平均值：

Hann 窗相干增益 = 0.5  → 幅度需乘 2 补偿

只做频谱可视化时可以忽略，但做绝对幅度测量（例如测响度、测谐波电平）必须补偿，否则结果偏低 6 dB。

3. 频率分辨率与栅栏效应

频率分辨率（bin 宽度）由窗长决定：

Δf = fs / N
N = 1024 @ 48 kHz  → Δf = 46.9 Hz
N = 4096 @ 48 kHz  → Δf = 11.7 Hz
N = 16384 @ 48 kHz → Δf = 2.9 Hz

3.1 栅栏效应与补零

FFT 只在 k·Δf 这些离散频率上给出值，中间的频率看不到——这就是栅栏效应。在峰值附近看不到真实最大值，会导致幅度测量偏低。

**补零（zero-padding）**在时域末尾补零把 N 增大，使频域采样点变密，从而"看"到峰值的真实形状：

X = np.fft.rfft(x * np.hanning(len(x)), n=8192)   # 4096 点补零到 8192

关键认知：补零不提高真实分辨率。它只是对已有的（由窗长决定的）连续频谱做更密的采样。分辨两个邻近频率的能力仍然由原始窗长决定。

3.2 分辨率的经验值

3.3 频率插值提高精度

不做长窗也能提高峰值频率精度：用峰值两侧的 bin 做抛物线插值，精度可提升 10~100 倍：

def interpolate_peak(mag, k):
    """用 log 幅度做抛物线插值，返回修正后的峰值位置"""
    if k <= 0 or k >= len(mag) - 1:
        return float(k)
    a, b, c = np.log(mag[k-1]), np.log(mag[k]), np.log(mag[k+1])
    denom = a - 2*b + c
    if abs(denom) < 1e-12:
        return float(k)
    return k + 0.5 * (a - c) / denom

这是吉他调音器、音高跟踪器等应用的标准技巧。

4. STFT 与短时分析

音乐信号是非平稳的：频率成分随时间变化。短时傅里叶变换（STFT）通过滑动窗把信号切成短段，对每段做 FFT，得到时频表示。

X[m, k] = Σ(n) x[n + m·H] · w[n] · e^(-j2πkn/N)

其中 H 是跳距（hop size），w[n] 是窗函数。

4.1 重叠与 COLA 条件

要能从 STFT 完美重建原信号，窗函数与跳距必须满足 COLA（Constant Overlap-Add） 条件：所有重叠窗之和为常数。

Hann 窗 + 50% 重叠（H = N/2）→ 满足 COLA
Hann 窗 + 75% 重叠（H = N/4）→ 满足 COLA（更平滑）
矩形窗 + 25% 重叠（H = N/4）→ 满足 COLA

Hann 窗 50% 重叠时，重叠窗之和恒为 1.0（在边界处归一化后）。这是最常用的配置。

4.2 重建

def stft(x, N=2048, H=512):
    w = np.hanning(N)
    frames = []
    for start in range(0, len(x) - N + 1, H):
        frames.append(np.fft.rfft(x[start:start+N] * w))
    return np.array(frames)

def istft(frames, N=2048, H=512, length=None):
    w = np.hanning(N)
    out = np.zeros(length or (len(frames) - 1) * H + N)
    wsum = np.zeros_like(out)
    for i, F in enumerate(frames):
        seg = np.fft.irfft(F, n=N) * w
        start = i * H
        out[start:start+N] += seg
        wsum[start:start+N] += w * w
    out /= np.maximum(wsum, 1e-8)      # 归一化，避免边界处幅度下降
    return out

忘记除 wsum 会导致重建信号在帧边界处幅度起伏（听感是周期性"抖动"）。

4.3 跳距与延迟

跳距越小，时间分辨率越高，但算力成倍增长。实时场景下，STFT 的延迟至少是一个窗长（要攒够 N 个样本才能做 FFT），因此：

N = 2048 @ 48 kHz → 至少 42.7 ms 延迟

这对实时效果器是不可接受的，所以实时处理更倾向用时域算法或小窗（N = 128~256）。

5. 相位、群延迟与相位展开

FFT 的相位输出被包裹在 (-π, π] 内，直接使用会看到跳变。

5.1 相位展开

def unwrap_phase(phase):
    return np.unwrap(phase)      # numpy 内置，处理 2π 跳变

5.2 群延迟

群延迟是相位对频率的负导数：

τ_g(ω) = -dφ(ω)/dω

它描述"某个频率成分延迟了多少"。对线性相位 FIR，群延迟恒定；对 IIR，群延迟随频率变化，这是 IIR 会"抹开瞬态"的数学解释。

5.3 相位在分析中的价值

• 音高检测：相位差可精确到亚 bin 精度。
• 瞬态检测：相位的一致性（phase coherence）可区分瞬态与稳态。
• 相位声码器：时间伸缩与变调的核心就是相位处理。

6. 倒谱与 Mel 频谱

6.1 倒谱（Cepstrum）

对频谱取对数再做一次 FFT（或 IFFT），得到倒谱：横轴是"quefrency"（类似时间的量）。

倒谱 = IFFT(log(|FFT(x)|))

用途：

• 基频检测：周期性信号的倒谱在基音周期处出现峰值。
• 语音分析：倒谱的低 quefrency 部分对应声道（谱包络），高 quefrency 部分对应声源（基频）。
• MFCC：语音识别的经典特征，就是 Mel 滤波器组 + 对数 + DCT。

6.2 Mel 频谱

Mel 刻度模拟人耳对频率的非线性感知：

mel(f) = 2595 · log10(1 + f/700)

Mel 滤波器组通常是 26~128 个三角滤波器，在低频密集、高频稀疏。计算流程：

def mel_spectrogram(x, sr=16000, n_fft=512, n_mels=40):
    frames = stft(x, N=n_fft, H=n_fft // 4)
    mag = np.abs(frames)
    fb = mel_filterbank(sr, n_fft, n_mels)      # (n_mels, n_fft//2+1)
    return np.log(fb @ mag.T + 1e-10)

Mel 频谱是语音识别、音频分类、音频指纹的通用输入特征。把这类特征送入神经网络（包括浏览器端的 WASM 推理）是当前的主流做法，见 wasm-ai-inference-browser 。

7. 实时频谱分析与可视化

7.1 用 AnalyserNode 做频谱

Web Audio 的 AnalyserNode 内置了 FFT 与窗函数：

const analyser = ctx.createAnalyser();
analyser.fftSize = 2048;                  // 频率分辨率 = sampleRate / fftSize
analyser.smoothingTimeConstant = 0.8;     // 时间平滑（0~1）
analyser.minDecibels = -90;
analyser.maxDecibels = -10;

const bins = new Uint8Array(analyser.frequencyBinCount);  // = fftSize / 2
function draw() {
  analyser.getByteFrequencyData(bins);    // 0~255 映射到 min~max dB
  // ... 绘制
  requestAnimationFrame(draw);
}

注意：

• getByteFrequencyData 返回的是非线性映射后的值（按 min/maxDecibels 映射到 0~255），要精确数值需用 getFloatFrequencyData。
• fftSize 必须是非零的 2 的幂，范围 32~32768。
• 平滑是通过对连续帧做指数平均实现的，smoothingTimeConstant 越大越平滑但响应越慢。

7.2 频谱图（Spectrogram）绘制

频谱图把每帧的频谱竖着画成一条，横向拼出时间轴。绘制性能是关键：

// 用 ImageData 直接写像素，避免每帧重建 canvas
const canvas = document.createElement('canvas');
const ctx2d = canvas.getContext('2d', { willReadFrequently: false });
const img = ctx2d.createImageData(width, height);
// 每帧：把新的一列写入，然后整体左移（或用滚动绘制）

更高效的做法是维护一个离屏的滚动缓冲，只更新变化的部分。

7.3 对数频率轴

线性频率轴在音乐场景下不好用（低频挤在一起）。改用对数轴：横轴坐标 x = log(f/f_min) / log(f_max/f_min)，需要按 bin 反查插值。

7.4 时间平滑与峰值保持

专业频谱仪提供三种显示模式：

• 瞬时：无平滑，能看瞬态但闪烁。
• 平均：多帧平均，噪声底稳定，适合测失真。
• 峰值保持：记录每 bin 的历史最大值，适合找共振峰。

8. 频域处理：相位声码器

相位声码器（Phase Vocoder）用 STFT 修改信号的时频结构，实现时间伸缩（time stretching）与变调（pitch shifting）。

8.1 基本流程

1. 分析：对每帧做 STFT，得到幅度 |X| 与相位 φ
2. 修改：按伸缩比改变跳距（分析跳距 ≠ 合成跳距）
3. 相位修正：累积相位差并展开，避免帧间相位不连续
4. 合成：用修改后的幅度与相位做 IFFT，重叠相加

核心是第 3 步：直接复制相位会导致帧间相位跳变，产生"机器人声"或"金属味"。正确做法是计算瞬时频率并累积：

def phase_vocoder(D, rate, N=2048, H=512):
    """D: STFT 复数矩阵；rate > 1 加速（变短）"""
    time_steps = np.arange(0, D.shape[0], rate)
    out = np.zeros((len(time_steps), D.shape[1]), dtype=complex)
    phase_acc = np.angle(D[0])
    for i, t in enumerate(time_steps):
        i0, i1 = int(t), min(int(t) + 1, D.shape[0] - 1)
        frac = t - i0
        mag = (1 - frac) * np.abs(D[i0]) + frac * np.abs(D[i1])
        dphi = np.angle(D[i1]) - np.angle(D[i0]) - H * 2*np.pi*np.arange(N//2+1)/N
        dphi = np.mod(dphi + np.pi, 2*np.pi) - np.pi     # 包裹到 [-π, π)
        phase_acc += dphi + H * 2*np.pi*np.arange(N//2+1)/N
        out[i] = mag * np.exp(1j * phase_acc)
    return out

8.2 变调与时间伸缩的组合

变调 = 时间伸缩 + 重采样。先把信号时间伸缩 r 倍，再以 1/r 倍速重采样回原长度，音高就变了 r 倍。

变调 +12 半音 = 时间伸缩 2× + 重采样 0.5×

8.3 质量与替代方案

相位声码器的固有问题是瞬态被涂抹（smearing）：鼓点会变得模糊。改进方案：

• 瞬态检测 + 相位重置：检测到瞬态时重置相位累积。
• WSOLA（波形相似叠加）：时域方法，瞬态保持更好，但变调质量不如相位声码器。
• Phase Vocoder with Identity Phase Locking：锁定峰值周围的相位关系，显著改善音乐质量。

9. 实现选型与性能优化

9.1 库选型

9.2 性能优化要点

• 复用 FFT 计划：FFTW 的 plan 创建开销很大，必须缓存复用。
• 预计算旋转因子与窗：不要在每次调用时重算。
• 实数 FFT：算力减半。
• 避免内存分配：实时路径上所有缓冲预分配，模式与 audio-worklet-realtime 中的实时安全规则一致。
• 批处理：多帧一次做 FFT 可以利用缓存局部性。

// FFTW 计划复用
static fftwf_plan plan = nullptr;
if (!plan) {
    plan = fftwf_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);   // 只创建一次
}
fftwf_execute_dft_r2c(plan, in, out);

9.3 精度问题

• float vs double：音频分析通常 float 足够（约 7 位有效数字），但做长累积（如平均 1000 帧）时 double 更稳。
• DC 与 Nyquist bin：实数 FFT 的 DC（k=0）与 Nyquist（k=N/2）bin 是实数，处理时要特殊对待，否则 IFFT 结果出现虚部残留。
• 溢出：定点 FFT 需要逐级缩放（每个 radix-2 级右移 1 位），否则中间结果溢出。

权衡取舍

常见坑清单

1. 不加窗直接 FFT：谱泄漏使纯正弦的能量扩散到全频段，看不出真实频谱。
2. 把补零当分辨率：补零只加密频域采样，分辨两个邻近频率的能力仍由原始窗长决定。
3. 忘记幅度归一化：N/2 因子遗漏导致 dB 值随窗长变化，无法跨配置比较。
4. 忘记窗的相干增益补偿：Hann 窗不补偿幅度会偏低 6 dB，做绝对测量时结论全错。
5. STFT 重建不除窗和：帧边界幅度起伏，听到周期性"抖动"，需除以 Σw²。
6. 相位直接复制：相位声码器不做相位累积会产生"机器人声"，必须计算瞬时频率。
7. 实时链路用大窗：2048 点窗 @ 48 kHz 意味着 42 ms 延迟，实时效果器不可用。
8. DC/Nyquist bin 当复数处理：实数 FFT 的这两个 bin 是实数，处理不当会残留虚部。

小结

FFT 的工程要点可以概括为三条：窗函数决定动态范围，窗长决定分辨率，跳距决定时间精度。这三者构成一个无法同时最优的三角，必须按任务取舍。补零、插值、平滑都是"显示层"的优化，不会改变物理分辨率。

频域处理（相位声码器）打开了时间伸缩、变调、降噪、音源分离的大门，但代价是瞬态被涂抹。选择时要在"频域灵活"与"时域保真"之间权衡，必要时用混合方案（瞬态检测 + 相位重置）。

继续深入建议读 audio-dsp-filters 把频域分析用于滤波器设计与验证，读 audio-worklet-realtime 在浏览器实时链路中落地 STFT，读 audio-quality-testing 了解如何用频谱指标做自动化音频回归。