经验模态分解去趋势：用 EMD 把趋势与周期自适应剥离

去趋势是几乎所有时间序列分析的预处理：剥掉慢趋势再看波动、看周期、做配对。但怎么剥？最常见的做法是简单移动平均——窗口设多大？这是个隐藏的”超参猜猜乐”，窗口偏小留周期、偏大滞后趋势，而且真实市场里周期本身就在漂移，根本没有一个固定窗能同时回答”剥掉 21 日与 63 日两层周期”和”跟住非平稳趋势”两个问题。

经验模态分解（EMD）走一条完全不同的路：它不做滤波，而是让数据自己把每个尺度拆出来。本文从零实现 EMD 的 sifting 算法，演示它在多周期+非线性慢趋势上的自适应分离，并诚实展示它与移动平均的真实差距——不少、但代价可以接受，关键是它不挑窗，对非平稳趋势也稳健。

一、EMD 在做什么#

EMD 把任意信号拆成一组**本征模态函数（IMF）**加一条单调残差。IMF 的定义是 Huang 等人 1998 年给出的，两条：（a）极值点个数与过零点个数相差不超过 1；（b）上下包络关于时间轴局部对称。这两条保证了 IMF 是”窄带”的振荡分量，可以独立做 Hilbert 变换求瞬时频率。

Sifting 是从一个候选信号里反复减去上下包络均值，直到它满足 IMF 条件：

import numpy as np
from scipy.interpolate import CubicSpline

def envelope(x, kind):
    """上下包络：局部极值点三次样条，端点 hold 防止外推发散。"""
    if kind == "max":
        idx = np.where((x[1:-1] > x[:-2]) & (x[1:-1] > x[2:]))[0] + 1
    else:
        idx = np.where((x[1:-1] < x[:-2]) & (x[1:-1] < x[2:]))[0] + 1
    if len(idx) < 2:
        return None
    val = x[idx]
    cs = CubicSpline(idx, val, bc_type="natural")
    t = np.arange(len(x))
    out = np.zeros_like(x, dtype=float)
    inter = (t >= idx[0]) & (t <= idx[-1])
    out[inter] = cs(t[inter])
    out[t < idx[0]] = val[0]
    out[t > idx[-1]] = val[-1]
    return out

def emd(x, max_imf=6, max_sift=15, sd_thr=0.25):
    """逐层 sifting：每次提取最高频 IMF，更新残差，残差近单调时停止。"""
    residue = x.astype(float).copy()
    imfs = []
    while len(imfs) < max_imf and np.std(residue) > 1e-7 * np.std(x):
        h = residue.copy()
        for _ in range(max_sift):
            up = envelope(h, "max"); lo = envelope(h, "min")
            if up is None or lo is None: break
            h_new = h - 0.5 * (up + lo)
            sd = np.mean((h - h_new) ** 2) / (np.mean(h ** 2) + 1e-12)
            h = h_new
            if sd < sd_thr: break
        imfs.append(h.copy())
        residue = residue - h
        if np.all(np.diff(residue) >= 0) or np.all(np.diff(residue) <= 0): break
    return imfs, residue

几个关键细节：包络用三次样条而不是线性插值，Huang 原文坚持用三次样条因为包络的光滑度直接决定 IMF 的频率纯度；端点用 hold 而不是外推，因为 CubicSpline 的 bc_type="natural" 在数据外仍会延伸，但容易在两端发散；停止条件除了 sd_thr 还加了一条”残差近单调就停”——这是 sifting 后期最常见的过度筛分陷阱，残差一旦变成直线就再榨不出 IMF 了，强行提取只会把趋势再分成几个伪 IMF。

max_imf 要根据信号里的真实尺度数设置——少到 4 通常足够，多到 8 容易把趋势拆碎；max_sift 和 sd_thr 是孪生参数，sd_thr 越小每层 IMF 越纯净但耗时越长，且越容易把慢成分当 IMF 抽走。

二、合成实验：非平稳慢趋势 + 双周期 + 噪声#

为了让 EMD 的”自适应”价值可证伪，我设计了一个让固定窗 MA 难受的信号：

np.random.seed(20260829)
N = 1024
t = np.arange(N, dtype=float)
# 慢趋势：S 形 regime 漂移（前后半段方向相反）+ 轻微线性趋势
trend_true = -0.25 * np.tanh((t - 512) / 80.0) + 0.10 * (t / N)
# 双周期：21d 快 + 63d 中
cycle_med = 0.30 * np.sin(2 * np.pi * t / 63.0)
cycle_fast = 0.15 * np.sin(2 * np.pi * t / 21.0)
x = trend_true + cycle_med + cycle_fast + 0.04 * np.random.randn(N)

numpy 把 tanh 设计成 S 形是有意为之——窗口在 63 天的 MA 在断点附近要滞后整整一个窗口长度，这正是固定窗 MA 最大的失败模式。

跑 imfs, residue = emd(x, max_imf=6)，得到 6 个 IMF 加 1 个残差（图 1）。从 IMF1（最高频）到 IMF5（最清晰的一条 ~252 日等值正弦）频谱自上而下，残差（绿色）紧贴真实 S 形慢趋势。

几个值得标记的观察：

• IMF1 几乎是白噪声——它吸收了原始噪声的高频部分。EMD 在 sifting 第一层时，上下包络完全由噪声驱动，结果就是一个高频 IMF。这一层很薄，所以”去噪”基本是 EMD 的副作用，不是主要功能。
• IMF2 抓住 21 日周期，IMF4 抓住 63 日周期——频谱逐层下移正是 EMD 的核心卖点。固定窗 MA 给不出这种”按尺度自动对齐”的分解。
• 残差（绿色）是非线性慢趋势。这与傅里叶和小波都不同——后两者需要预先选基函数周期，EMD 完全没有这步。

三、EMD 残差对真实趋势的相关性 0.958#

图 2 把残差（EMD 认定的趋势）单独拉出来和真实趋势叠在一起：

def safe_corr(a, b):
    if np.std(a) < 1e-12 or np.std(b) < 1e-12:
        return float("nan")
    return np.corrcoef(a, b)[0, 1]
corr_trend = safe_corr(residue, trend_true)         # 0.9579
corr_fast  = max(safe_corr(im, cycle_fast) for im in imfs_keep)  # 0.684
corr_med   = max(safe_corr(im, cycle_med)  for im in imfs_keep)  # 0.657

corr(residue, true_trend) = 0.958——EMD 在没有任何先验信息（既不知道周期也不知道趋势形状）的情况下，用一个 6 层 sifting 就恢复了非线性慢趋势的 95.8% 相关结构。剩下的 4.2% 主要来自 IMF 之间的边界泄漏：21 日 IMF 和残差之间的能量有少量串扰，这是 sifting 的固有问题，不是 bug。

IMF 与注入周期的相关性（0.68 / 0.66）看起来不高，原因是 EMD 不保证 IMF 严格正交，21 日分量被 IMF2 抓住了一部分、又被 IMF3 残留下了一部分。这是 EMD 区别于”硬正交分解”（PCA、SSA）的地方——它追求频率局部化而非正交性，相关性自然不是 1。

四、固定窗 MA 的两难 + EMD 不挑窗#

公平起见，我做了完整窗口扫描，把 MA(5), MA(6), …, MA(120) 对真实趋势的 MSE 全部画出来（图 3 左），并把 EMD 的 MSE 画成绿色虚线：

ws = np.arange(5, 121)
mse_sweep = [np.mean((np.convolve(x, np.ones(w)/w, mode="same") - trend_true) ** 2) for w in ws]
best_w = ws[np.argmin(mse_sweep)]; best_mse = min(mse_sweep)
mse_emd = np.mean((residue - trend_true) ** 2)
mse_ma15  = np.mean((np.convolve(x, np.ones(15)/15, mode="same") - trend_true) ** 2)
mse_ma120 = np.mean((np.convolve(x, np.ones(120)/120, mode="same") - trend_true) ** 2)

数据是诚实的：

关键观察：

• MSE 呈 U 形——MA 的失败模式是结构性的：窗口小留周期，窗口大滞后趋势，且真实数据里两条边界都未知，扫窗口找最优其实是数据窥探（lookahead）。
• EMD 的 MSE 不如最优 MA——这是真的。在已知”信号含 63 日周期”这个隐含知识的前提下，MA63 是被精心调到位的；EMD 不知道这件事，自然无法对齐。但 EMD 不知道这件事本来就是它的设计——它从未被要求过要预知周期。
• EMD 在 100x MA15 和 1.2x 最优 MA 之间——也就是说”无窗”换来的是 ~100 倍的鲁棒性，付出的是 ~100 倍的精度（在已知周期的前提下）。对未知的真实市场，前者更值钱。

五、Hilbert-Huang 时频谱：EMD 真正的杀手锏#

EMD 单独做趋势提取并不是它最不可替代的部分。真正让它区别于傅里叶和小波的是 Hilbert-Huang 谱（HHT）——对每个 IMF 做 Hilbert 变换，提取瞬时频率与瞬时幅值，把”周期随时变化”这件事画出来：

from scipy.signal import hilbert
periods = np.logspace(np.log10(8), np.log10(400), 120)
H = np.zeros((len(periods), 200))
for im in imfs_keep:
    an = hilbert(im)
    phase = np.unwrap(np.angle(an))
    inst_p = 1.0 / np.diff(phase) / (2 * np.pi)
    amp = np.abs(an[1:])
    valid = (inst_p > 8) & (inst_p < 400)
    pe = np.clip(np.searchsorted(periods, inst_p) - 1, 0, len(periods) - 1)
    te = np.clip(((np.arange(1, N) / N) * 200).astype(int), 0, 199)
    for k in range(N - 1):
        if valid[k]:
            H[pe[k], te[k]] += amp[k] ** 2

Y 轴是对数周期（21 / 63 / 126 / 252 日四档标记），X 轴是交易日，颜色是能量 ∑|IMF|²。注意几个事：

• 能量集中在 60-100 日区间——这正好是注入的 63 日周期。HHT 自动把它定位到正确的频带，不需要预先指定。
• 强度随时变化——这是傅里叶根本做不出来的：傅里叶给一个固定频带，HHT 告诉你”63 日这个周期在第 200-300 天特别强、第 600 天开始衰减”。这等于免费给了一个周期强度时序信号，可以直接喂给 regime detector 或自适应择时。
• 21 日周期很弱——HHT 看不到 IMF1 的 21 日能量，因为 IMF1 主要是噪声，瞬时频率噪声大。HHT 对低信噪比 IMF 的瞬时频率估计本身就是个研究问题，本文用 valid = (inst_p > 8) & (inst_p < 400) 做了简单清洗。

六、EMD 不是银弹：三个必须知道的坑#

第一，mode mixing。当信号里有间歇性高频（如价格冲击），EMD 会把同一段时间内的高频能量分散到多个 IMF 里，让 IMF 失去”单一频率”的物理意义。Wu 和 Huang 2009 年的 EEMD（集合经验模态分解）通过加白噪声平均解决这个问题，但代价是计算量翻几个数量级。本文 numpy 从零实现的就是经典 EMD，没做抗 mode mixing 处理。

第二，端点效应。三次样条包络在两端只能 hold，极值不够时会偏离真实包络，结果就是 IMF 在两端失真。延拓方法（镜像、AR 预测、神经网络延拓）可以缓解但无法根除。本文合成数据 N=1024 端点占 0.2%，对整体指标影响可忽略；真实短序列要警惕。

第三，参数敏感性。max_imf、max_sift、sd_thr 三个参数都影响结果，没有统一默认值。本文的 max_imf=6, max_sift=15, sd_thr=0.25 是为这个特定数据调的；换数据要重新选。max_imf 设大了会把慢趋势拆碎成伪 IMF，残差相关会掉到 0.8 以下；设小了会漏掉中间尺度。

七、和上一篇小波的关系#

上一篇 多尺度波动率分解 用 MODWT(à trous) 也是把信号拆到多个尺度。两者最大的区别是 EMD 的尺度是数据自适应的、每个 IMF 的频率带宽由 sifting 决定；MODWT 的尺度是预定义的二进频带。换句话说，MODWT 知道”我要看 16-32 日这个频带”，EMD 不知道——它只知道”这个 IMF 大概在某个频率上振荡”。

对金融数据，两者经常互补：EMD 用来发现”哪些周期是真实的”（data-driven 周期识别），MODWT 用来在已知周期上做带通滤波（fixed-band energy estimation）。HHT 则是 EMD 独占的副产品。

八、结论#

EMD 在含 21d/63d 双周期 + S 形非平稳趋势的 1024 天合成数据上：6 层 sifting 抽出 6 个 IMF，残差对真实趋势相关 0.958；MA 窗口扫描 MSE 呈 U 形（最优 MA63 mse 0.00004，EMD mse 0.004），EMD 不必选窗即落在合理区间；Hilbert-Huang 谱把能量自适应定位到 60-100 日频带并显示其随时间的强度变化。

诚实地讲，EMD 不是去趋势的最优解——在已知周期时，固定窗 MA 的 MSE 仍然更低。EMD 的真正价值在于当周期未知且趋势非平稳时仍能给出一个像样的趋势估计，并且额外提供 Fourier/wavelet 给不出的瞬时频率信息。这两个特点决定了 EMD 在 regime detection、周期漂移建模、冲击事件分析这三个场景里仍无可替代。