9 views (last 30 days)

Show older comments

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

Mathieu NOE
on 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

Rik
on 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.

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

Start Hunting!