%% Noise Removal from ECG Signal using VMD
% The signals labeled in this example are from the MIT-BIH Arrhythmia Database 3. The signal in the database was sampled at 360 Hz.
% Load the MIT database signals corresponding to record 200.

load mit200
Fs = 360;

% Window the signal using a Gaussian key window of the same length as the ECG signal.

w = gausswin(length(ecgsig),5);
sigWin = w.*ecgsig;

% Plot the windowed and original signal.

plot(tm,sigWin,tm,ecgsig,':')
legend('Windowed Signal','Original Signal')

% The windowed signal contains spikes driven by the rhythm of the heartbeat and an oscillating low frequency pattern. The distinct spokes of the ECG signal create important higher order harmonics.
%
% Calculate nine intrinsic mode functions of the windowed signal. Visualize the IMFs.

imf = vmd(sigWin,'NumIMF',9);

tiledlayout(3,3)

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

%
% 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,sigWin,tm,cleanECG)
legend('Original ECG','Clean ECG')

% Copyright 2019 The MathWorks, Inc.