简介这份资源面向具备一定Python基础、希望上手非线性信号处理的开发者与算法学习者提供变分模态分解VMD的可运行代码与配套测试用例可用于声音分析、故障诊断、金融时序等场景的信号分解实践。压缩包共2个文件包含1个py脚本与1个txt说明整体约1KB体量轻巧py文件承载VMD分解函数与测试逻辑txt则给出依赖库的安装指引便于快速搭建numpy、scipy环境。目前已有5452人学习下载热度较高。通过研读代码读者可理解模态个数K、迭代次数与尺度因子α等参数如何影响分解效果掌握将非平稳信号拆解为多个简谐模态、并检验模态正交性与频谱特性的完整流程也能借鉴其测试思路把VMD作为预处理环节接入机器学习或后端数据分析管道提升时序特征提取与预测精度。1. 变分模态分解到底解决了什么信号难题如果你做过旋转机械故障诊断、脑电节律提取或者电力谐波分析大概率被同一个问题折磨过一段非平稳信号里混着好几个中心频率不同的分量想拆开看用傅里叶变换会糊成一片用经验模态分解EMD又经常模态混叠端点还发散。变分模态分解VMD就是冲着这个痛点来的。它把「分解」这件事从递归筛分改写成变分求解先假设信号由 K 个有限带宽的模态叠加而成再让每个模态围绕自己的中心频率最紧凑同时所有模态加起来能重构原信号。这套思路的好处是抗噪、抗混叠模态数和中心频率都能提前约束结果比 EMD 稳定得多。这篇笔记不讲空理论直接给能跑的 python 代码、测试用例以及我调参时踩过的坑适合已经会写 python、想把这套方法落到自己数据上的工程师。2. VMD 的数学骨架与 python 实现选型2.1 从维纳滤波到 ADMMVMD 为什么能拆得干净VMD 的核心是把每个模态定义成一个调幅调频信号写成 $u_k(t) A_k(t)\cos(\phi_k(t))$然后对每个模态做希尔伯特变换得到解析信号再乘以 $e^{-j\omega_k t}$ 把频谱搬到基带。搬完之后用 L2 范数平方去度量这个基带信号的带宽带宽越小说明这个模态越「纯」。于是整个分解变成一个带约束的优化问题所有模态之和等于原信号。约束不好直接解就引入二次惩罚项和拉格朗日乘子转成无约束问题再用交替方向乘子法ADMM迭代。每一轮迭代里模态、中心频率、拉格朗日乘子轮流更新模态在频域有闭式解中心频率是模态频谱的重心。这就是为什么 VMD 比 EMD 稳它不是靠极值点筛出来的而是靠频域迭代收敛出来的。理解这一点对调参很关键。惩罚因子 α 控制的是带宽约束的强弱α 越大每个模态被逼得越窄容易把有用分量切碎α 越小模态越宽容易混叠。模态数 K 直接决定拆成几份K 给小了会欠分解两个分量挤在一个模态里K 给大了会过分解同一个物理分量被劈成两半中心频率还会靠得很近。这两个参数是 VMD 的命门后面会专门讲怎么定。2.2 用 vmdpy 跑通第一段信号最小可复现脚本python 里现成的 VMD 实现不多常见做法是用vmdpy这个包它把 ADMM 迭代封装好了接口干净。先装依赖再跑一段合成信号验证。pip install vmdpy numpy matplotlibimport numpy as np import matplotlib.pyplot as plt from vmdpy import VMD # 构造合成信号三个不同频率分量 噪声 fs 1000 # 采样率 Hz t np.arange(0, 1, 1/fs) # 1 秒时间轴 f1, f2, f3 30, 80, 150 # 三个分量中心频率 sig (np.cos(2*np.pi*f1*t) 0.6*np.cos(2*np.pi*f2*t) 0.4*np.cos(2*np.pi*f3*t)) sig 0.1 * np.random.randn(len(t)) # 加一点高斯噪声 # VMD 参数 alpha 2000 # 带宽约束中等强度 tau 0 # 噪声容限0 表示无松弛 K 3 # 模态数已知三个分量 DC 0 # 不含直流分量 init 1 # 中心频率均匀初始化 tol 1e-7 # 收敛容差 u, u_hat, omega VMD(sig, alpha, tau, K, DC, init, tol) # u: 各模态时域波形形状 (K, N) # u_hat: 各模态频域谱 # omega: 每次迭代的中心频率最后一行是收敛结果 print(收敛后中心频率(Hz):, omega[-1] * fs) plt.figure(figsize(10, 6)) for i in range(K): plt.subplot(K, 1, i1) plt.plot(t, u[i]) plt.title(fMode {i1}) plt.tight_layout() plt.show()这段代码的逻辑是先造一个频率成分已知的信号方便对照分解结果对不对。VMD返回三个东西u是时域模态u_hat是频域omega记录每轮迭代的中心频率。最后一行omega[-1] * fs把归一化频率还原成 Hz正常应该接近 30、80、150。参数上alpha2000是个常用起点K3是因为我们事先知道分量数init1表示中心频率均匀铺开初始化比全零初始化收敛更稳。跑完如果三个模态的频谱峰值分别落在 30、80、150 附近说明这套流程通了。2.3 参数怎么定K 和 alpha 的实操判据新手最容易卡在 K 和 alpha 上。我的经验是先定 K 再调 alpha。定 K 有个便宜办法对信号做 FFT看主峰个数主峰几个 K 就取几个再往上加一两个试。更严谨的做法是看分解后各模态的中心频率有没有靠得太近如果两个模态中心频率差不到主频的 10%基本就是过分解了。另一个判据是重构误差把所有模态加起来和原信号比误差应该很小如果 K 加大后误差没明显下降说明加多了。alpha 的取值和采样率、信号带宽有关。采样率高、分量窄的时候 alpha 要往大调常见范围是 1000 到 5000。判断标准是看模态频谱的带宽如果模态谱拖得很宽、和邻居重叠alpha 调大如果模态谱被切得只剩一根尖峰、时域波形失真alpha 调小。我一般会写个小循环扫几个 alpha画重构误差和模态带宽的曲线挑拐点。for a in [500, 1000, 2000, 4000, 8000]: u_tmp, _, om_tmp VMD(sig, a, 0, 3, 0, 1, 1e-7) recon u_tmp.sum(axis0) err np.linalg.norm(sig - recon) / np.linalg.norm(sig) print(falpha{a}, 重构相对误差{err:.4f}, 中心频率{om_tmp[-1]*fs})这段扫描代码能帮你快速看出 alpha 对结果的影响。注意重构误差不是越小越好过小的 alpha 会让模态变宽、误差也小但分解没意义。要结合中心频率是否稳定一起看。3. 测试用例怎么写从合成信号到真实数据的验证3.1 合成信号测试用例已知答案的回归验证给 VMD 写测试用例第一层是合成信号因为答案已知能当回归测试用。核心断言有三条分解出的模态数等于 K、中心频率和设定值偏差在容差内、重构误差低于阈值。用 pytest 组织方便每次改参数后重跑。import numpy as np import pytest from vmdpy import VMD def make_signal(fs1000, dur1.0, freqs(30, 80, 150), noise0.1): t np.arange(0, dur, 1/fs) sig sum(np.cos(2*np.pi*f*t) for f in freqs) rng np.random.default_rng(42) # 固定随机种子保证可复现 sig sig noise * rng.standard_normal(len(t)) return t, sig def test_vmd_recovers_known_frequencies(): fs 1000 freqs (30, 80, 150) t, sig make_signal(fsfs, freqsfreqs) u, _, omega VMD(sig, 2000, 0, 3, 0, 1, 1e-7) est np.sort(omega[-1] * fs) # 中心频率偏差不超过 5 Hz assert np.allclose(est, sorted(freqs), atol5), f频率偏差过大: {est} def test_vmd_reconstruction_error(): t, sig make_signal() u, _, _ VMD(sig, 2000, 0, 3, 0, 1, 1e-7) recon u.sum(axis0) err np.linalg.norm(sig - recon) / np.linalg.norm(sig) assert err 0.05, f重构误差过大: {err}这里两个用例分别盯频率恢复和重构精度。default_rng(42)固定种子是关键否则每次噪声不同断言阈值不好定。atol5是频率容差采样率 1000 Hz 下 5 Hz 已经够松如果实际偏差超过这个数说明 K 或 alpha 选错了。重构误差阈值 0.05 是经验值噪声大时可以放宽到 0.1。3.2 边界与异常用例短信号、单频、强噪合成信号跑通不代表代码健壮。真实数据里经常遇到信号很短、只有一个频率成分、信噪比很低的情况这些都要单独写用例。短信号的问题是 ADMM 迭代次数不够中心频率还没收敛就停了单频信号如果 K 设成 3会强行拆出三个模态其中两个是噪声这时候要断言「多余模态的能量占比很低」。def test_vmd_short_signal(): # 只有 100 个采样点验证不报错且能返回 K 个模态 t, sig make_signal(fs1000, dur0.1) u, _, _ VMD(sig, 2000, 0, 3, 0, 1, 1e-7) assert u.shape[0] 3 def test_vmd_single_tone_extra_modes_low_energy(): # 单频信号K3多出的模态能量应很低 fs 1000 t np.arange(0, 1, 1/fs) sig np.cos(2*np.pi*50*t) 0.05*np.random.randn(len(t)) u, _, _ VMD(sig, 2000, 0, 3, 0, 1, 1e-7) energy np.sum(u**2, axis1) energy_ratio np.sort(energy) / energy.sum() # 最小的两个模态能量占比之和应小于 0.2 assert energy_ratio[:2].sum() 0.2, f多余模态能量过高: {energy_ratio}短信号用例只断言形状因为 100 点下频率分辨率只有 10 Hz硬卡频率会误报。单频用例用能量占比判断过分解比看频率更稳。这两个用例能挡住大部分「参数乱设也能跑但结果没意义」的情况。3.3 真实数据接入以轴承振动信号为例真实数据接入时第一件事是确认采样率和单位。轴承振动信号常见采样率是 12 kHz 或 25.6 kHz直接丢进 VMD 之前最好先做归一化否则 alpha 的量纲对不上。我一般先把信号减均值、除标准差再送进 VMD。另外真实信号往往有趋势项低频漂移会干扰第一个模态可以先高通滤波或者去趋势。from scipy.signal import detrend def preprocess(x): x detrend(x) # 去线性趋势 x (x - x.mean()) / x.std() # 标准化 return x # 假设 raw 是从 csv 读入的一维振动信号 # raw np.loadtxt(bearing.csv, delimiter,)[:, 1] # sig preprocess(raw) # u, _, omega VMD(sig, 2000, 0, 5, 0, 1, 1e-7)预处理这步别省。我见过直接把原始振动信号丢进去的结果第一个模态全是趋势项后面几个模态挤在一起调半天参数都没用。标准化之后 alpha 的取值也更好迁移同一套参数换个数据集不至于完全失效。4. 避坑与排查VMD 落地时最容易翻车的五件事4.1 模态数 K 设大了中心频率挤在一起现象分解出 5 个模态但打印中心频率发现有两个只差几 Hz时域波形看着像同一个分量被劈开。原因K 超过信号实际分量数ADMM 会把一个宽带模态拆成两个窄的。解决先做 FFT 数主峰K 从主峰数开始试或者画中心频率随 K 变化的曲线找频率不再明显分开的拐点。我一般会跑 K2 到 8看中心频率分布出现两个模态频率差小于主频 10% 就停。4.2 alpha 太小导致模态混叠现象模态时域波形里能看到明显的拍频频谱上两个峰连在一起。原因alpha 太小带宽约束太弱ADMM 允许模态变宽去拟合更多成分。解决把 alpha 往上调常见从 2000 起每次翻倍试。注意 alpha 太大又会让模态谱变成尖峰、时域失真所以要两边夹。判断标准是模态频谱的 3dB 带宽应该和该分量的物理带宽量级一致。4.3 信号没预处理低频趋势吃掉第一个模态现象第一个模态几乎是一条缓慢变化的曲线后面模态才是有用的振荡分量。原因原始信号有直流或线性趋势VMD 把它当成一个低频模态。解决进 VMD 前做去趋势和标准化。如果趋势是非线性的先用高通滤波或者经验方法去掉。这一步不做后面参数怎么调都是白费。4.4 短信号迭代不收敛结果每次不一样现象信号长度只有几百点跑两次 VMD 结果中心频率差很多。原因ADMM 迭代次数上限到了还没收敛而初始化对短信号影响大。解决短信号要么加长数据窗要么固定init参数并放宽tol同时把最大迭代次数调大。如果还是不稳考虑对信号做镜像延拓再分解取中间段结果。4.5 把 VMD 当万能分解器什么信号都上现象对冲击成分为主的信号做 VMD分解结果没有物理意义。原因VMD 假设模态是窄带准平稳的冲击信号是宽带瞬态不符合假设。解决先判断信号类型窄带多分量用 VMD宽带瞬态用其他方法或者先做包络分析。工具没有好坏只有合不合适硬套只会得到一堆没法解释的模态。5. 进阶技巧用中心频率稳定性反推最优 K调 K 最烦的是没有客观标准。我后来固定用一个技巧对同一个信号K 从 2 扫到 10每个 K 跑多次改变初始化或加微小扰动记录每次收敛后的中心频率看哪些频率在所有 K 下都稳定出现。真正属于信号的频率会反复出现过分解产生的虚假频率会随 K 跳来跳去。这个方法比单看一次分解结果靠谱得多。def stable_frequencies(sig, fs, k_rangerange(2, 11), repeats5): from collections import defaultdict freq_count defaultdict(int) for K in k_range: for r in range(repeats): # 每次加极小扰动模拟初始化差异 noise 1e-6 * np.random.randn(len(sig)) _, _, omega VMD(sig noise, 2000, 0, K, 0, 1, 1e-7) for f in omega[-1] * fs: # 把频率归到 5 Hz 的桶里统计 bucket round(f / 5) * 5 freq_count[bucket] 1 # 出现次数多的频率就是稳定分量 stable sorted([f for f, c in freq_count.items() if c repeats], keylambda x: -freq_count[x]) return stable, freq_count # stable, counts stable_frequencies(sig, fs) # print(稳定中心频率:, stable)这段代码把频率按 5 Hz 分桶统计出现次数达到重复次数的频率才算稳定。repeats5是每个 K 跑 5 次k_range覆盖 2 到 10。跑完看stable列表前几个就是信号的真实分量频率K 取这个个数就行。这个方法对噪声有一定容忍度因为噪声产生的虚假频率不会稳定复现。参数上分桶宽度 5 Hz 要按采样率调采样率低的时候桶要宽一点否则同一个频率会被分到相邻桶里。repeats越大越稳但越慢一般 5 次够用。如果信号很长可以只取一段做这个扫描定好 K 再对全长信号分解。我自己现在的习惯是拿到新信号先跑这个稳定性扫描把 K 定下来再单独调 alpha最后写两个 pytest 用例把参数锁住。这样换数据、换同事接手结果都能复现。VMD 不难难的是参数有依据、结果可验证把这两件事做扎实它就能稳定产出可解释的模态。希望帮到你。本文还有配套的精品资源点击获取