Homec4science
Diffusion Simon Botrylloides Regeneration Group (master) Browse

SiBoRG/libraries/local_extrema.m
dbb5add19cbamaster

Last Change
Raw File

local_extrema.m
View Options

This document is not UTF8. It was detected as Shift JIS and converted to UTF8 for display.
function [xmax,imax,xmin,imin] = local_extrema(x, nneigh)
% LOCAL_EXTREMA Gets the local extrema points from a time series.
%
% [XMAX,IMAX,XMIN,IMIN] = LOCAL_EXTREMA(X) returns the local minima
% and maxima points of the vector X ignoring NaN's, where
% XMAX - maxima points
% IMAX - indexes of the XMAX
% XMIN - minima points
% IMIN - indexes of the XMIN
%
% [...] = LOCAL_EXTREMA(X, NNEIGHBORS) returns only the maximal/minimal
% extrama values in a [-NNEIGHBORS +NNEIGHBORS] local window.
%
% Based on EXTRAM by
% Lic. on Physics Carlos Adri疣 Vargas Aguilera
% Physical Oceanography MS candidate
% UNIVERSIDAD DE GUADALAJARA
% Mexico, 2004
%
% Wilson lab, University of Otago
% Simon Blanchoud
% 25.06.2015
% No local average
if (nargin < 2)
nneigh = 1;
end
xmax = [];
imax = [];
xmin = [];
imin = [];
% Vector input?
Nt = numel(x);
if Nt ~= length(x)
error('Entry must be a vector.')
end
% NaN's:
inan = find(isnan(x));
indx = 1:Nt;
if ~isempty(inan)
indx(inan) = [];
x(inan) = [];
Nt = length(x);
end
% Difference between subsequent elements:
x = x(:);
dx = differentiator(x, 'super');
% Local maxima and minima for windowing
nneigh = ceil(nneigh);
if (nneigh > 1)
window = 2*nneigh + 1;
indxs = [1:nneigh nneigh+2:window].';
ssize = [window, 1];
lmax = colfilt(x, ssize, 'sliding', @(y)(max(y(indxs,:), [], 1)));
if (size(x, 1) == 1)
lmax = lmax.';
end
% Only if needed
if (nargout > 2)
lmin = colfilt(x, ssize, 'sliding', @(y)(min(y(indxs,:), [], 1)));
if (size(x, 1) == 1)
lmin = lmin.';
end
end
else
lmax = x-1;
lmin = x+1;
end
% Is an horizontal line?
if ~any(dx)
return
end
% Flat peaks? Put the middle element:
a = find(dx~=0); % Indexes where x changes
lm = find(diff(a)~=1) + 1; % Indexes where a do not changes
d = a(lm) - a(lm-1); % Number of elements in the flat peak
a(lm) = a(lm) - floor(d/2); % Save middle elements
a(end+1) = Nt;
% Peaks?
xa = x(a); % Serie without flat peaks
b = (diff(xa) > 0); % 1 => positive slopes (minima begin)
% 0 => negative slopes (maxima begin)
xb = diff(b); % -1 => maxima indexes (but one)
% +1 => minima indexes (but one)
a = a(1:end-1);
% Compensate for the missing datapoints after diff
if (size(x,1)==1)
xb = [0 xb];
else
xb = [0; xb];
end
imax = find((xb == -1) & (x(a) > lmax(a))); % maxima indexes
imax = a(imax);
% Only if needed
if (nargout > 2)
imin = find((xb == +1) & (x(a) < lmin(a))); % minima indexes
imin = a(imin);
else
imin = [];
end
nmaxi = length(imax);
nmini = length(imin);
% Maximum or minumim on a flat peak at the ends?
if (nmaxi==0) && (nmini==0)
if x(1) > x(Nt)
xmax = x(1);
imax = indx(1);
xmin = x(Nt);
imin = indx(Nt);
elseif x(1) < x(Nt)
xmax = x(Nt);
imax = indx(Nt);
xmin = x(1);
imin = indx(1);
end
return
end
xmax = x(imax);
xmin = x(imin);
% NaN's:
if (nargout > 1)
if ~isempty(inan)
imax = indx(imax);
imin = indx(imin);
end
% Same size as x:
imax = reshape(imax,size(xmax));
imin = reshape(imin,size(xmin));
end
return;
end
c4science · Help