How to detect point before a positive ramp

26 Ansichten (letzte 30 Tage)
Francesco Lucarelli
Francesco Lucarelli am 15 Jul. 2021
Kommentiert: Rik am 16 Jul. 2021
Hi all,
I have a stimulation signal (current mA) that is like the picture below:
So it is made by peaks. Please find below a zoom of one pulse:
My goal is to find the point in the red circle, so the last point before the positve ramp of the peak
How can i do that?
Thanks a lot for your help and time
KR

Akzeptierte Antwort

Mathieu NOE
Mathieu NOE am 15 Jul. 2021
hello
I have combined threshold crossing detection with second order derivative
this is the result on a noisy pulse train ....
the entire code , including my prefered functions , is below :
clc
clearvars
%% dummy data
% repetition frequency of 3 Hz and a sawtooth width of 0.1 sec.
% The signal is to be 1 second long with a sample rate of 1kHz.
Fs = 1e3;
t = 0 : 1/Fs : 1; % 1 kHz sample freq for 1 sec
D = 0 : 1/5 : 1; % 5 Hz repetition freq
y1 = pulstran(t,D,'rectpuls',0.02);
y1 = y1 + 0.03*randn(size(y1)); % add some noise
%% step 1 : threshold crossing points - get an approximative time stamps where to look for
threshold = 0.5; % your value here
[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y1,t,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
%% step 2 : second derivative peaks to find in the neiborhood of first selection of threshold crossing points
[dy, ddy] = firstsecondderivatives(t,y1);
sel = (max(ddy)-min(ddy))/4;
% sel - The amount above surrounding data for a peak to be
% identified (default = (max(x0)-min(x0))/4). Larger values mean
% the algorithm is more selective in finding peaks.
thresh = threshold*Fs^2;
extrema = 1; %if maxima are desired
[peak_ind] = peakfinder(ddy,sel,thresh,extrema) ; % or findpeaks if you like it...
time_peak_ind = t(peak_ind);
%% step 3 : select only peaks that are closest to crossing points (in first section)
for ci = 1:numel(t0_pos1);
dist = abs(time_peak_ind-t0_pos1(ci));
[M,I] = min(dist);
t0_pos1_selected(ci) = time_peak_ind(I);
end
y1_selected = interp1(t,y1,t0_pos1_selected);
figure(1)
subplot(211),plot(t,y1,'b',t0_pos1,s0_pos1,'dg',t0_pos1_selected,y1_selected,'dr','linewidth',2,'markersize',12);
legend('signal 1','positive slope crossing points',' "kirk" points' );
xlabel('Time(s)');
subplot(212),plot(t,ddy,'r','linewidth',2,'markersize',12);
legend('second derivative' );
xlabel('Time(s)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
%
% [ind,t0,s0] = ... also returns the data vector corresponding to
% the t0 values.
%
% [ind,t0,s0,t0close,s0close] additionally returns the data points
% closest to a zero crossing in the arrays t0close and s0close.
%
% This version has been revised incorporating the good and valuable
% bugfixes given by users on Matlabcentral. Special thanks to
% Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
% Zach Lewis for their input.
% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27 revised version
% Copyright (c) Steffen Brueckner, 2002-2007
% brueckner@sbrs.net
% M Noe
% added positive or negative slope condition
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) > eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) > eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
end
function [dy, ddy] = firstsecondderivatives(x,y)
% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.
n = length (x);
dy = zeros;
ddy = zeros;
% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1) = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1))); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative
for i = 2:n-1
dy(i) = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end
%%%%%%%%%%%%%%%%
function varargout = peakfinder(x0, sel, thresh, extrema, include_endpoints)
%PEAKFINDER Noise tolerant fast peak finding algorithm
% INPUTS:
% x0 - A real vector from the maxima will be found (required)
% sel - The amount above surrounding data for a peak to be
% identified (default = (max(x0)-min(x0))/4). Larger values mean
% the algorithm is more selective in finding peaks.
% thresh - A threshold value which peaks must be larger than to be
% maxima or smaller than to be minima.
% extrema - 1 if maxima are desired, -1 if minima are desired
% (default = maxima, 1)
% include_endpoints - If true the endpoints will be included as
% possible extrema otherwise they will not be included
% (default = true)
% OUTPUTS:
% peakLoc - The indicies of the identified peaks in x0
% peakMag - The magnitude of the identified peaks
%
% [peakLoc] = peakfinder(x0) returns the indicies of local maxima that
% are at least 1/4 the range of the data above surrounding data.
%
% [peakLoc] = peakfinder(x0,sel) returns the indicies of local maxima
% that are at least sel above surrounding data.
%
% [peakLoc] = peakfinder(x0,sel,thresh) returns the indicies of local
% maxima that are at least sel above surrounding data and larger
% (smaller) than thresh if you are finding maxima (minima).
%
% [peakLoc] = peakfinder(x0,sel,thresh,extrema) returns the maxima of the
% data if extrema > 0 and the minima of the data if extrema < 0
%
% [peakLoc, peakMag] = peakfinder(x0,...) returns the indicies of the
% local maxima as well as the magnitudes of those maxima
%
% If called with no output the identified maxima will be plotted along
% with the input data.
%
% Note: If repeated values are found the first is identified as the peak
%
% Ex:
% t = 0:.0001:10;
% x = 12*sin(10*2*pi*t)-3*sin(.1*2*pi*t)+randn(1,numel(t));
% x(1250:1255) = max(x);
% peakfinder(x)
%
% Copyright Nathanael C. Yoder 2011 (nyoder@gmail.com)
% Perform error checking and set defaults if not passed in
error(nargchk(1,5,nargin,'struct'));
error(nargoutchk(0,2,nargout,'struct'));
s = size(x0);
flipData = s(1) < s(2);
len0 = numel(x0);
if len0 ~= s(1) && len0 ~= s(2)
error('PEAKFINDER:Input','The input data must be a vector')
elseif isempty(x0)
varargout = {[],[]};
return;
end
if ~isreal(x0)
warning('PEAKFINDER:NotReal','Absolute value of data will be used')
x0 = abs(x0);
end
if nargin < 2 || isempty(sel)
sel = (max(x0)-min(x0))/4;
elseif ~isnumeric(sel) || ~isreal(sel)
sel = (max(x0)-min(x0))/4;
warning('PEAKFINDER:InvalidSel',...
'The selectivity must be a real scalar. A selectivity of %.4g will be used',sel)
elseif numel(sel) > 1
warning('PEAKFINDER:InvalidSel',...
'The selectivity must be a scalar. The first selectivity value in the vector will be used.')
sel = sel(1);
end
if nargin < 3 || isempty(thresh)
thresh = [];
elseif ~isnumeric(thresh) || ~isreal(thresh)
thresh = [];
warning('PEAKFINDER:InvalidThreshold',...
'The threshold must be a real scalar. No threshold will be used.')
elseif numel(thresh) > 1
thresh = thresh(1);
warning('PEAKFINDER:InvalidThreshold',...
'The threshold must be a scalar. The first threshold value in the vector will be used.')
end
if nargin < 4 || isempty(extrema)
extrema = 1;
else
extrema = sign(extrema(1)); % Should only be 1 or -1 but make sure
if extrema == 0
error('PEAKFINDER:ZeroMaxima','Either 1 (for maxima) or -1 (for minima) must be input for extrema');
end
end
if nargin < 5 || isempty(include_endpoints)
include_endpoints = true;
else
include_endpoints = boolean(include_endpoints);
end
x0 = extrema*x0(:); % Make it so we are finding maxima regardless
thresh = thresh*extrema; % Adjust threshold according to extrema.
dx0 = diff(x0); % Find derivative
dx0(dx0 == 0) = -eps; % This is so we find the first of repeated values
ind = find(dx0(1:end-1).*dx0(2:end) < 0)+1; % Find where the derivative changes sign
% Include endpoints in potential peaks and valleys is desired
if include_endpoints
x = [x0(1);x0(ind);x0(end)];
ind = [1;ind;len0];
else
x = x0(ind);
end
% x only has the peaks, valleys, and possibly endpoints
len = numel(x);
minMag = min(x);
if len > 2 % Function with peaks and valleys
% Set initial parameters for loop
tempMag = minMag;
foundPeak = false;
leftMin = minMag;
if include_endpoints
% Deal with first point a little differently since tacked it on
% Calculate the sign of the derivative since we taked the first point
% on it does not neccessarily alternate like the rest.
signDx = sign(diff(x(1:3)));
if signDx(1) <= 0 % The first point is larger or equal to the second
if signDx(1) == signDx(2) % Want alternating signs
x(2) = [];
ind(2) = [];
len = len-1;
end
else % First point is smaller than the second
if signDx(1) == signDx(2) % Want alternating signs
x(1) = [];
ind(1) = [];
len = len-1;
end
end
end
% Skip the first point if it is smaller so we always start on a
% maxima
if x(1) > x(2)
ii = 0;
else
ii = 1;
end
% Preallocate max number of maxima
maxPeaks = ceil(len/2);
peakLoc = zeros(maxPeaks,1);
peakMag = zeros(maxPeaks,1);
cInd = 1;
% Loop through extrema which should be peaks and then valleys
while ii < len
ii = ii+1; % This is a peak
% Reset peak finding if we had a peak and the next peak is bigger
% than the last or the left min was small enough to reset.
if foundPeak
tempMag = minMag;
foundPeak = false;
end
% Make sure we don't iterate past the length of our vector
if ii == len
break; % We assign the last point differently out of the loop
end
% Found new peak that was lager than temp mag and selectivity larger
% than the minimum to its left.
if x(ii) > tempMag && x(ii) > leftMin + sel
tempLoc = ii;
tempMag = x(ii);
end
% ii = ii+1; % Move onto the valley
% % Come down at least sel from peak
% if ~foundPeak && tempMag > sel + x(ii)
% foundPeak = true; % We have found a peak
% leftMin = x(ii);
% peakLoc(cInd) = tempLoc; % Add peak to index
% peakMag(cInd) = tempMag;
% cInd = cInd+1;
% elseif x(ii) < leftMin % New left minima
% leftMin = x(ii);
% end
jj = ii+1; % Move onto the valley
% Come down at least sel from peak
if ~foundPeak && tempMag > sel + x(jj)
foundPeak = true; % We have found a peak
leftMin = x(jj);
peakLoc(cInd) = tempLoc; % Add peak to index
peakMag(cInd) = tempMag;
cInd = cInd+1;
elseif x(jj) < leftMin % New left minima
leftMin = x(jj);
end
end
% Check end point
if x(end) > tempMag && x(end) > leftMin + sel
peakLoc(cInd) = len;
peakMag(cInd) = x(end);
cInd = cInd + 1;
elseif ~foundPeak && tempMag > minMag % Check if we still need to add the last point
peakLoc(cInd) = tempLoc;
peakMag(cInd) = tempMag;
cInd = cInd + 1;
end
% Create output
peakInds = ind(peakLoc(1:cInd-1));
peakMags = peakMag(1:cInd-1);
else % This is a monotone function where an endpoint is the only peak
[peakMags,xInd] = max(x);
if peakMags > minMag + sel
peakInds = ind(xInd);
else
peakMags = [];
peakInds = [];
end
end
% Apply threshold value. Since always finding maxima it will always be
% larger than the thresh.
if ~isempty(thresh)
m = peakMags>thresh;
peakInds = peakInds(m);
peakMags = peakMags(m);
end
% Rotate data if needed
if flipData
peakMags = peakMags.';
peakInds = peakInds.';
end
% Change sign of data if was finding minima
if extrema < 0
peakMags = -peakMags;
x0 = -x0;
end
% Plot if no output desired
if nargout == 0
if isempty(peakInds)
disp('No significant peaks found')
else
figure;
plot(1:len0,x0,'.-',peakInds,peakMags,'ro','linewidth',2);
end
else
varargout = {peakInds,peakMags};
end
end
  9 Kommentare
Mathieu NOE
Mathieu NOE am 16 Jul. 2021
well
ok , if it get's so complicated, I just stop using my time trying to help people
seems the world is now just a matter of laws , and endless discussion about what belongs to who.
I could also decide that what I post should be protected , which isn't the case.
But I guess we should stop arguing here because this is not a forum about legal fight
Rik
Rik am 16 Jul. 2021
I don't know how this is complicated. You posted something with a comment that claimed copyright over a function. Someone put in enough effort to want to control how others can use that. Now you took that control away. I wouldn't like it if people just took my code and copied it here. I do like it if people get it directly from me.
I don't see why that would mean you should no longer help people. I'm glad you do help people here, and I hope you will continue to do so for a long time.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Matt J
Matt J am 15 Jul. 2021
Bearbeitet: Matt J am 15 Jul. 2021
Using this File Exchange submission
threshold=2.5; %to separate background from signal
locations = groupLims(groupTrue(yourSignal>threshold), 1 )-1
  9 Kommentare
Mathieu NOE
Mathieu NOE am 15 Jul. 2021
what is that matlab function : stim ?
stim(locp(i))>0.5
Francesco Lucarelli
Francesco Lucarelli am 15 Jul. 2021
sorry stim is stimulation that is my signal.
In my code i called it stim and not signal, but to make it easier for you i changed name in signal

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by