morphological operation and filtering

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

[Adam Danz]

We can't run your code due to missing variables (ie, fs, Tw, etc). I realize that some of those variable may be extracted from the txt file but it's better if you set up that code so we can easily copy-paste it and we can be certain that we're working with the same values.

[ah sho]

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

[Image Analyst]

Answer 1

[Max Murphy]

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

[Max Murphy]

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

Answer 2

[Image Analyst]

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

[Max Murphy]

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.