Filter amplitude data as a function of time, calculating the time difference

@Star Strider

Along with greet,

I wanted to ask for help, please. I need to identify the wave with the highest amplitude (maximum peak), and intersect it with a threshold value, in this case 1.5, and identify the intersection points, that is, the nodes where the maximum peak curve intersects with the threshold value In order to identify node 1 and node 2, and calculate their difference (subtract both values), I attach an illustrative image of the problem. I attach the data, and a picture of the problem. Greetings.

close all, clear all, clc
data=load('data.txt');
t=0:10:4*3600;
plot(t,data)
xlabel('Time (s)'); ylabel('Amplitude (m)')
hold on
plot(t,1.5,'r')

[Star Strider]

@ignacio bobadilla tapia — I appreciate your bringing this to my attention, however while my code worked with the original data, I could not make it work with any of the files added later, since they appear to have little in common with the first file. I added an Answer, however then deleted it for that reason.

[ignacio bobadilla tapia]

@Star Strider

First of all, thank you, a query: when obtaining a range of values, that is, a time interval, obtained from a time vs. amplitude graph, how can I extract an interval from another data series, in this case of speed in the range time already obtained previously. Thanks.

close all, clear all, clc
amplitud=load('data');
velocidad=load('vel');
tiempo=0:10:4*3600;
th=1.5; &threshold

figure, plot(t,amplitud), hold on, plot(t,th);
figure, plot(t,velocidad)

Answer 1

[Mathieu NOE]

hello Ignacio

see figure(2) results

result is displayed in command window :

time distance dt = 440.9509 s at threshold value : 1.5

code is :

clc
clearvars

y=load('data.txt');
n=length(y);
dt = 10; % 
Fs = 1/dt;
x=(1:n)*dt;

threshold = 1.5; % your value here

[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(y,x,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"

figure(1)
plot(x,y,t0_pos,s0_pos,'+r',t0_neg,s0_neg,'+g','linewidth',2,'markersize',12);grid on
legend('signal','positive slope crossing points','negative slope crossing points');
xlabel('Time (s)');

% select pos and neg crossing points closest to the max peak
[val,ind] = max(y); % max peak
x_peak = x(ind);

[minValue, closestIndex] = min(abs(t0_pos - x_peak));
closestValue_pos = t0_pos(closestIndex);

[minValue, closestIndex] = min(abs(t0_neg - x_peak));
closestValue_neg = t0_neg(closestIndex);

figure(2)
plot(x,y,x_peak,val,'dr',closestValue_pos,threshold,'+r',closestValue_neg,threshold,'+g','linewidth',2,'markersize',12);grid on
legend('signal','signal max peak','positive slope crossing points','negative slope crossing points');
xlabel('Time (s)');

dt = closestValue_neg - closestValue_pos;
disp(['time distance dt = ' num2str(dt) ' s at threshold value : ' num2str(threshold)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
% [e-mail]

% 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 ~isempty(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);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Addition:
% % Some people like to get the data points closest to the zero crossing,
% % so we return these as well
% [CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1); 
% ind2 = ind + (II-2); %update indices 
% 
% t0close = t(ind2);
% s0close = S(ind2);
end

[ignacio bobadilla tapia]

Estimated, if I wanted to implement this code for several data series, in a for loop, where it is stored in the workspace, how should I do it? Thanks.

Attached data series, there are 7 files.

All the data are in cm, you have to divide it by 100 to convert them to m, and be able to execute the code so that the dt are stored, thank you.

[Mathieu NOE]

hello again

this is a simple for loop

the main code has been converted into a function to make the whole thing easier to read

you can disable the figure plot if you don't need it;

I saw that at the beginning of some records there is a bit of noise / oscillations , so I wonder if we need to do a bit of smoothing first ? I leave it to you

so final code :

clc
clearvars

threshold = 1.5; % your value here

for ci = 1:7
    
    y=load(['data' num2str(ci) '.txt']);
    y = y/100; % convert cm to m
    dt(ci) = compute_dt(y,threshold);
    disp(['time distance dt = ' num2str(dt(ci)) ' s at threshold value : ' num2str(threshold)]);

end

function dt = compute_dt(y,threshold)

n=length(y);
dt = 10; % 
Fs = 1/dt;
x=(1:n)*dt;

% threshold = 1.5; % your value here

[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(y,x,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"

% select pos and neg crossing points closet to the max peak
[val,ind] = max(y); % max peak
x_peak = x(ind);

[minValue, closestIndex] = min(abs(t0_pos - x_peak));
closestValue_pos = t0_pos(closestIndex);

[minValue, closestIndex] = min(abs(t0_neg - x_peak));
closestValue_neg = t0_neg(closestIndex);

figure()
plot(x,y,x_peak,val,'dr',closestValue_pos,threshold,'+r',closestValue_neg,threshold,'+g','linewidth',2,'markersize',12);grid on
legend('signal','signal max peak','positive slope crossing points','negative slope crossing points');
xlabel('Time (s)');

dt = closestValue_neg - closestValue_pos;
% disp(['time distance dt = ' num2str(dt) ' s at threshold value : ' num2str(threshold)]);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
% [e-mail]

% 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 ~isempty(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);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Addition:
% % Some people like to get the data points closest to the zero crossing,
% % so we return these as well
% [CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1); 
% ind2 = ind + (II-2); %update indices 
% 
% t0close = t(ind2);
% s0close = S(ind2);
end

[ignacio bobadilla tapia]

@Mathieu NOE

Hello estimated,

My query is for example when obtaining the time nodes, that is, node 1 and node 2 of the intersection of the threshold value with the amplitude series, how can I extract in another database (speed base attached in this question) the speed range defined in the interval obtained in the time series vs. amplitude. Thanks.