
How to detect point before a positive ramp
    8 visualizaciones (últimos 30 días)
  
       Mostrar comentarios más antiguos
    
    Francesco Lucarelli
 el 15 de Jul. de 2021
  
    
    
    
    
    Comentada: Rik
      
      
 el 16 de Jul. de 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
0 comentarios
Respuesta aceptada
  Mathieu NOE
      
 el 15 de Jul. de 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 comentarios
  Mathieu NOE
      
 el 16 de Jul. de 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
      
      
 el 16 de Jul. de 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.
Ver también
Categorías
				Más información sobre Vibration Analysis en Help Center y File Exchange.
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!








