sktime 中的 vmdpy 库：变分模态分解（VMD）的 Python 实现与实战指南

sktime 中的 vmdpy 库变分模态分解VMD的 Python 实现与实战指南【免费下载链接】sktimeA unified framework for machine learning with time series项目地址: https://gitcode.com/GitHub_Trending/sk/sktime导读本文以 sktime 仓库内置的vmdpy库位于 sktime/libs/vmdpy为讲解主体深入介绍变分模态分解Variational Mode Decomposition, VMD方法的原理、安装方式、核心 API 参数含义以及一套可直接运行的三分量合成信号分解示例。读完本文你将掌握如何调用VMD函数对一维时间序列进行模态分解、正确理解每个超参数的作用并了解如何基于该底层库使用 sktime 的高层封装VmdTransformer将 VMD 无缝嵌入分解—预测—重构的完整流水线。一、vmdpy 是什么VMD 方法的 Python 实现vmdpy是一个用于对信号执行**变分模态分解VMD**的 Python 库其算法源自 Dragomiretskiy 与 Zosso 于 2014 年发表的经典论文Variational Mode DecompositionIEEE Transactions on Signal Processing, vol. 62, no. 3, pp. 531–544, 2014。该包是原版VMD MATLAB 工具箱由 Dominique Zosso 维护的 Python 移植版本。在 sktime 中vmdpy是 sktime/libs 目录下的“随 sktime 分发并受官方维护的库”之一。根据 sktime/libs/README.md 的说明vmdpy自 2023 年 8 月起成为 sktime 官方维护的 fork从 sktime/libs/vmdpy/init.py 的模块文档可以看出其演进历史首版创建于 2019 年 2 月0.1 版本发布于 2019 年 4 月 9 日0.2 版本发布于 2020 年 8 月 11 日2023 年 8 月迁移进入 sktime作为官方 fork 持续维护。该库既可以被sktime的估计器内部使用例如 sktime/transformations/vmd.py 中的VmdTransformer也可以脱离 sktime 直接作为独立函数库使用。其实现与测试的核心文件如下文件作用sktime/libs/vmdpy/vmdpy.pyVMD 算法的完整 Python 实现VMD函数与辅助函数_safe_averagesktime/libs/vmdpy/init.py包入口导出VMD包含许可证与版本历史sktime/libs/vmdpy/tests/test_readme_example.py对 README 示例的回归测试验证示例代码可运行sktime/transformations/vmd.py面向 sktime 用户的VmdTransformer高层封装sktime/transformations/tests/test_vmd.pyVmdTransformer的流水线与序列长度测试二、安装随 sktime 一起分发由于vmdpy已经作为内部库随 sktime 分发用户无需单独安装vmdpy只需安装 sktime 即可使用可通过pip或conda安装pip install sktime或conda install sktime安装完成后即可直接导入from sktime.libs.vmdpy import VMD需要注意的适用前提VMD函数的核心依赖仅为numpy实现中使用math与numpy而 README 示例中的可视化部分需要额外安装matplotlib。测试代码 test_readme_example.py 中正是用_check_soft_dependencies(matplotlib, severitynone)对绘图部分做了软依赖检查缺省matplotlib时仅跳过绘图。三、核心 APIVMD函数签名与参数详解VMD的函数签名为见 sktime/libs/vmdpy/vmdpy.pydef VMD(f, alpha, tau, K, DC, init, tol):3.1 参数含义参数类型含义说明farray_like待分解的时域信号一维函数内部会执行f np.array(f)做类型转换alphafloat数据保真约束的平衡参数带宽约束控制各模态的带宽值越大模态带宽越窄对噪声更敏感taufloat对偶上升dual ascent的时间步长取 0 表示噪声松弛noise-slack不强制严格保真Kint需要恢复的模态数量即分解出的本征模态函数IMF个数DCbool是否强制第一个模态保持在 0 频DC 分量为真时第一个模态的中心频率被锁定为 0initintomega中心频率的初始化方式0 所有 omega 从 0 开始1 均匀分布初始化2 随机初始化tolfloat收敛判据容差典型取值约1e-6示例中取1e-73.2 返回值VMD返回三个对象返回值含义u分解得到的各模态集合K 行 × 信号长度列即 K 个 IMFu_hat各模态的频谱复数频谱omega估计得到的各模态中心频率在 README 示例中u的形状为(K, len(f))因此绘制模态时使用plt.plot(u.T)将每个模态作为一条曲线。四、算法原理从源码看 VMD 的迭代流程理解VMD的内部实现有助于正确选择参数。结合 sktime/libs/vmdpy/vmdpy.py 的源码其计算流程可分为以下几个阶段4.1 信号镜像延拓与频谱构造为避免边界效应算法先将原信号做镜像延拓mirror extension再进入频域处理vmdpy.pyfs 1.0 / len(f) fMirr np.array(np.flip(f[0 : math.ceil(T / 2)])) fMirr np.append(fMirr, f) fMirr np.append(fMirr, np.flip(f[math.ceil(T / 2) :]))随后对延拓后的信号做 FFT 并进行fftshift得到单边频谱f_hat_plus将负频部分置零这是 VMD 处理实信号的标准做法。4.2 中心频率 omega 的初始化init参数在这里生效vmdpy.pyinit 1omega_plus[0, i] (0.5 / K) * i即在(0, 0.5)频段内均匀分布init 2在(fs, 0.5)区间内按对数均匀采样并排序引入随机性init 0全部置 0若DC为真则第一个 omega 被强制设为 0。4.3 交替方向乘子法ADMM主循环核心求解采用交替方向乘子法ADMM主循环持续迭代直到收敛判据uDiff tol或达到最大迭代次数Niter 500vmdpy.py更新模态频谱每个模态通过残差信号的Wiener 滤波1.0 Alpha[k] * (freqs - omega_plus[n, k]) ** 2作为分母更新更新中心频率以模态频谱的功率谱密度作为权重对频率做加权平均通过辅助函数_safe_average对偶上升lambda_hat按tau步长更新实现约束的逐步满足收敛检查计算相邻两次迭代模态频谱之差的能量小于tol即终止。4.4 后处理去镜像、重构频谱迭代收敛后算法丢弃因镜像延拓产生的多余部分u u[:, T//4 : 3*T//4]并基于共轭对称关系重构完整双边频谱最后通过逆 FFT 得到时域模态u。值得一提的细节是辅助函数_safe_averagevmdpy.py当所有权重之和为 0 时退化为普通算术平均避免除零异常——这种边界处理保证了算法在极端初始化下的数值稳定性。五、完整示例三分量合成信号的 VMD 分解README 提供了一个可完整运行的示例构造一个由 3 个不同频率分量叠加并加入高斯噪声的合成信号然后分解回 3 个模态并可视化。该示例同时被 test_readme_example.py 作为回归测试固化保证了示例代码与实现的一致性。#%% Simple example: generate signal with 3 components noise import numpy as np import matplotlib.pyplot as plt from sktime.libs.vmdpy import VMD # Time Domain 0 to T T 1000 fs 1 / T t np.arange(1, T 1) / T freqs 2 * np.pi * (t - 0.5 - fs) / (fs) # center frequencies of components f_1 2 f_2 24 f_3 288 # modes v_1 np.cos(2 * np.pi * f_1 * t) v_2 1 / 4 * (np.cos(2 * np.pi * f_2 * t)) v_3 1 / 16 * (np.cos(2 * np.pi * f_3 * t)) f v_1 v_2 v_3 0.1 * np.random.randn(v_1.size) # some sample parameters for VMD alpha 2000 # moderate bandwidth constraint tau 0.0 # noise-tolerance (no strict fidelity enforcement) K 3 # 3 modes DC 0 # no DC part imposed init 1 # initialize omegas uniformly tol 1e-7 # Run VMD u, u_hat, omega VMD(f, alpha, tau, K, DC, init, tol) # Visualize decomposed modes plt.figure() plt.subplot(2, 1, 1) plt.plot(f) plt.title(Original signal) plt.xlabel(time (s)) plt.subplot(2, 1, 2) plt.plot(u.T) plt.title(Decomposed modes) plt.xlabel(time (s)) plt.legend([Mode %d % m_i for m_i in range(u.shape[0])]) plt.tight_layout()5.1 示例参数选型解读alpha 2000中等带宽约束示例注释为“moderate bandwidth constraint”。alpha越大各模态频带越窄、分离越激进但过大会导致模态过度分割或振荡tau 0.0噪声松弛不强制执行保真约束允许分解结果保留一定的噪声容差适合含噪信号K 3模态数与信号真实的 3 个频率分量2、24、288 Hz匹配。模态数需要与信号结构匹配这是 VMD 使用中最关键的先验参数DC 0不强制直流分量信号中无恒定的 DC 偏置因此不需要把第一个模态锁定在 0 频init 1均匀初始化中心频率在频段内均匀散布有助于稳定收敛tol 1e-7严格收敛容差比论文推荐的典型值1e-6更严格换取更高的分解精度。5.2 运行结果解读运行后u为(3, 1000)的数组对应 3 个分解模态omega记录各模态中心频率的迭代轨迹最终收敛值应接近真实分量频率 2、24、288。下图上部为原始信号下部为分解出的 3 个模态曲线。六、从底层函数到 sktime 封装VmdTransformerVMD底层函数解决的是“给定 K如何分解”的问题而 sktime 进一步将其封装为符合估计器接口的转换器VmdTransformer见 sktime/transformations/vmd.py解决了“K 未知时如何自动确定模态数”这一实际难题。6.1 关键特性K 的自动搜索在 VmdTransformer 中当KNone时__runVMDUntilCoefficientThreshold方法vmd.py会从K1开始逐步递增每次分解后计算能量损失系数energy_loss_coef np.linalg.norm((data - reconstruct), 2) ** 2 / np.linalg.norm(data, 2)当重构信号各模态求和与原始信号的相对能量误差低于阈值默认energy_loss_coefficient0.01或 K 达到上限kMax30时停止取满足条件的最小 K 作为分解模态数。6.2 核心参数对照VmdTransformer完整继承了VMD的参数并补充了高层控制项参数默认值说明KNone模态数为None时自动搜索kMax30自动搜索时模态数上限energy_loss_coefficient0.01自动搜索的收敛阈值能量损失系数alpha2000带宽约束对应底层alphatau0.0噪声容差/对偶上升步长对应底层tauDC0是否强制 DC 分量对应底层DCinit1omega 初始化方式对应底层inittol1e-7收敛容差对应底层tolreturned_decompu返回内容u返回模态、u_hat返回模态频谱绝对值、u_both返回列拼接的模态频谱6.3 在预测流水线中使用VMD 的典型应用模式是“分解—预测—重构”先用 VMD 把复杂序列分解为若干更易学习的 IMF逐个预测再求和还原。这一模式可直接用 sktime 的流水线语法实现vmd.pyfrom sktime.transformations.vmd import VmdTransformer from sktime.datasets import load_solar from sktime.forecasting.trend import TrendForecaster y load_solar() transformer VmdTransformer() modes transformer.fit_transform(y) pipe VmdTransformer() * TrendForecaster() pipe.fit(y, fh[1, 2, 3]) y_pred pipe.predict()测试 test_vmd.py 验证了这一流水线模式VmdTransformer与TrendForecaster组合成TransformedTargetForecaster后可以正常fit与predict。同时VmdTransformer支持inverse_transform各模态按行求和即可还原序列并具备多变量capability:multivariate: True能力标签定义见 vmd.py。另有测试test_vmd.py确保无论输入序列长度为 1000偶数还是 1001奇数分解结果长度都与输入一致——这对应了__runVMDUntilCoefficientThreshold中对奇数长度序列做末尾补值的处理逻辑。七、引用与参与贡献如果 vmdpy 在你的研究中发挥作用README 建议引用其方法论论文Carvalho, Moraes, Braga, Mendes.Evaluating five different adaptive decomposition methods for EEG signal seizure detection and classification.Biomedical Signal Processing and Control, Volume 62, 2020, 102073, ISSN 1746-8094, doi:10.1016/j.bspc.2020.102073.算法本身的理论出处为 Dragomiretskiy 与 Zosso2014的 VMD 原始论文。功能建议、问题反馈可通过 sktime 的 issue tracker 或讨论区提交可联系作者vrcarva对包的新功能或修复可以以 PR 形式提交到 sktime 仓库的sktime.libs.vmdpy模块。八、总结与适用边界vmdpy为 sktime 提供了高质量的 VMD 底层实现其价值体现在两个层面独立函数层VMD(f, alpha, tau, K, DC, init, tol)接口简洁仅依赖 numpy可脱离 sktime 单独使用适合信号处理场景中的快速模态分解估计器层VmdTransformer将其无缝融入 sktime 生态自动确定模态数、支持逆变换与流水线组合适合时间序列预测与特征工程场景。使用时的核心注意事项K模态数是影响分解质量的最关键先验alpha控制带宽约束的松紧二者需要结合信号自身的频谱结构与噪声水平进行调优当模态数不确定时优先使用VmdTransformer的自动搜索机制并通过energy_loss_coefficient控制信息损失的可接受程度。【免费下载链接】sktimeA unified framework for machine learning with time series项目地址: https://gitcode.com/GitHub_Trending/sk/sktime创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考