fs = 4e3;
t = 0:1/fs:0.5-1/fs;

ver = [697 770 852 941];
hor = [1209 1336 1477];

tones = [];

for k = 1:length(ver)
    for l = 1:length(hor)
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones = [tones;tone;zeros(size(tone))];
    end
end

% To hear, type soundsc(tones,fs)

S = timetable(tones,'SampleRate',fs);

绘制该时间表的变分模态分解。

vmd(S)

MATLAB figure

多分量信号的 VMD

打开实时脚本

生成一个多分量信号，该信号由频率为 2 Hz、10 Hz 和 30 Hz 的三个正弦波组成。这些正弦波以 1 kHz 的采样率进行采样，采样时间为 1 秒。将信号嵌入方差为 0.01² 的高斯白噪声中。

Fs = 1e3;
t = (1/Fs:1/Fs:1)';
x = cos(2*pi*[2 10 30].*t)*[1 2 4]' + 0.01*randn(size(t));

计算含噪信号的 IMF，并在三维绘图中进行可视化。

imf = vmd(x);
[p,q] = ndgrid(t,1:size(imf,2));
plot3(p,q,imf)
grid on
xlabel("Time Values")
ylabel("Mode Number")
zlabel("Mode Amplitude")

Figure contains an axes object. The axes object with xlabel Time Values, ylabel Mode Number contains 5 objects of type line.

使用计算出的 IMF 绘制多分量信号的希尔伯特频谱。将频率范围限制为 [0, 40] Hz。

hht(imf,Fs,FrequencyLimits=[0,40])

Figure contains an axes object. The axes object with title Hilbert Spectrum, xlabel Time (s), ylabel Frequency (Hz) contains 5 objects of type patch.

分段信号的 VMD

打开实时脚本

生成一个分段合成信号，该信号由二次趋势、啁啾以及在 t = 0.5 处两个恒定频率之间具有急剧过渡的余弦组成。

$x (t) = 6 t^{2} + \cos (4 π t + 10 π t^{2}) + {\begin{array}{ll} \cos (60 π t), & t \leq 0.5, \\ \cos (100 π t - 10 π), & t > 0.5 . \end{array}$

该信号在 1 kHz 下采样 1 秒。绘制每个单独分量和合成信号。

fs = 1e3;
t = 0:1/fs:1-1/fs;

x = 6*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...
    [cos(60*pi*(t(t<=0.5))) cos(100*pi*(t(t>0.5)-10*pi))];

tiledlayout('flow')
nexttile
plot(t,[zeros(1,length(t)/2) cos(100*pi*(t(length(t)/2+1:end))-10*pi)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,[cos(60*pi*(t(1:length(t)/2))) zeros(1,length(t)/2)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,cos(4*pi*t+10*pi*t.^2))
xlabel('Time (s)')
ylabel('Chirp')


nexttile
plot(t,6*t.^2)
xlabel('Time (s)')
ylabel('Quadratic trend')

nexttile(5,[1 2])
plot(t,x)
xlabel('Time (s)')
ylabel('Signal')

执行变分模态分解以计算四个固有模态函数。信号的四个不同分量被还原。

[imf,res] = vmd(x,'NumIMFs',4);

tiledlayout('flow')

for i = 1:4
    nexttile
    plot(t,imf(:,i))
    txt = ['IMF',num2str(i)];
    ylabel(txt)
    xlabel('Time (s)')
    grid on
end

通过将模态函数和残差相加来重新构造信号。绘制并比较原始信号和重新构造的信号。

sig = sum(imf,2) + res;

nexttile(5,[1 2])
plot(t,sig,'LineWidth',1.5)

hold on

plot(t,x,':','LineWidth',2)
xlabel('Time (s)')
ylabel('Signal')
hold off
legend('Reconstructed signal','Original signal', ...
       'Location','northwest')

计算原始信号和重新构造的信号之间差值的范数。

norm(x-sig',Inf)

ans = 
0

使用 VMD 从 ECG 信号中去除噪声

打开实时脚本

加载对应于 MIT-BIH 心律失常数据库的第 200 条记录的信号 [3]。该信号以 360 Hz 的采样率进行采样。绘制信号。

load mit200

plot(tm,ecgsig)
xlabel("Time (s)")
xlim([tm(1) tm(end)])

Figure contains an axes object. The axes object with xlabel Time (s) contains an object of type line.

ECG 信号包含由心跳节律驱动的尖峰以及振荡的低频模式。ECG 的不同尖峰会产生重要的高阶谐波。

计算加窗信号的九个固有模态函数。可视化 IMF。

[imf,residual] = vmd(ecgsig,NumIMFs=9);

t = tiledlayout(3,3, ...
    TileSpacing="compact",Padding="compact");
for n = 1:9
    nexttile(t)
    plot(tm,imf(:,n))
    xlim([tm(1) tm(end)])
    title("IMF "+n)
    xlabel("Time (s)")
end
title(t,"Variational Mode Decomposition")

通过对除第一个和最后一个之外的所有 VMD 模态求和来构造干净的 ECG 信号（第一个模态包含大部分高频噪声，最后一个模态包含低频基线振荡）。

cleanECG = sum(imf(:,2:8),2);

figure
plot(tm,ecgsig,":",tm,cleanECG)
legend("Original ECG","Clean ECG")
xlabel("Time (s)")
xlim([tm(1) tm(end)])

Figure contains an axes object. The axes object with xlabel Time (s) contains 2 objects of type line. These objects represent Original ECG, Clean ECG.

输入参数

全部折叠

`x` — 均匀采样的时域信号
向量 | 时间表

均匀采样的时域信号，指定为向量或时间表。

如果 x 是时间表，则它必须包含递增的有限行时间。时间表必须仅包含一个具有有限加载值的数值数据向量。

注意
如果时间表有缺失或重复的时间点，您可以使用清理包含缺失、重复或不均匀时间的时间表中的提示进行修复。
仅支持对 RowTimes 使用 duration 或 datetime 向量的时间表。

名称-值参数

全部折叠

将可选参量对组指定为 Name1=Value1,...,NameN=ValueN，其中 Name 是参量名称，Value 是对应的值。名称-值参量必须出现在其他参量之后，但对各个参量对组的顺序没有要求。

示例: imf = vmd(x,MaxIterations=200,NumIMFs=10)

`AbsoluteTolerance` — 模态收敛绝对容差
`5e-6` (默认) | 正实数标量

模态收敛绝对容差，指定为正实数标量。AbsoluteTolerance 是优化的停止条件之一，即当 IMF 在两次连续迭代中朝向收敛的平均平方绝对改进小于 AbsoluteTolerance 时，优化停止。

`RelativeTolerance` — 模态收敛相对容差
`AbsoluteTolerance*1e3` (默认) | 正实数标量

模态收敛相对容差，指定为正实数标量。RelativeTolerance 是优化的停止条件之一，即当 IMF 在两次连续迭代中朝向收敛的平均相对改进小于 RelativeTolerance 时，优化停止。

注意

当同时满足 AbsoluteTolerance 和 RelativeTolerance 时，优化过程停止。

`MaxIterations` — 最大优化迭代次数
`500` (默认) | 正标量整数

最大优化迭代次数，指定为正标量整数。MaxIterations 是优化的停止条件之一，即当迭代次数大于 MaxIterations 时，优化停止。

MaxIterations 只能使用正整数指定。

`NumIMFs` — 提取的 IMF 数
`5` (默认) | 正标量整数

提取的 IMF 数，指定为正标量整数。

`CentralFrequencies` — 初始中心 IMF 频率
向量

初始中心 IMF 频率，指定为长度是 NumIMFs 的向量。向量值必须在 [0, 0.5] 周期/采样内，这表明真实频率在 [0, f_s/2] 内，其中 f_s 是采样率。

`InitialIMFs` — 初始 IMF
零矩阵 (默认) | 实矩阵

初始 IMF，指定为实矩阵。行对应于时间采样，列对应于模态。

`PenaltyFactor` — 罚分因子
`1000` (默认) | 正实数标量

罚分因子，指定为正实数标量。此参量确定重新构造保真度。使用的罚分因子值越小，获得的数据保真度越严格。

`InitialLM` — 初始拉格朗日乘数
由零组成的复数向量 (默认) | 复数向量

初始拉格朗日乘数，指定为复数向量。初始拉格朗日乘数在频域中的范围为 [0, 0.5] 周期/采样。该乘数强制执行重新构造约束。乘数的长度取决于输入大小。

`LMUpdateRate` — 拉格朗日乘数的更新速率
`0.01` (默认) | 实数标量

每次迭代中拉格朗日乘数的更新速率，指定为正实数标量。速率越高，收敛越快，但会增大优化过程陷入局部最优的机会。

`InitializeMethod` — 初始化中心频率的方法
`"peaks"` (默认) | `"random"` | `"grid"`

初始化中心频率的方法，指定为 "peaks"、"random" 或 "grid"。

"peaks" - 将中心频率初始化为信号在频域中的峰值位置（默认值）。如果 NumIMFs 大于频域峰值数，vmd 会随机初始化其余中心频率。
"random" - 将中心频率初始化为在区间 [0,0.5] 周期/采样内均匀分布的随机数。
"grid" - 将中心频率初始化为在区间 [0,0.5] 周期/采样内均匀采样的网格。

`Display` — 切换命令行窗口中的进度显示
`false` 或 0 (默认) | `true` 或 1

切换命令行窗口中的进度显示，指定为 "true"（或 1）或 "false"（或 0）。如果指定 "true"，该函数每 20 次迭代就输出模态和中心频率的平均绝对和相对改进，并显示最终停止信息。

将 Display 指定为 1 以显示该表，或指定为 0 以隐藏该表。

输出参量

全部折叠

`imf` — 固有模态函数
矩阵 | 时间表

固有模态函数，以矩阵或时间表形式返回。每个 imf 是一个振幅和频率调制信号，具有正的且缓慢变化的包络。每个模态具有非递减、缓慢变化并集中在中心值附近的瞬时频率。使用 imf 应用希尔伯特-黄变换以对信号执行频谱分析。

imf 以如下形式返回：

矩阵，其中每列是 imf（当 x 是向量时）。
时间表（当 x 是单数据列时间表时）。

`residual` — 残差信号
列向量 | 单数据列时间表

残差信号，以列向量或单数据列时间表形式返回。residual 表示原始信号 x 中未被 vmd 分解的部分。

residual 以如下形式返回：

列向量（当 x 是向量时）。
单数据列时间表（当 x 是单数据列时间表时）。

`info` — 用于诊断的附加信息
结构体

用于诊断的附加信息，以具有以下字段的结构体形式返回：

ExitFlag - 终止标志。值 0 表示算法在达到最大迭代次数时停止。值 1 表示算法在满足绝对和相对容差时停止。
CentralFrequencies - IMF 的中心频率。
NumIterations - 迭代总数。
AbsoluteImprovement - 最后两次迭代之间 IMF 朝向收敛的平均平方绝对改进。
RelativeImprovement - 最后两次迭代之间 IMF 朝向收敛的平均相对改进。
LagrangeMultiplier - 最后一次迭代时的频域拉格朗日乘数。

详细信息

全部折叠

固有模态函数

vmd 函数将信号 x(t) 分解为少量 K 个窄带固有模态函数 (IMF)：

$x (t) = \sum_{k = 1}^{K} u_{k} (t) .$

这些 IMF 具有以下特征：

每个模态 u_k 是一个振幅和频率调制信号，形式为

$u_{k} (t) = A_{k} (t) \cos (ϕ_{k} (t)),$
其中 ϕ_k(t) 是该模态的相位，A_k(t) 是其包络。
这些模态具有正的且缓慢变化的包络。
每个模态具有瞬时频率 ϕ'_k(t)，该瞬时频率是非递减的、缓慢变化的，并且集中在中心值 f_k 附近。

变分模态分解方法同时计算所有模态波形及其中心频率。该过程包括找到一组 u_k(t) 和 f_k(t)，使约束变分问题最小化。

优化

为了计算 u_k 和 f_k，该过程找到增广拉格朗日的最优值

$\begin{matrix} L (u_{k} (t), f_{k}, λ (t)) & = & α \sum_{k = 1}^{K} {‖ \frac{d}{d t} [(δ (t) + \frac{j}{π t}) * u_{k} (t)] e^{- j 2 π f_{k} t} ‖}_{2}^{2} & + & {‖ x (t) - \sum_{k = 1}^{K} u_{k} (t) ‖}_{2}^{2} & + & 〈 λ (t), x (t) - \sum_{k = 1}^{K} u_{k} (t) 〉 \\ (i) & (i i) & (i i i) \end{matrix},$

其中内积 $〈 p (t), q (t) 〉 = \int_{- \infty}^{\infty} p^{*} (t) q (t) d t$ ，2-范数 ${‖ p (t) ‖}_{2}^{2} = 〈 p (t), p (t) 〉$ 。正则化项 (i) 包括以下步骤：

使用希尔伯特变换计算与每个模态相关联的解析信号，其中 * 表示卷积。这会导致每个模态具有纯正频谱。
通过将解析信号乘以复指数来将其解调到基带。
通过计算解调解析信号梯度的平方 2-范数来估计带宽。

项 (ii) 和 (iii) 通过施加二次惩罚项并引入拉格朗日乘数来强制执行约束 $x (t) = \sum_{k = 1}^{K} u_{k} (t)$ 。PenaltyFactor α 测量 (i) 与 (ii) 和 (iii) 相比的相对重要性。

该算法使用 [1] 中所述的交替方向乘数法求解优化问题。

算法

vmd 函数在频域中计算 IMF，以 U_k(f) = DFT{u_k(t)} 形式重新构造 X(f) = DFT{x(t)}。为了去除边沿效果，该算法在信号两端分别镜像扩展其长度的一半。

优化中引入的拉格朗日乘数具有傅里叶变换 Ʌ(f)。拉格朗日乘数向量的长度是扩展后的信号的长度。

除非在 InitialIMFs 中另有指定，否则 IMF 在零处进行初始化。使用在 InitializeMethod 中指定的方法之一初始化 CentralFrequencies。vmd 以迭代方式更新模态，直到满足以下条件之一：

同时满足 $\sum_{k} {‖ u_{k}^{n + 1} (t) - u_{k}^{n} (t) ‖}_{2}^{2} / {‖ u_{k}^{n} (t) ‖}_{2}^{2} < ε_{r}$ 和 $\sum_{k} {‖ u_{k}^{n + 1} (t) - u_{k}^{n} (t) ‖}_{2}^{2} < ε_{a}$ ，其中 ε_r 和 ε_a 分别使用 RelativeTolerance 和 AbsoluteTolerance 指定。
算法超出在 MaxIterations 中指定的最大迭代次数。

对于第 (n + 1) 次迭代，该算法执行以下步骤：

循环访问信号的 K 个模态（使用 NumIMFs 指定）以计算：
1. 每个模态的频域波形，方法是使用
  
  $U_{k}^{n + 1} (f) = \frac{X (f) - \sum_{i < k} U_{k}^{n + 1} (f) - \sum_{i > k} U_{k}^{n} (f) + \frac{Λ^{n}}{2} (f)}{1 + 2 α {2 π (f - f_{k}^{n})}^{2}},$
  其中 $U_{k}^{n + 1} (f)$ 是在第 (n + 1) 次迭代中计算的第 k 个模态的傅里叶变换。
2. 第 k 个中心频率 $f_{k}^{n + 1}$ ，方法是使用
  
  $f_{k}^{n + 1} = \frac{\int_{0}^{\infty} {| U_{k}^{n + 1} (f) |}^{2} f d f}{\int_{0}^{\infty} {| U_{k}^{n + 1} (f) |}^{2} d f} \approx \frac{\sum f {| U_{k}^{n + 1} (f) |}^{2}}{\sum {| U_{k}^{n + 1} (f) |}^{2}} .$
使用 $Λ^{n + 1} (f) = Λ^{n} (f) + τ (X (f) - \sum_{k} U_{k}^{n + 1} (f)),$ 更新拉格朗日乘数，其中 τ 是拉格朗日乘数的更新速率，使用 LMUpdateRate 指定。

参考

[1] Boyd, Stephen, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. "Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers." Foundations and Trends^® in Machine Learning. Vol 3, Number 1, 2011, pp. 1–122.

[2] Dragomiretskiy, Konstantin, and Dominique Zosso. "Variational Mode Decomposition." IEEE^® Transactions on Signal Processing. Vol. 62, Number 3, 2014, pp. 531–534.

[3] Moody, George B., and Roger G. Mark. "The impact of the MIT-BIH Arrhythmia Database." IEEE Engineering in Medicine and Biology Magazine. Vol. 20, No. 3, May–June 2001, pp. 45–50.

扩展功能

全部展开

C/C++ 代码生成
使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

用法说明和限制：

如果 x 是时间表，则可选输入 NumIMFs 必须为常量。

GPU 代码生成
使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

请参阅“C/C++ 代码生成”部分中的用法说明和限制。同样的用法说明和限制也适用于 GPU 代码生成。

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

vmd 函数支持基于线程的环境，但存在以下用法说明和限制：

不支持不带输出参量的语法。

有关详细信息，请参阅在基于线程的环境中运行 MATLAB 函数。

GPU 数组
通过使用 Parallel Computing Toolbox™ 在图形处理单元 (GPU) 上运行来加快代码执行。

vmd 函数完全支持 GPU 数组。要在 GPU 上运行该函数，请将输入数据指定为 gpuArray (Parallel Computing Toolbox)。有关详细信息，请参阅在 GPU 上运行 MATLAB 函数 (Parallel Computing Toolbox)。

版本历史记录

在 R2020a 中推出

全部展开

R2026b: 使用基于线程的环境

您可以在基于线程的环境中使用 vmd。有关详细信息，请参阅在基于线程的环境中运行 MATLAB 函数。

R2026b: 使用 `gpuArray` 对象

vmd 函数支持 gpuArray 对象。您必须有 Parallel Computing Toolbox™ 才能使用此功能。

另请参阅

vmd

语法

描述

示例

拨号音信号的变分模态分解

多分量信号的 VMD

分段信号的 VMD

使用 VMD 从 ECG 信号中去除噪声

输入参数

x — 均匀采样的时域信号 向量 | 时间表

名称-值参数

AbsoluteTolerance — 模态收敛绝对容差 5e-6 (默认) | 正实数标量

RelativeTolerance — 模态收敛相对容差 AbsoluteTolerance*1e3 (默认) | 正实数标量

MaxIterations — 最大优化迭代次数 500 (默认) | 正标量整数

NumIMFs — 提取的 IMF 数 5 (默认) | 正标量整数

CentralFrequencies — 初始中心 IMF 频率 向量

InitialIMFs — 初始 IMF 零矩阵 (默认) | 实矩阵

PenaltyFactor — 罚分因子 1000 (默认) | 正实数标量

InitialLM — 初始拉格朗日乘数 由零组成的复数向量 (默认) | 复数向量

LMUpdateRate — 拉格朗日乘数的更新速率 0.01 (默认) | 实数标量

InitializeMethod — 初始化中心频率的方法 "peaks" (默认) | "random" | "grid"

Display — 切换命令行窗口中的进度显示 false 或 0 (默认) | true 或 1

输出参量

imf — 固有模态函数 矩阵 | 时间表

residual — 残差信号 列向量 | 单数据列时间表

info — 用于诊断的附加信息 结构体

详细信息

固有模态函数

优化

算法

参考

扩展功能

C/C++ 代码生成 使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成 使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

基于线程的环境 使用 MATLAB® backgroundPool 在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool 加快代码运行速度。

GPU 数组 通过使用 Parallel Computing Toolbox™ 在图形处理单元 (GPU) 上运行来加快代码执行。

版本历史记录

R2026b: 使用基于线程的环境

R2026b: 使用 gpuArray 对象

另请参阅

App

函数

主题

`x` — 均匀采样的时域信号
向量 | 时间表

`AbsoluteTolerance` — 模态收敛绝对容差
`5e-6` (默认) | 正实数标量

`RelativeTolerance` — 模态收敛相对容差
`AbsoluteTolerance*1e3` (默认) | 正实数标量

`MaxIterations` — 最大优化迭代次数
`500` (默认) | 正标量整数

`NumIMFs` — 提取的 IMF 数
`5` (默认) | 正标量整数

`CentralFrequencies` — 初始中心 IMF 频率
向量

`InitialIMFs` — 初始 IMF
零矩阵 (默认) | 实矩阵

`PenaltyFactor` — 罚分因子
`1000` (默认) | 正实数标量

`InitialLM` — 初始拉格朗日乘数
由零组成的复数向量 (默认) | 复数向量

`LMUpdateRate` — 拉格朗日乘数的更新速率
`0.01` (默认) | 实数标量

`InitializeMethod` — 初始化中心频率的方法
`"peaks"` (默认) | `"random"` | `"grid"`

`Display` — 切换命令行窗口中的进度显示
`false` 或 0 (默认) | `true` 或 1

`imf` — 固有模态函数
矩阵 | 时间表

`residual` — 残差信号
列向量 | 单数据列时间表

`info` — 用于诊断的附加信息
结构体

C/C++ 代码生成
使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成
使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

GPU 数组
通过使用 Parallel Computing Toolbox™ 在图形处理单元 (GPU) 上运行来加快代码执行。

R2026b: 使用 `gpuArray` 对象