function [ y ] = ATVMF( x, interp_method, operator )
% Algorithm name: Adaptive Time-Varying Morphological Filtering (ATVMF)
%
% Algorithm description: This method can achieve adaptive morphological filtering,
% in which a time-varying structure element (TVSE) is adaptively designed
% based on the characteristics of a signal and is no longer fixed.
%
% Input:
% x: signal to be analyzed (a vector)
% interp_method: selected interpllation method, such as 'spline', 'pchip',
% 'linear' and 'nearest', in which 'spline' is recommended.
% 'spline' -- cubic spline interpolation
% 'pchip' -- cubic Hermitian interpolation
% 'linear' -- piecewise linear interpolation
% 'nearest' -- nearest neighbor interpolation
% operator: selected morphological operator, see sub-function 'MF_operator'
%
% Output:
% y: morphological filtered signal
x = x(:)-mean(x); % a vector
N = length(x);
[indmin, indmax] = extreme_points(x); % Determine the location of local minima and maxima
% indmin -- the position of the local minimum point in the sequence x
% indmax -- the position of the local maximum point in the sequence x
tmin = indmin;
tmax = indmax;
xmin = x(tmin); % The magnitude of the local minimum point
xmax = x(tmax); % The magnitude of the local maximum point
% Sorting of local minimum and maximum points
textra = zeros(1,length(tmin)+length(tmax));
xextra = zeros(1,length(xmin)+length(xmax));
if tmin(1) < tmax(1) % The first extreme point is the minimum point
if tmin(end) > tmax(end) % The last extreme point is the minimum point
textra(1) = tmin(1);
xextra(1) = xmin(1);
for i = 1:length(tmax)
textra(2*i) = tmax(i);
textra(2*i+1) = tmin(i+1);
xextra(2*i) = xmax(i);
xextra(2*i+1) = xmin(i+1);
end
else % The last extreme point is the maximum point
for i = 1:length(tmax)
textra(2*i-1) = tmin(i);
textra(2*i) = tmax(i);
xextra(2*i-1) = xmin(i);
xextra(2*i) = xmax(i);
end
end
else % The first extreme point is the maximum point
if tmin(end) < tmax(end) % The last extreme point is the maximum point
textra(1) = tmax(1);
xextra(1) = xmax(1);
for i = 1:length(tmin)
textra(2*i) = tmin(i);
textra(2*i+1) = tmax(i+1);
xextra(2*i) = xmin(i);
xextra(2*i+1) = xmax(i+1);
end
else % The last extreme point is the minimum point
for i = 1:length(tmin)
textra(2*i-1) = tmax(i);
textra(2*i) = tmin(i);
xextra(2*i-1) = xmax(i);
xextra(2*i) = xmin(i);
end
end
end
% Selection of 'interp_method'
env = interp1(textra,xextra,textra(1):textra(end),interp_method);
delta = textra(1)-1;
S = length(indmin)-1; % number of SE
y = []; % output initialization
for s = 1:S
xnew = x(indmin(s)+1:indmin(s+1));
g = env(indmin(s)+1-delta:indmin(s+1)-delta);
g = g-min(g);
% the morphological filtering result
ynew = MF_operator( xnew, g, operator );
y = [y; ynew];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sub-function
function d = dilation(f,g)
% Morphological dilation operation
N = length(f);
M = length(g);
dtmp = f;
for i = 1:N
for j = 1:M
if (i-j) >= 1 && (i-j) <= N
tmp = f(i-j) + g(j);
if tmp > dtmp(i)
dtmp(i) = tmp;
end
end
end
end
d = dtmp;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sub-function
function e = erosion(f,g)
% Morphological erosion operation
N = length(f);
M = length(g);
dtmp = f;
for i = 1:N
for j = 1:M
if (i+j) >= 1 && (i+j) <= N
tmp = f(i+j) - g(j);
if tmp < dtmp(i)
dtmp(i) = tmp;
end
end
end
end
e = dtmp;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sub-function
function y = MF_operator( x, g, operator )
% Morphological operators
%
a1 = dilation(x,g); % dilation
a2 = erosion(a1,g); % closing
a3 = erosion(a2,g);
a4 = dilation(a3,g); % closing-opening
%
b1 = erosion(x,g); % erosion
b2 = dilation(b1,g); % opening
b3 = dilation(b2,g);
b4 = erosion(b3,g); % opening-closing
if strcmp(operator,'Gde') == 1
y = a1-b1;
elseif strcmp(operator,'Gco') == 1
y = a2-b2;
elseif strcmp(operator,'Gcooc') == 1
y = a4-b4;
elseif strcmp(operator,'AHde') == 1
y = x-(a1+b1)/2;
elseif strcmp(operator,'AHco') == 1
y = x-(a2+b2)/2;
elseif strcmp(operator,'AHcooc') == 1
y = x-(a4+b4)/2;
elseif strcmp(operator,'MGPO1') == 1
y = (a1-b1).*(a2-b2);
elseif strcmp(operator,'MGPO2') == 1
y = (a1-b1).*(a4-b4);
elseif strcmp(operator,'MGPO3') == 1