5.3 谱分析与 EMD：从信号里挖出周期

5.3 谱分析与 EMD：从信号里挖出周期

本节摘要：周期长度不固定时，简单的 periodogram 难以稳定识别。Welch 估计通过分段平均提高稳健性，经验模态分解（EMD）则把信号自适应地拆为多个本征模态函数（IMF），每个 IMF 都有自己的瞬时频率。本节用一段合成 GDP 周期数据演示两种方法。

本节目标

阅读完本节，你应当能够：

1. 解释为什么简单 periodogram 在周期性识别中不够稳，以及 Welch 方法如何改进。
2. 写出 EMD 的工作流，理解"筛分"（sifting）过程的物理直觉。
3. 在 Python 中跑一次 EMD，把 GDP 长波分解为多个 IMF。

本节是"识别工具"专项：周期图太脆、Welch 用分辨率换稳健性、EMD 自适应分解为多个 IMF。读完后你能从一段不规则长波数据中挖出每个独立的振荡成分。

一、为什么简单 periodogram 容易失灵

4.2 节介绍过周期图（periodogram）——它把序列拆成不同频率的复指数波，统计每个频率的能量。问题在哪？

• 频率分辨率受窗口限制：长度为 N 的序列，频率分辨率约为 1/N。N=240（月度 20 年）时，分辨率是 1/240 ≈ 0.004，对应 240 个月 = 20 年。要分辨 5 年（60 月）和 7 年（84 月）周期，分辨率至少要 1/(84-60) ≈ 0.04。N=240 时勉强够。
• 噪声敏感：原始 periodogram 是直接对所有数据做 FFT，每个频率点的方差很大——峰值可能只是噪声放大，不是真实信号。
• 不适合非平稳信号：周期图假设"信号在整段观测期平稳"；周期性经常带漂移（5 年周期可能逐渐拉长到 6 年），periodogram 在这种情况下不适用。

二、Welch 方法：用分段平均提高稳健性

Welch 估计把序列切成多段重叠的子序列，分别算 periodogram 再平均：

import numpy as np
from scipy.signal import welch, periodogram
import matplotlib.pyplot as plt

# 合成 30 年月度 GDP：长趋势 + 5 年周期 + 噪声
rng = np.random.default_rng(0)
n = 360  # 30 年
t = np.arange(n)
trend = 0.05 * t
cycle = 8 * np.sin(2 * np.pi * t / 60 + 0.3)
y = trend + cycle + rng.normal(0, 1, n)
s = pd.Series(y)

# 简单 periodogram
freqs_p, power_p = periodogram(s.values, fs=1.0)

# Welch：分段长度 120 月（10 年），50% 重叠
freqs_w, power_w = welch(s.values, fs=1.0, nperseg=120, noverlap=60)
​

Welch 的代价是频率分辨率降低（分段后频率点更稀疏），但每个点的方差更小、估计更稳。工程上 Welch 几乎总是优于简单 periodogram。

💡 关键直觉：Welch 是"用分辨率换稳健性"。Nperseg 选 1/4 到 1/2 序列长度是常见经验。

三、EMD：自适应分解

经验模态分解（Empirical Mode Decomposition, EMD）是一种数据自适应的分解方法。它把任意信号拆成多个本征模态函数（Intrinsic Mode Function, IMF），每个 IMF 都满足两个条件：

1. 极值点数量与过零点数量之差 ≤ 1。
2. 上下包络的均值处处为零。

直觉上，IMF 是"信号里最简单的振荡"——有自己的频率、振幅，但不要求全局平稳。

EMD 的工作流

1. 找出信号的所有局部极大值和极小值。
2. 用三次样条插值得到上下包络。
3. 计算包络均值。
4. 原始信号减包络均值 = 候选 IMF。
5. 重复 1-4 直到候选 IMF 满足上述两个条件。
6. 把这个 IMF 从原信号中减掉，剩余部分继续分解，直到剩余信号无明显振荡。

from PyEMD import EMD
import numpy as np

# 信号：长趋势 + 5 年周期 + 1 年周期 + 噪声
t = np.arange(360)
cycle5 = 8 * np.sin(2 * np.pi * t / 60)
cycle1 = 3 * np.sin(2 * np.pi * t / 12)
y = 0.05 * t + cycle5 + cycle1 + rng.normal(0, 1, 360)

emd = EMD()
imfs = emd(y, t)   # 返回 (n_imfs, n_samples) 数组
print(f"分解出 {imfs.shape[0]} 个 IMF + 1 个残差")
​

输出通常是 4–6 个 IMF + 1 个残差：

• IMF 1：最高频（噪声 + 短周期）
• IMF 2：中频（1 年周期）
• IMF 3：低频（5 年周期）
• 残差：单调趋势

每个 IMF 都可以单独做频谱分析或建模。

四、EMD 的局限与变体

• 模态混叠：相邻 IMF 频率太接近时会出现"两个 IMF 共享一个频率"的情况。**EEMD（Ensemble EMD）**通过在原信号上加白噪声多次分解再平均来缓解。
• 边界效应：EMD 用样条插值估计包络，序列两端往往不准。CEEMDAN等改进方法用不同噪声注入策略。
• 理论保证弱：EMD 的数学收敛性没严格证明，工程上靠经验使用。

⚠️ 常见坑：不要把 EMD 当成"无监督的 SARIMA"。EMD 输出的 IMF 数量由数据决定，每次跑结果可能略有不同；要保证结果可复现，可以固定随机种子或用 EEMD。

五、一个完整案例：GDP 长波分解

下面用一段合成 GDP 月度数据演示 EMD + Welch 的完整工作流：

import pandas as pd
import numpy as np
from PyEMD import EMD
from scipy.signal import welch
import matplotlib.pyplot as plt

# 1) 合成 30 年月度 GDP 增速：长趋势 + 5 年周期 + 季度噪声
n = 360
t = np.arange(n)
trend = 5 + 0.02 * t
cycle = 1.5 * np.sin(2 * np.pi * t / 60 + 0.5)
quarterly = 0.3 * np.sin(2 * np.pi * t / 3)
y = trend + cycle + quarterly + np.random.default_rng(0).normal(0, 0.2, n)

# 2) EMD 分解
emd = EMD()
imfs = emd(y, t)

# 3) 对每个 IMF 做 Welch 谱估计
fig, axes = plt.subplots(len(imfs) + 1, 1, figsize=(10, 8))
axes[0].plot(t, y, lw=0.6); axes[0].set_title("原始")
for i, imf in enumerate(imfs):
    axes[i+1].plot(t, imf, lw=0.6)
    freqs, power = welch(imf, fs=1.0, nperseg=120)
    peak = freqs[np.argmax(power[1:]) + 1]
    period = 1 / peak if peak > 0 else np.inf
    axes[i+1].set_title(f"IMF{i+1}（主周期 ≈ {period:.1f} 月）")
plt.tight_layout()
​

输出会展示：

• IMF 1：高频噪声（主周期几月）
• IMF 2：季度噪声（主周期 ≈ 3 月）
• IMF 3：年度周期（主周期 ≈ 12 月）
• IMF 4：5 年长波（主周期 ≈ 60 月）
• 残差：单调趋势

每个 IMF 都可以进一步分析：长波 IMF 4 对应"经济周期"，年度 IMF 3 对应"季节性"（如果数据本身没这么强季节性，可能要靠业务判断）。

六、回到侦探视角

周期性识别就像从一段混合录音里分辨出"小提琴""大提琴""定音鼓"各自的声音：

• periodogram 是"一次性 FFT"——能听出大致频率，但分辨率有限。
• Welch 是"分段听"——更稳但精度降低。
• EMD 是"自适应滤波器"——自动找到每个乐器的频率，但有模态混叠风险。

实战中常常三种方法一起用：先用 Welch 看大致频段，再对每个频段用 EMD 拆细节，最后用业务知识确认"这个 IMF 是不是真的经济周期"。

本节要点回顾

• 简单 periodogram 局限：频率分辨率受窗口限制、噪声敏感、不适合非平稳信号。
• Welch 改进：分段平均提高稳健性，工程上几乎总是优于简单 periodogram。
• EMD 优势：自适应分解，不需要预设周期长度，能处理非平稳信号。
• EMD 局限：模态混叠、边界效应、理论保证弱；EEMD / CEEMDAN 是改进。
• 工程组合拳：Welch 看大致频段 + EMD 拆细节 + 业务知识确认。
• 下一章线索：第六章进入"建模与评估"——把分解工具和预测模型串起来，引入 STL、Prophet、ARIMA、SARIMA 等综合方法。

图：EMD 分解示意