Dilatation and Erosion as custom function for signal processing
显示 更早的评论
I want to define the two function for accomplish the morphological operation of erosion and dilatation on a 1-D signal (in this case an ecg signal). I know that I can use functions like imerode and imdilate, but I want to define them by myself, obtaining a result comparable to the functions mentioned before.
To design the two function I have to followed the example propose below:
Based on this knowledge I have defined two functions: dilatation and erosion (below). However when I apply the algorithm (show below in the code) the result is something that I do not understand. Based on the methods presented above I can't understand what I is missing or what is wrong in the two functions.
My Code:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
% Signal info.
G = 200; % Gain
Fs = 360; % [Hz]
L = 3600; % lenght of ECG signals
T = linspace(0,L/Fs,L); % time axis
F = linspace(-Fs/2, Fs/2, L); % Frequency axis
% Load Signal
load("100m.mat");
ecg = val/G;
% Plot the ecg signal
figure(1); plot(T, ecg); grid on;
title("ECG Signal","FontSize",fontSize);
xlabel("Time (sec)", "FontSize", fontSize);
ylabel("voltage [mV]", "FontSize", fontSize);
% Apply MMF algorithm
% Declaration of the two Structuring Elements
B1 = [0 1 5 1 0]; B2 = [1 1 1 1 1];
% modified morphological filtering algorithm
d1 = dilatation(ecg, B1); % d1 = imdilate(ecg, B1);
e1 = erosion(d1,B2); % e1 = imerode(ecg, B2);
e2 = erosion(ecg,B1); % e2 = imerode(ecg, B1);
d2 = dilatation(e2,B2); % d2 = imdilate(e2, B2);
mmf = (e1+d2)/2; % average
% Plot
figure(2); plot(T,mmf,"b-"); grid on;
title("MMF Denoised Ecg signal", "FontSize", fontSize);
xlabel("Time (sec)", "FontSize", fontSize);
ylabel("Voltage (Hz)", "FontSize", fontSize);
Dilatation Function:
function [output] = dilatation(inputSignal,structuringElement)
% Dilatation operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
output = zeros(length(inputSignal));
m = length(structuringElement);
n = length(inputSignal);
for i = abs((m-1)/2):n-abs((m+1)/2)-1
for j=1:m
output(i) = max(output(i-1), inputSignal(i-abs((m-1)/2)+j)+structuringElement(j));
end
end
end
Erosion Function:
function [output] = erosion(inputSignal,structuringElement)
% Erosion operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
output = zeros(length(inputSignal));
m = length(structuringElement);
n = length(inputSignal);
for i = abs((m-1)/2):n-abs((m+1)/2)-1
for j=1:m
output(i) = min(output(i-1), inputSignal(i-abs((m-1)/2)+j)-structuringElement(j));
end
end
end
采纳的回答
The loop part can be cleaned up, but a big part of why the output looks messed up starts here:
output = zeros(length(inputSignal));
This will create a square matrix, not a vector. So then this assignment inside the loop:
output(i) = min(...);
... will populate the first column. Then if you feed the result to plot(), it plots each row as a series.
So just make sure to preallocate a vector.
I'm going to assume that this is closer to what you're looking for.
insignal = [0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0];
se = [1 1 1];
outsignal = erosion(insignal,se);
subplot(2,1,1)
plot(insignal)
subplot(2,1,2)
plot(outsignal)
function [output] = erosion(inputSignal,structuringElement)
% Erosion operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
% get sizes
M = length(structuringElement);
N = length(inputSignal);
% this relies on the strel being odd width
if mod(M,2)==0
error('strel length must be odd')
end
% perform dilation
hw = floor(M/2); % half-width of strel
output = zeros(1,N);
for n = ceil(M/2):N-ceil(M/2)+1
output(n) = min(inputSignal(n-hw:n+hw) - structuringElement);
end
end
Bear in mind that the behavior at the ends is due to the way the operation is constrained there. The alternative to avoiding the ends of the signal vector is to add enough padding to them that you can filter the whole vector.
9 个评论
Sabaudian
2022-7-11
DGM
2022-7-11
The output array is initialized with zeros, but the calculated values do not extend all the way to the end of the output vector due to the width of the strel. As a result, the output vector has (because the strel was width 3) one element at each end that still has the original zero in it.
Assume that the strel was width 5. There would be two elements at each end which don't get processed:
The alternative is to pad the array with floor(M/2) extra values at each end. That way you can process all the values in the signal vector. The padding can be done either by simple concatenation, or by using padarray(). By adjusting the loop indices accordingly, the output array size remains the same.
There are other alternatives that can be imagined, and what's appropriate depends on the intent. Just pay attention to how the strel overlaps the signal vector.
Something like this.
insignal = [0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0];
se = [1 1 1];
outsignal = erosion_padded(insignal,se);
subplot(2,1,1)
plot(insignal)
subplot(2,1,2)
plot(outsignal)
function [output] = erosion_padded(inputSignal,structuringElement)
% Erosion operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
% get sizes
M = length(structuringElement);
N = length(inputSignal);
% this relies on the strel being odd width
if mod(M,2)==0
error('strel length must be odd')
end
% pad vector
hw = floor(M/2); % half-width of strel
inputSignal = padarray(inputSignal,[0 hw],0,'both');
% perform dilation
output = zeros(1,N); % same length as input
for n = hw+1:hw+N % n = (1:N)+hw (indexes over the same length as the input, but offset by +hw)
% the indexing on the LHS removes that +hw offset
output(n-hw) = min(inputSignal(n-hw:n+hw) - structuringElement);
end
end
Also bear in mind that everything so far assumes that the input signal and strel are row vectors.
Sabaudian
2022-7-13
Ok, I have tried your solution and it seems to give a better result but not optimal lets say, there is still some part of the signal that is not properly included into the operations (erosion and dilatation), causing a "distortion" in the output at the start and at the end, as it is possible to see in the image below
So to resolve this, what I am supposed to to? Enlarge the Padding? or something else?
DGM
2022-7-13
If you're still padding prior to calling the erosion/dilation functions, don't. The padding is handled within the functions themselves.
Yes I have already delete that part when I applied your solution. the only padding is in the functions.
The code is the same as before (minus the padding part) and as you can see the problem persist.
% Signal info.
G = 200; % Gain
Fs = 360; % [Hz]
L = 3600; % lenght of ECG signals
T = linspace(0,L/Fs,L); % time axis
F = linspace(-Fs/2, Fs/2, L); % Frequency axis
% Load Signal
load("100m.mat");
ecg = val/G;
% Distort the Baseline aadding drift
X = linspace(0,20*pi,L);
A = 2; % ampl. of upward or downward drift
N = 30; % severity of baseline roll
sn = A*cos(X./N);
minDriftOffset = 1;
maxDriftOffset = 5;
slantedLine = linspace(minDriftOffset,maxDriftOffset,L);
drift = sn + slantedLine;
% Apply drift to all signals
dft_ecg = ecg + drift;
%% Baseline correction
close all;
% Defining the two structuring element
Bo = ones(1,73);
Bc = ones(1,145);
% Detection and Correction of the Wandering Baseline
peaksSuppression = dilatation(erosion(dft_ecg,Bo),Bo);
pitsSuppression = erosion(dilatation(peaksSuppression,Bc),Bc);
detectedDrift = pitsSuppression;
% Correction subtracting the drift from signals
Correction = dft_ecg - detectedDrift;
% Final Baseline result
finalBaseline = erosion(dilatation(dilatation(erosion(Correction,Bo),Bo),Bc),Bc);
% Plot Corrected Signal and Baseline
figure(1); subplot(2,1,1); hold on;
plot(T,detectedDrift,"g-","LineWidth",3);
plot(T,dft_ecg, "b-","LineWidth",0.5);
title("Detected Baseline Drift");
% legend("baseline","signal");
xlabel("Time (sec)");
ylabel("Voltage (Hz)");
grid on; hold off;
subplot(2,1,2); hold on;
plot(T,finalBaseline,"g-","LineWidth",3);
plot(T,Correction, "b-","LineWidth",0.5);
title("Drift Correction");
% legend("baseline","signal");
xlabel("Time (sec)");
ylabel("Voltage (Hz)");
grid on; hold off;
Augh. It's picking up the padding anyway. Try doing this instead of padding with zeros:
input = padarray(input,[0 hw],'replicate','both');
I am kind of confused by the behavior of this method of dilation/erosion. I'm used to doing this with images, where the strel is multiplied instead of added/subtracted. It's just not a convention I'm used to thinking about. At least with the above change, it seems to work okay. See what you think.
Sabaudian
2022-7-14
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Bartlett 的更多信息
Translated by 
