# morphological operation and filtering

3 vues (au cours des 30 derniers jours)
ah sho le 5 Jan 2020
Modifié(e) : Max Murphy le 23 Fév 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 commentairesAfficher 1 commentaire plus ancienMasquer 1 commentaire plus ancien
ah sho le 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');
Image Analyst le 22 Fév 2020

Connectez-vous pour commenter.

### Réponses (2)

Max Murphy le 5 Jan 2020
Modifié(e) : Max Murphy le 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;
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)
% --> 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]);
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]);
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]);
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]);
'Units','Normalized','Position',[0.05 0.60 0.9 0.35]);
'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 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Max Murphy le 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)));

Connectez-vous pour commenter.

Image Analyst le 22 Fév 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;
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 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Max Murphy le 23 Fév 2020
Modifié(e) : Max Murphy le 23 Fév 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.

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Get Started with Signal Processing Toolbox dans Help Center et File Exchange

R2015b

### Community Treasure Hunt

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

Start Hunting!

Translated by