morphological operation and filtering

3 views (last 30 days)
ah sho
ah sho on 5 Jan 2020
Edited: Max Murphy on 23 Feb 2020
Hello, I wrote this code trying to correct baseline in ECG signal but it didn't work, any help?
B0 = strel('line',0.2*fs*Tw,0);
Bc = strel('line',1.5*0.2*fs*Tw,0);
fb = imopen(f0, B0);
fb = imclose(fb, Bc);
fbc = f0 - fb;
figure(2)
plot(fbc);
title('baseline correction');
  3 Comments
ah sho
ah sho on 5 Jan 2020
f0 = val;
fs = 360;
Tw = 10;
figure(1)
plot(f0);
title('noisy signal');
B0 = strel('line',0.2*fs*Tw,0);
Bc = strel('line',1.5*0.2*fs*Tw,0);
fb = imopen(f0, B0);
fb = imclose(fb, Bc);
fbc = f0 - fb;
figure(2)
plot(fbc);
title('baseline correction');

Sign in to comment.

Answers (2)

Max Murphy
Max Murphy on 5 Jan 2020
Edited: Max Murphy on 5 Jan 2020
Here are some ways of looking at it (although by only specifying baseline I am not sure I know what part of the signal you would like to remove, as that is application-specific):
clear; clc; close all force;
%% Load data
in = load('100m.mat','val');
f0 = in.val; % Units? I am assuming mV
% Compute values
fs = 360; % Sample rate (Hz)
t = (0:(numel(f0)-1))/fs; % Time (seconds)
% Tw - max(t) (seconds); adjusted B0 and Bc values to reflect
%% Processing attempt 1 (original baseline-removal procedure)
% Trying to understand your algorithm here (my comments)
% --> Just note from my perspective "baseline removal" is a confusing term.
% Baseline could refer to many things, such as some period prior to
% your experiment that you are using as a control for the signal you
% observe in a 10-second window, for example. So I am interpreting as
% "offset removal" instead, since it seems you want to remove the drift
% between ECG spikes in your signal. In the fourth figure I do the
% opposite and remove the spikes instead of the slow part.
% Create 'line' structuring elements of fixed "timescales"
B0 = strel('line',2*fs,0); % 2 seconds ("fast timescale")
Bc = strel('line',3*fs,0); % 3 seconds ("slow timescale")
fb_fast = imopen(f0, B0); % Morphologically opens "image"
fb_slow = imclose(fb_fast, Bc); % Morphologically closes "opened image"
fbc_orig = f0 - fb_slow; % Subtracts computed reference value
% Plot
fig_orig = figure('Name','ECG Baseline Removal Overlay: Original',...
'Color','w',...
'Units','Normalized',...
'Position',[0.1 0.5 0.3 0.3]);
ax = axes(fig_orig,'NextPlot','add');
plot(ax,t,f0,'r','LineWidth',2,'DisplayName','Noisy');
plot(ax,t,fb_fast,'k--','LineWidth',1,'DisplayName','Ref_{fast}');
plot(ax,t,fb_slow,'k:','LineWidth',1,'DisplayName','Ref_{slow}');
plot(ax,t,fbc_orig,'b','LineWidth',2,'DisplayName','Corrected');
title(ax,'Offset Removal (Posted Algorithm)',...
'FontName','Arial','Color','k');
xlabel(ax,'Time (sec)','FontName','Arial','Color','k');
ylabel(ax,'Amplitude (mV)','FontName','Arial','Color','k');
legend(ax,{'Noisy','Ref_{fast}','Ref_{slow}','Corrected'});
%% Processing attempt 2 (remove median - simplest)
% Yours is still better than this one, but just to illustrate:
fb_med = median(f0); % Compute flat median of entire signal
fbc_med = f0 - fb_med; % Simple median subtraction
% Plot
fig_med = figure('Name','ECG Baseline Removal Overlay: Original',...
'Color','w',...
'Units','Normalized',...
'Position',[0.1 0.1 0.3 0.3]);
ax = axes(fig_med,'NextPlot','add');
plot(ax,t,f0,'r','LineWidth',2,'DisplayName','Noisy');
plot(ax,t,fb_med,'k:','LineWidth',1,'DisplayName','Ref_{median}');
plot(ax,t,fbc_med,'b','LineWidth',2,'DisplayName','Corrected');
title(ax,'Offset Removal (Median Algorithm)',...
'FontName','Arial','Color','k');
xlabel(ax,'Time (sec)','FontName','Arial','Color','k');
ylabel(ax,'Amplitude (mV)','FontName','Arial','Color','k');
legend(ax,{'Noisy','Ref_{median}','Corrected'});
%% Processing attempt 3 (remove median in sliding window)
% Maybe this is better, you can tweak nSeconds for width of window:
nSeconds = 1; % 1-second "timescale" for median
fb_medw = medfilt1(f0,nSeconds*fs); % get median in sliding window
fbc_medw = f0 - fb_medw; % subtract median from signal
% Reducing nSeconds will remove median on a faster time-scale, effectively
% increasing the high-pass filter effect.
% Plot
fig_medw = figure('Name','ECG Baseline Removal Overlay: Original',...
'Color','w',...
'Units','Normalized',...
'Position',[0.5 0.5 0.3 0.3]);
ax = axes(fig_medw,'NextPlot','add');
plot(ax,t,f0,'r','LineWidth',2,'DisplayName','Noisy');
plot(ax,t,fb_medw,'k:','LineWidth',1,'DisplayName','Ref_{median}');
plot(ax,t,fbc_medw,'b','LineWidth',2,'DisplayName','Corrected');
title(ax,'Offset Removal (Window-Median Algorithm)',...
'FontName','Arial','Color','k');
xlabel(ax,'Time (sec)','FontName','Arial','Color','k');
ylabel(ax,'Amplitude (mV)','FontName','Arial','Color','k');
legend(ax,{'Noisy','Ref_{window-median}','Corrected'});
%% Processing attempt 4 (remove "spiky" part of the signal?)
% Maybe you mean the Baseline is the spiky part and you would like to
% remove that? The "spiky" part is higher in frequency-content than the
% majority of your signal. You are limited by Nyquist frequency to only
% observe signals of 160 Hz or less. The fluctuations in "non-spiky" part
% are in much lower frequency domain than that. Let's remove the
% high-frequency content from our signal and see how it looks:
% Filter design
fc = 5; % Highpass cutoff frequency (Hz) - application specific
Wn = fc / (fs/2); % Normalized frequency
% If you need "more of the spike" removed, then reduce fc
[b,a] = butter(4,Wn,'high'); % 'high'-->Keep "fast" component of signal
% filter -> causal implementation; filtfilt -> no phase offset
fb_despike = filtfilt(b,a,f0);
% Note that fb_despike could be used as an endpoint for the previous
% attempts. Or you could achieve this result using a different filter
% specification
% >> [b,a] = butter(4,Wn,'low'); % 'low'-->Keep "slow" component of signal
fbc_despike = f0 - fb_despike;
% Plot
fig_despike = figure('Name','ECG Baseline Removal Overlay: Original',...
'Color','w',...
'Units','Normalized',...
'Position',[0.5 0.1 0.3 0.3]);
ax_top = axes(fig_despike,'NextPlot','add',...
'Units','Normalized','Position',[0.05 0.60 0.9 0.35]);
ax_bot = axes(fig_despike,'NextPlot','add',...
'Units','Normalized','Position',[0.05 0.10 0.9 0.35]);
plot(ax_top,t,f0,'r','LineWidth',2,'DisplayName','Noisy');
plot(ax_top,t,fb_despike,'k:','LineWidth',1,'DisplayName','Ref_{de-spike}');
plot(ax_top,t,fbc_despike,'b','LineWidth',2,'DisplayName','Corrected');
title(ax_top,'Offset Removal (De-Spike Algorithm)',...
'FontName','Arial','Color','k');
xlabel(ax_top,'Time (sec)','FontName','Arial','Color','k');
ylabel(ax_top,'Amplitude (mV)','FontName','Arial','Color','k');
legend(ax_top,{'Noisy','Ref_{de-spike}','Corrected'});
copyobj(get(ax_top,'Children'),ax_bot);
xlim(ax_bot,[0.5 2.5]); % zoom in
title(ax_bot,'Zoomed View',...
'FontName','Arial','Color','k');
  1 Comment
Max Murphy
Max Murphy on 10 Jan 2020
% Assuming T is all the samples in your signal
tStart = 1;
T = numel(fbc_orig);
% (note: it could be an arbitrary length that is shorter)
% tStart = sample;
% T = 20; % --> Window that is 20 samples long
% If fbc_orig is a vector:
BCR = sum(abs(fbc_orig(tStart:(tStart+T-1)))) / ...
sum(abs(f0(tStart:(tStart+T-1))));
% If fbc_orig is a matrix (e.g. sample n has multiple data points):
% (note: this assumes that data is organized as an N x M matrix,
% with N data samples (time points) and M variables)
BCR = sum(sqrt(sum(fbc_orig(tStart:(tStart+T-1),:).^2,2))) / ...
sum(sqrt(sum(f0(tStart:(tStart+T-1),:).^2,2)));

Sign in to comment.


Image Analyst
Image Analyst on 22 Feb 2020
Try this. Adjust findpeaks() parameters as needed.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clearvars;
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 18;
s = load('100m.mat')
y = s.val;
x = 1 : length(y);
subplot(2, 1, 1);
plot(x, y, 'b-', 'LineWidth', 2);
xlabel('Index', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
grid on;
[valleyValues, valleyIndexes] = findpeaks(-y, 'MinPeakheight', 60, 'MinPeakDistance', 30);
valleyValues = -valleyValues;
hold on
plot(valleyIndexes, valleyValues, 'rv');
title('Original Signal, with Negative-going peaks shown', 'FontSize', fontSize);
% Get interpolated baseline
baseLine = interp1(valleyIndexes, valleyValues, x);
plot(x, baseLine, 'm-', 'LineWidth', 2);
% Subtract the baseline to get the baseline corrected signal.
correctedSignal = y - baseLine;
subplot(2, 1, 2);
plot(x, correctedSignal, 'b-', 'LineWidth', 2);
grid on;
xlabel('Index', 'FontSize', fontSize);
ylabel('y Corrected', 'FontSize', fontSize);
title('Baseline-Corrected Signal', 'FontSize', fontSize);
  1 Comment
Max Murphy
Max Murphy on 23 Feb 2020
Edited: Max Murphy on 23 Feb 2020
Nice!
You may want to increase the 'MinPeakDIstance' parameter a little bit, as it looks like the baseline correction may introduce a spurious effect at around index 1300 (eye-balling it). I do like the approach though.

Sign in to comment.

Categories

Find more on Get Started with Signal Processing Toolbox in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!