信号与频谱：傅里叶级数、傅里叶变换与 FFT

1. 时域与频域

先看一个信号：

画成波形是一条弯弯曲曲的线，很难一眼看出它由什么组成；换成频谱，只有三根竖线：0 Hz（直流）处 0.2，50 Hz 处 1.0，120 Hz 处 0.3。

为什么偏偏拆成正弦波？因为线性系统对正弦输入的稳态输出还是同频率的正弦，只是幅值和相位变了（第5章 5.1 节会证明）。把输入拆成一堆正弦，分别算出每个正弦经过系统后变成什么，再加起来，就是输出。这就是频域方法的出发点。

2. 傅里叶级数：周期信号 = 直流 + 基波 + 谐波

周期为 \(T\) 的信号（基波频率 \(f_0=1/T\)，角频率 \(\omega_0=2\pi f_0\)）可以写成

\(a_0\) 是直流分量（平均值），\(n=1\) 的项叫基波，\(n\ge2\) 的项叫 \(n\) 次谐波，频率是 \(nf_0\)。系数这样算：

同频率的 cos 和 sin 可以合成一个带初相位的正弦：\(a_n\cos n\omega_0t+b_n\sin n\omega_0t=A_n\sin(n\omega_0t+\varphi_n)\)，幅值 \(A_n=\sqrt{a_n^2+b_n^2}\)。把 \(A_n\) 画在频率 \(nf_0\) 处，就是幅值谱：周期信号的频谱只在基波的整数倍上有值，是一根根离散的谱线。

例：方波

幅值 ±1 的方波：前半个周期为 +1，后半个周期为 −1。平均值为 0，所以 \(a_0=0\)；它是奇函数（\(x(-t)=-x(t)\)），所以所有 \(a_n=0\)。只需要算 \(b_n\)：

基波幅值 \(4/\pi\approx1.273\)，比方波本身的幅值 1 还大；3 次谐波 0.424，5 次 0.255，7 次 0.182。在下面的图里把谐波一项一项加上去：

两个和工作有关的例子：

• PWM：占空比为 \(D\) 的 PWM 电压 = 直流 \(DU\) + PWM 频率和它的各次谐波。电机绕组的电感对高频的阻抗很大，把谐波挡在外面，电流里只剩下直流部分和一点纹波。所以 PWM 能当「可调电压」用。
• 描述函数（第12章）：继电器输出幅值为 \(M\) 的方波，基波幅值就是 \(4M/\pi\)。描述函数 \(N(A)=\dfrac{4M}{\pi A}\) 只保留了基波，因为高次谐波大多被对象的低通特性滤掉了。

复指数形式（为 FFT 做准备）

由欧拉公式 \(\cos\theta=\frac12(e^{j\theta}+e^{-j\theta})\)、\(\sin\theta=\frac1{2j}(e^{j\theta}-e^{-j\theta})\)，级数可以改写成更紧凑的形式：

这里出现了「负频率」\(n<0\)。它不是什么物理现象，只是记账方式：一个实正弦被拆成 \(+nf_0\) 和 \(-nf_0\) 两个复指数，各分走一半幅值。对实信号，\(c_{-n}\) 和 \(c_n\) 共轭，\(|c_n|=\frac12A_n\)（\(n\ne0\)），\(c_0=a_0\)。记住这个「一半」，后面 FFT 的幅值换算要用。

3. 傅里叶变换：不是周期的信号

一个脉冲、一次冲击、一段阶跃响应都不是周期信号。办法是：把它看成周期 \(T\to\infty\) 的周期信号。

以宽度 \(\tau\)、高度 1 的矩形脉冲串为例，周期为 \(T\)。按上一节的公式算出（\(\operatorname{sinc}x=\dfrac{\sin\pi x}{\pi x}\)）：

谱线间隔 \(1/T\)。周期越长，谱线越密，每根也越矮。但如果把纵轴乘以 \(T\)，\(Tc_n=\tau\operatorname{sinc}(f\tau)\big|_{f=n/T}\)，所有谱线都落在同一条曲线上，和 \(T\) 无关。拖动周期试试：

\(T\to\infty\) 时只剩下一个脉冲，谱线连成那条曲线，求和变成积分，这就是傅里叶变换：

用频率 \(f\) 写也一样：\(X(f)=\int x(t)e^{-j2\pi ft}dt\)。

例　宽 \(\tau\)、高 1、中心在 0 的单个矩形脉冲：

第一个零点在 \(f=1/\tau\)：脉冲越窄，频谱越宽。1 ms 的脉冲，能量主要在 1 kHz 以内；1 µs 的脉冲，频谱一直铺到 1 MHz。极限情况是冲激 \(\delta(t)\)，它的傅里叶变换在所有频率上都等于 1。所以用小锤敲一下机械结构（近似冲激），就能一次激励出它所有的共振频率，这是测机械共振的常用办法。

单位不一样：级数的系数 \(c_n\) 和信号的单位相同（比如 V）；\(X(f)\) 的单位是「信号 × 秒」（V·s，也就是 V/Hz），它是密度。连续谱上单独一个频率没有「幅值」可言，一段频带才有能量。噪声就是这种连续谱，所以要用功率谱密度描述（第5节）。

口头上「傅里叶变换」常被用来泛指这几种。工作中说「做个傅里叶变换看看」，几乎都是指对采样数据做 FFT。

4. 和拉氏变换、频率特性的关系

拉氏变换 \(F(s)=\int_0^\infty f(t)e^{-st}dt\)，\(s=\sigma+j\omega\)。对 \(t<0\) 时为 0 的信号，令 \(\sigma=0\)，它就变成了傅里叶变换。拉氏变换多出来的 \(e^{-\sigma t}\) 是一个衰减因子，让阶跃、斜坡甚至发散的信号也能积分；它还能自然地带上初始条件。所以分析控制系统（关心过渡过程和稳定性）用拉氏变换，分析信号（关心含有哪些频率）用傅里叶变换。

信号经过系统：时域里输出是输入和脉冲响应的卷积（基础2 的卷积定理），频域里就是相乘：

\(G(j\omega)=G(s)\big|_{s=j\omega}\) 就是第5章的频率特性：输入里频率为 \(\omega\) 的分量，幅值被放大 \(|G(j\omega)|\) 倍，相位移动 \(\angle G(j\omega)\)。

例　100 Hz、±1 的方波，经过截止频率 100 Hz 的一阶低通 \(G(s)=\dfrac{1}{\tau s+1}\)（\(\tau=\frac{1}{2\pi\cdot100}\)）。第 \(n\) 次谐波的频率是 \(n\times100\) Hz，被乘上 \(|G|=\dfrac{1}{\sqrt{1+n^2}}\)，相位滞后 \(\arctan n\)：

高次谐波被压得比基波厉害得多，输出的方波被「磨圆」了。用时域的微分方程算这个输出要费不少功夫，在频域里几行就算完了。

反过来，同时记录一个系统的输入和输出，分别做 FFT，\(Y/X\) 就是 \(G(j\omega)\) 的估计。这是实测频率特性的基础（第10章 10.5 节）。

5. 计算机里的傅里叶：DFT 与 FFT

计算机拿到的不是连续的 \(x(t)\)，而是每隔 \(T_s=1/f_s\) 采一次、一共 \(N\) 个点的数据，记录长度 \(T=N/f_s\)。采样以后只能看到 \(f_s/2\)（奈奎斯特频率）以下的频率，更高的会混叠成假的低频（第8章 8.2 节），所以 A/D 之前要有抗混叠滤波。

对这 \(N\) 个点做离散傅里叶变换（DFT）：

第 \(k\) 个频点对应的频率是 \(f_k=k\,f_s/N\)。FFT（快速傅里叶变换）是算 DFT 的快速算法，算出来的数和按定义算完全一样，只是快得多：按定义算要约 \(N^2\) 次乘法，FFT 约 \(N\log_2N\) 次；\(N=1024\) 时是一百万次和一万次的差别。

怎么读 FFT 的结果

• 只看前一半：实信号的 \(X[k]\) 和 \(X[N-k]\) 共轭，后一半是负频率，只看 \(k=0\dots N/2\)（numpy 的 rfft 直接只算这一半）。
• 幅值换算：由第2节，DFT 的求和约等于 \(N\) 倍的级数系数，\(c_k\approx X[k]/N\)；再把负频率那一半加回来，正弦分量的幅值是
  $$A_k=\frac{2|X[k]|}{N}\quad(k\ne0,\ N/2),\qquad \text{直流}=\frac{|X[0]|}{N}$$
• 例：\(f_s=1000\) Hz，\(N=1000\)，\(\Delta f=1\) Hz。对第1节的信号做 FFT，第 0、50、120 个频点换算后正好是 0.200、1.000、0.300，其他频点都是 0。这是理想情况：每个频率在 1 s 里都是整数个周期。

在实验里改一改频率和点数：

泄漏与加窗

把 50 Hz 改成 50.5 Hz，1 s 里就有 50.5 个周期。不加窗时，峰值只剩 0.640（实际是 1.0，低了 36%），而且能量漏得到处都是，离它 10 Hz 远的 60 Hz 处还有 0.031。旁边如果有一个 0.02 的小信号，就被淹没了。

解决办法是加窗：做 FFT 之前，把数据乘上一个两头为 0、中间为 1 的窗函数，让首尾平滑地接起来，没有跳变。最常用的是汉宁窗 \(w[n]=0.5-0.5\cos\dfrac{2\pi n}{N-1}\)。加窗后，50.5 Hz 的峰值是 0.849，60 Hz 处的泄漏降到 0.0004，小了约 80 倍。

代价有两个：

• 主峰变宽：一个正好在频点上的正弦，不加窗只占 1 个频点，加汉宁窗后占 3 个（中间 1.0，两边各 0.5）；
• 幅值要重新换算：窗把数据「压小」了，公式里的 \(N\) 要换成窗函数的和 \(\sum w[n]\)。汉宁窗的和约为 \(N/2\)，相当于再乘 2。

分辨率：两个靠得很近的频率

两个频率要相差约 3 个 \(\Delta f\) 以上，才能可靠地看出两个峰。例：50 Hz（幅值 1.0）和 53 Hz（幅值 0.5），加汉宁窗：

• \(N=250\)（\(\Delta f=4\) Hz）：两个频率相差不到 1 个 \(\Delta f\)，只能看到 48 Hz 处一个 0.746 的峰，分不开；
• \(N=1000\)（\(\Delta f=1\) Hz）：相差 3 个 \(\Delta f\)，50 Hz 处 1.000、53 Hz 处 0.500，两个峰清清楚楚。

补零不能提高分辨率。把 1000 个点后面补上 3000 个 0 再做 FFT，频点间隔变成 0.25 Hz，曲线画得更细了，但 50 Hz 和 50.5 Hz 还是只有一个峰（50 Hz 处 1.087）。补零只是在原来的频点之间插值，数据里的信息还是那 1 s 的。

噪声

白噪声的能量分摊在所有频点上，\(N\) 越大，每个频点分到的越少，噪声底按 \(1/\sqrt N\) 下降；正弦分量的幅值不随 \(N\) 变。所以想从噪声里找出一个小的振动，就采长一点。第6节有一个例子：记录长度变成 16 倍，噪声底降到约 1/4。

反过来，正因为噪声的 FFT 幅值随 \(N\) 变，要描述噪声本身有多大，不能用幅值谱，要用功率谱密度（PSD，单位 V²/Hz），它不随 \(N\) 变。scipy.signal.welch 把长数据切成很多段，每段做 FFT 后再平均，得到的 PSD 曲线平滑得多。

6. 用 numpy 做 FFT

代码在 numpy 2.5.4、scipy 1.18.1 上运行过，注释里写的是实际输出，按页面顺序运行。

import numpy as np

fs, N = 1000, 1000                        # 采样频率 1 kHz，采 1000 点，记录长度 1 s
t = np.arange(N) / fs
x = 0.2 + 1.0 * np.sin(2 * np.pi * 50 * t) + 0.3 * np.sin(2 * np.pi * 120 * t)


def amp_spectrum(x, fs, window="hann"):
    """单边幅值谱：返回频率和每个频率上正弦分量的幅值（单位和 x 相同）"""
    n = len(x)
    w = np.hanning(n) if window == "hann" else np.ones(n)
    X = np.fft.rfft(x * w)                # 只算 0 ~ fs/2 这一半
    A = 2 * np.abs(X) / w.sum()           # 除以 N（加窗时除以窗的和），再乘 2 把负频率的一半加回来
    A[0] /= 2                             # 直流没有负频率，不乘 2
    if n % 2 == 0:
        A[-1] /= 2                        # fs/2 处也不乘 2
    return np.fft.rfftfreq(n, 1 / fs), A


f, A = amp_spectrum(x, fs, window=None)
print(f[1] - f[0])                        # 1.0：频率分辨率 fs/N
for k in (0, 50, 120):
    print(f"{f[k]:.0f} Hz: {A[k]:.3f}")   # 0 Hz: 0.200 / 50 Hz: 1.000 / 120 Hz: 0.300

泄漏和加窗：

x2 = np.sin(2 * np.pi * 50.5 * t)         # 1 s 里有 50.5 个周期，不是整数
for win in (None, "hann"):
    f, A = amp_spectrum(x2, fs, window=win)
    k = A.argmax()
    print(f"{win}: 峰在 {f[k]:.0f} Hz，幅值 {A[k]:.3f}；60 Hz 处 {A[60]:.4f}")
# None: 峰在 50 Hz，幅值 0.640；60 Hz 处 0.0307
# hann: 峰在 51 Hz，幅值 0.849；60 Hz 处 0.0004

要准确读幅值，换平顶窗：

from scipy.signal.windows import flattop

w = flattop(N)
A = 2 * np.abs(np.fft.rfft(x2 * w)) / w.sum()
print(f"{A.max():.3f}")                   # 0.999

噪声里找小振动：直流 3.0、噪声标准差 0.1，里面藏着一个幅值只有 0.02 的 180 Hz 振动：

rng = np.random.default_rng(0)
for n in (1000, 16000):
    tn = np.arange(n) / fs
    y = 3.0 + 0.02 * np.sin(2 * np.pi * 180 * tn) + 0.1 * rng.standard_normal(n)
    f, A = amp_spectrum(y - y.mean(), fs)  # 先去掉均值
    k = np.argmin(abs(f - 180))
    print(f"N = {n}: 180 Hz 处 {A[k]:.4f}，噪声底（中位数）{np.median(A):.4f}")
# N = 1000: 180 Hz 处 0.0140，噪声底（中位数）0.0063
# N = 16000: 180 Hz 处 0.0213，噪声底（中位数）0.0016

\(N=1000\) 时，振动只比噪声底高 2 倍，读数也被噪声带偏了（0.0140，实际是 0.02）；\(N=16000\) 时高出 13 倍，读数也准了。

分析实测数据的步骤

1. 确认采样频率：用时间戳算实际的 \(f_s\)，它不准，所有频率都会跟着错；
2. 定记录长度：要分开相差 \(\delta f\) 的两个频率，记录至少要 \(3/\delta f\) 长；
3. 去均值（有缓慢漂移就用 scipy.signal.detrend 去趋势）；
4. 加汉宁窗，做 rfft，按上面的公式换算幅值；
5. 看峰：读频率和幅值；要准确的幅值换平顶窗；要比较噪声大小用 PSD；
6. 在单片机上做 FFT 时，常用的库（比如 ARM 的 CMSIS-DSP）一般要求 \(N\) 是 2 的幂；numpy 没有这个限制。

本章小结

练习