vmd

Variational mode decomposition

collapse all in page

Syntax

imf = vmd(x)
[imf,residual] = vmd(x)
[imf,residual,info] = vmd(x)
[___] = vmd(x,Name,Value)
vmd(___)

Description

example

imf = vmd(x) returns the variational mode decomposition of x. Use vmd to decompose and simplify complicated signals into a finite number of intrinsic mode functions (IMFs) required to perform Hilbert spectral analysis.

[imf,residual] = vmd(x) also returns the residual signal residual corresponding to the variational mode decomposition of x.

example

[imf,residual,info] = vmd(x) returns additional information info on IMFs and the residual signal for diagnostic purposes.

example

[___] = vmd(x,Name,Value) performs the variational mode decomposition with additional options specified by one or more Name,Value pair arguments.

example

vmd(___) plots the original signal, IMFs, and the residual signal as subplots in the same figure.

Examples

collapse all

Open Live Script

Create a signal, sampled at 4 kHz, that resembles dialing all the keys of a digital telephone. Save the signal as a MATLAB® timetable.

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);

Plot the variational mode decomposition of the timetable.

vmd(S)

Open Live Script

Generate a multicomponent signal consisting of three sinusoids of frequencies 2 Hz, 10 Hz, and 30 Hz. The sinusoids are sampled at 1 kHz for 2 seconds. Embed the signal in white Gaussian noise of variance 0.01².

fs = 1e3;
t = 1:1/fs:2-1/fs;
x = cos(2*pi*2*t) + 2*cos(2*pi*10*t) + 4*cos(2*pi*30*t) + 0.01*randn(1,length(t));

Compute the IMFs of the noisy signal and visualize them in a 3-D plot.

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')

Use the computed IMFs to plot the Hilbert spectrum of the multicomponent signal. Restrict the frequency range to [0, 40] Hz.

hht(imf,fs,'FrequencyLimits',[0,40])

Open Live Script

Generate a piecewise composite signal consisting of a quadratic trend, a chirp, and a cosine with a sharp transition between two constant frequencies at t = 0.5.

x(t)=6t2 +cos(4πt+10πt2)+{cos(60πt),t≤0.5,cos(100πt-10π),t>0.5.

The signal is sampled at 1 kHz for 1 second. Plot each individual component and the composite signal. 

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')

Perform variational mode decomposition to compute four intrinsic mode functions. The four distinct components of the signal are recovered.

[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

Reconstruct the signal by adding the mode functions and the residual. Plot and compare the original and reconstructed signals.

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')

Calculate the norm of the difference between the original and the reconstructed signals.

norm(x-sig',Inf)
ans = 0
Open Live Script

The signals labeled in this example are from the MIT-BIH Arrhythmia Database [3] (Signal Processing Toolbox). The signal in the database was sampled at 360 Hz.

Load the MIT database signals corresponding to record 200 and plot the signal.

load mit200
Fs = 360;
plot(tm,ecgsig)
ylabel('Time (s)')
xlabel('Signal')

The ECG signal contains spikes driven by the rhythm of the heartbeat and an oscillating low frequency pattern. The distinct spokes of the ECG create important higher order harmonics.

Calculate nine intrinsic mode functions of the windowed signal. Visualize the IMFs. 

[imf,residual] = vmd(ecgsig,'NumIMF',9);
t = tiledlayout(3,3,'TileSpacing','compact','Padding','compact');
for n = 1:9
    ax(n) = nexttile(t);
    plot(tm,imf(:,n)')
    xlim([tm(1) tm(end)])
    txt = ['IMF',num2str(n)];
    title(txt)
    xlabel('Time (s)')
end
title(t,'Variational Mode Decomposition')

The first mode contains the most noise, and the second mode oscillates at the frequency of the heartbeat. Construct a clean ECG signal by summing all but the first and last VMD modes, thus discarding the low frequency baseline oscillation and most of the high frequency noise.

cleanECG = sum(imf(:,2:8),2);
figure
plot(tm,ecgsig,tm,cleanECG)
legend('Original ECG','Clean ECG')
ylabel('Time (s)')
xlabel('Signal')

Input Arguments

collapse all

Uniformly sampled time-domain signal, specified as either a vector or a timetable. If x is a timetable, then it must contain increasing finite row times.. The timetable must contain only one numeric data vector with finite load values.

Note

If a timetable has missing or duplicate time points, you can fix it using the tips in Clean Timetable with Missing, Duplicate, or Nonuniform Times.

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'NumIMF',10

Mode convergence absolute tolerance, specified as the comma-separated pair consisting of 'AbsoluteTolerance' and a positive real scalar. AbsoluteTolerance is one of the stopping criteria for optimization, that is, optimization stops when the average squared absolute improvement toward convergence of IMFs, in two consecutive iterations, is less than AbsoluteTolerance.

Mode convergence relative tolerance, specified as the comma-separated pair consisting of 'RelativeTolerance' and a positive real scalar. RelativeTolerance is one of the stopping criteria for optimization, that is, optimization stops when the average relative improvement toward convergence of IMFs, in two consecutive iterations, is less than RelativeTolerance.

Note

The optimization process stops when 'AbsoluteTolerance' and 'RelativeTolerance' are jointly satisfied.

Maximum number of optimization iterations, specified as the comma-separated pair consisting of 'MaxIterations' and a positive scalar integer. MaxIterations is one of the stopping criteria for optimization, that is, optimization stops when the number of iterations is greater than MaxIterations.

MaxIterations can be specified using only positive whole numbers.

Number of IMFs extracted, specified as the comma-separated pair consisting of 'NumIMF' and a positive scalar integer.

Initial central IMF frequencies, specified as the comma-separated pair consisting of 'CentralFrequencies' and a vector of length NumIMFs. Vector values must be within [0, 0.5] cycles/sample, which indicates that the true frequencies are within [0, fs/2], where fs is the sample rate.

Initial IMFs, specified as the comma-separated pair consisting of 'InitialIMFs' and a real matrix. The rows correspond to time samples and columns correspond to modes.

Penalty factor, specified as the comma-separated pair consisting of 'PenaltyFactor' and a positive real scalar. This argument determines the reconstruction fidelity. Use smaller values of penalty factor to obtain stricter data fidelity.

Initial Lagrange multiplier, specified as the comma-separated pair consisting of 'InitialLM' and a complex vector. The range of the initial Lagrange multiplier in the frequency domain is [0, 0.5] cycles/sample. The multiplier enforces the reconstruction constraint. The length of the multiplier depends on the input size.

Update rate for the Lagrange multiplier in each iteration, specified as the comma-separated pair consisting of 'LMUpdateRate' and a positive real scalar. A higher rate results in faster convergence, but increases the chance of the optimization process getting stuck in a local optimum.

Method to initialize the central frequencies, specified as the comma-separated pair consisting of 'InitializeMethod' and either 'peaks', 'random', or 'grid'.

Specify InitializeMethod as:

  • 'peaks' to initialize the central frequencies as the peak locations of the signal in the frequency domain (default).

  • 'random' to initialize the central frequencies as random numbers distributed uniformly in the interval [0,0.5] cycles/sample.

  • 'grid' to initialize the central frequencies as a uniformly sampled grid in the interval [0,0.5] cycles/sample.

Toggle progress display in the command window, specified as the comma-separated pair consisting of 'Display' and either 'true' (or 1) or 'false' (or 0). If you specify 'true', the function displays the average absolute and relative improvement of modes and central frequencies every 20 iterations, and show the final stopping information.

Specify Display as 1 to show the table or 0 to hide the table.

Output Arguments

collapse all

Intrinsic mode functions, returned as a matrix or timetable. Each imf is an amplitude and frequency modulated signal with positive and slowly varying envelopes. Each mode has an instantaneous frequency that is nondecreasing, varies slowly, and is concentrated around a central value. Use imf to apply Hilbert-Huang transform to perform spectral analysis on the signal.

imf is returned as:

  • A matrix where each column is an imf, when x is a vector

  • A timetable, when x is a single data column timetable

Residual signal, returned as a column vector or a single data column timetable. residual represents the portion of the original signal x not decomposed by vmd.

residual is returned as:

  • A column vector, when x is a vector.

  • A single data column timetable, when x is a single data column timetable.

Additional information for diagnostics, returned as a structure with these fields:

  • ExitFlag – Termination flag. A value of 0 indicates the algorithm stopped when it reached the maximum number of iterations. A value of 1 indicates the algorithm stopped when it met the absolute and relative tolerances.

  • CentralFrequencies – Central frequencies of the IMFs.

  • NumIterations – Total number of iterations.

  • AbsoluteImprovement – Average squared absolute improvement toward convergence of the IMFs between the final two iterations.

  • RelativeImprovement – Average relative improvement toward convergence of the IMFs between the final two iterations.

  • LagrangeMultiplier – Frequency-domain Lagrange multiplier at the last iteration.

More About

collapse all

Intrinsic Mode Functions

The vmd function decomposes a signal x(t) into a small number K of narrowband intrinsic mode functions (IMFs):

x(t)=∑k=1Kuk(t).

The IMFs have these characteristics:

  1. Each mode uk is an amplitude and frequency modulated signal of the form

    uk(t)=Ak(t)cos(ϕk(t)),

    where ϕk(t) is the phase of the mode and Ak(t) is its envelope.

  2. The modes have positive and slowly varying envelopes.

  3. Each mode has an instantaneous frequency ϕ'k(t) that is nondecreasing, varies slowly, and is concentrated around a central value fk.

The variational mode decomposition method simultaneously calculates all the mode waveforms and their central frequencies. The process consists of finding a set of uk(t) and fk(t) that minimize the constrained variational problem.

Optimization

To calculate uk and fk, the procedure finds an optimum of the augmented Lagrangian

L(uk(t),fk,λ(t))=α∑k=1K‖ddt[(δ(t)+jπt)∗uk(t)]e−j2πfkt‖22+‖x(t)−∑k=1Kuk(t) ‖22+〈λ(t),x(t)−∑k=1Kuk(t)〉(i)(ii)(iii),

where the inner product 〈p(t),q(t)〉=∫−∞∞p∗(t) q(t) dt and the 2-norm ‖p(t)‖22=〈p(t),p(t)〉. The regularization term (i) includes these steps:

  1. Use the Hilbert transform to calculate the analytic signal associated with each mode, where * denotes convolution. This results in each mode having a purely positive spectrum.

  2. Demodulate the analytic signal to baseband by multiplying it with a complex exponential.

  3. Estimate the bandwidth by calculating the squared 2-norm of the gradient of the demodulated analytic signal.

Terms (ii) and (iii) enforce the constraint x(t)=∑k=1Kuk(t) by imposing a quadratic penalty and incorporating a Lagrange multiplier. The PenaltyFactor α measures the relative importance of (i) compared to (ii) and (iii).

The algorithm solves the optimization problem using the alternating direction method of multipliers described in [1] (Signal Processing Toolbox).

Algorithms

The vmd function calculates the IMFs in the frequency domain, reconstructing X(f) = DFT{x(t)} in terms of Uk(f) = DFT{uk(t)}. To remove edge effects, the algorithm extends the signal by mirroring half its length on either side.

The Lagrange multiplier introduced in Optimization (Signal Processing Toolbox) has the Fourier transform Ʌ(f). The length of the Lagrange multiplier vector is the length of the extended signal.

Unless otherwise specified in 'InitialIMFs', the IMFs are initialized at zero. Initialize 'CentralFrequencies' using one of the methods specified in 'InitializeMethod'. vmd iteratively updates the modes until one of these conditions is met:

  • ∑k‖ukn+1(t)−ukn(t)‖22/‖ukn(t)‖22<εr and ∑k‖ukn+1(t)−ukn(t)‖22<εa are jointly satisfied, where εr and εa are specified using 'RelativeTolerance' and 'AbsoluteTolerance', respectively.

  • The algorithm exceeds the maximum number of iterations specified in 'MaxIterations'.

For the (n + 1) -th iteration, the algorithm performs these steps:

  1. Iterate over the K modes of the signal (specified using 'NumIMF') to compute:

    1. The frequency-domain waveforms for each mode using

      Ukn+1(f)=X(f)−∑i<kUkn+1(f)−∑i>kUkn(f)+Λn2(f)1+2α{2π(f−fkn)}2,

      where Ukn+1(f) is the Fourier transform of the kth mode calculated in the (n + 1)-th iteration.

    2. The kth central frequency fkn+1 using

      fkn+1=∫0∞|Ukn+1(f)|2 f df∫0∞|Ukn+1(f)|2 df≈∑f|Ukn+1(f)|2 ∑|Ukn+1(f)|2 .

  2. Update the Lagrange multiplier using Λn+1(f)=Λn(f)+τ(X(f)−∑kUkn+1(f)), where τ is the update rate of the Lagrange multiplier, specified using 'LMUpdateRate'.

References

[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.

Extended Capabilities

See Also

|

Introduced in R2020a

Wavelet Toolbox Documentation

Support