morphological operation and filtering

조회 수: 3 (최근 30일)
ah sho 2020년 1월 5일
편집: Max Murphy 2020년 2월 23일
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이전 댓글 1개 표시이전 댓글 1개 숨기기
ah sho 2020년 1월 5일
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 2020년 2월 22일

댓글을 달려면 로그인하십시오.

답변 (2개)

Max Murphy 2020년 1월 5일
편집: Max Murphy 2020년 1월 5일
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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Max Murphy 2020년 1월 10일
% 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)));

댓글을 달려면 로그인하십시오.

Image Analyst 2020년 2월 22일
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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Max Murphy 2020년 2월 23일
편집: Max Murphy 2020년 2월 23일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Get Started with Signal Processing Toolbox에 대해 자세히 알아보기

R2015b

Community Treasure Hunt

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

Start Hunting!

Translated by