![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/685123/image.png)
How to detect point before a positive ramp
14 次查看(过去 30 天)
显示 更早的评论
Hi all,
I have a stimulation signal (current mA) that is like the picture below:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/684988/image.png)
So it is made by peaks. Please find below a zoom of one pulse:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/684993/image.png)
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
0 个评论
采纳的回答
Mathieu NOE
2021-7-15
hello
I have combined threshold crossing detection with second order derivative
this is the result on a noisy pulse train ....
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/685123/image.png)
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 个评论
Mathieu NOE
2021-7-16
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
2021-7-16
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!