calculate speed and distance from gear wheel pulses sensor

Dear all,

I want to calculate de distance travelled and the instantaneous velocity of a vehicle. Trough a sensor and a dented gear wheel a pulse signal has been recorded. It is known that the vehicle travels 0.1092 m between 2 pulses.

The velocity varies considerably thus the pulses width vary acoordingly, and occasionally, the signal suffers from interference from another device, causing a temporary amplitude increase. A segment is exemplified in figure.1 below:

My approach so far has not fulfilled the analysis objective (distance and velocity calculation/plot).

I've tried to smooth and bandpass the signal as to get rid of HF and be able to identify each pulse as a peak, demonstrated in figure.2 below:

pulseSignalAvg = medfilt1(pulseSignal,60,'truncate');
flo=50;
fhi=210;
order=[48 48];
phase=0;
pulseSignalAvg_Bp = BANDPASS_fFP(pulseSignalAvg, flo, fhi, order, Fs, phase);
[pks,idx_pks]=findpeaks(pulseSignalAvg_Bp,'MinPeakProminence',0.01,'MinPeakDistance',300);

figure(3)
plot(pulseSignalAvg_Bp,'b-');
grid on
hold on
plot(idx_pks,pks,'rx') 
plot(pulseSignalAvg,'r-');
xlabel('samples');
ylabel('Amplitude');
title ('bandpass signal Peaks');
legend('Bpass pulses','peaks','smoothed pulses');

In the figure above stretch, the approach works well. I can identify the peaks and calculate the distance, and from the samples between them calculate the velocity. Though this is possible because the velocity at the exemplified stretch is fairly constant at 80km/h, as such the parameters from findpeaks 'MinPeakProminence', 'MinPeakDistance' and bandPass choosen work well. But, for stretches of different velocities results are not suitable. Figure.3 below, shows the same bandpass and findpeaks paremeters results, but for an accelerating stretch, where no longer I could get one peak per pulse.

With this approach, finding a combination of parameters that work for variable velocities seems not achievable...I guess a different approach is required to tackle this problem.

Please find figures 1, 2 and 3 attached, as well as an example pulse signal file (.mat)

I would kindly ask the community for some tips on solving this issue.

Thank you.

regards,

- Fred

Answer 1

[Mathieu NOE]

hello Frederico

Had to guess the sampling rate as it's not defined in your post... so I estimated it around 100 kHz.

I did a few things on your signal to get a "good" result as much as possible :

1/ select "shifted" data and shifted back so that the this buffer has mean value matching the full data mean value. Could be seen as a high pass effect without the phase distorsion

2/ bandpass filtered - as you already implemented - just a bit different coding.

3/ then used my zero crossing function (crossing_V7) to get the time index.

benefit vs findpeaks (or alike) : the zero crossing time values given by crossing_V7 are linearly interpolated between samples => better time accuracy vs any method that will rely on finding min or max values among sampled data.

here a plot of the result :

and the plot ot the measured time values

code below :

clc
clearvars

%%  data
load('pulseSignal.mat');
samples = length(pulseSignal);
Fs = 1e5;  % guessed !! 
dt = 1/Fs;
t = (0 : samples-1)*dt;  % 

%%% vertical shift of shifted data 
thr = max(abs(pulseSignal))/2;

ind = find(pulseSignal>thr);

mean1 = mean(pulseSignal);
mean2 = mean(pulseSignal(ind));
shift = mean1-mean2;

pulseSignal2 = pulseSignal;
pulseSignal2(ind) = pulseSignal(ind)+shift;

%%%%%%%%%%%%%%%%
NN = 2;
flo = 50;
fhi = 210;
% flo = 80;
% fhi = 200;
Wn = 2/Fs*[flo fhi];
[B,A] = butter(NN,Wn);

pulseSignals = filter(B,A,pulseSignal2);

%% zero crossing points time data (linear interpolation of time values)
threshold = 0; % your value here

[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(pulseSignals,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"


figure(1)
plot(t,pulseSignal,'b',t,pulseSignal2,'r',t,pulseSignals,'k',t0_pos1,s0_pos1,'+c','linewidth',2,'markersize',12);grid on
legend('Raw','Shifted ','Shifted & Smoothed','PS 0 crossing points');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with butterworth BP' ]);
xlabel('Time(s)');

figure(2)
plot(t0_pos1);
xlabel('Time index (sample #)');
ylabel('Time(s)');

title('plot of t0 pos1 (time values of positive slope crossing points)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 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

[Frederico Pereira]

Hi, very much appreciated for your time on answering my question.

Your suggestions are great!

Apologies I did'nt give the samplerate for the pulsesignal before. It is 65536 Hz.

Although, running your code I've noticed that in some regions (regions of low velocities), it still gives more than one "hit" per pulse. I.e., there are more than one positive zero-crossings. This happens in regions of low velocities (for examole, around 0.3s onwards), as shown in figure below:

So my calculation would count between the two hits a travelled distance of 2 pulses= 2x 0.1092 m.)

But in fact, lowering the Band Pass low cut to 10 Hz has solved this quite nicely!

The figure below shows calculated velocities ( signal detected zero crossings are shown in red +) :

Once again, thank you.

-Fred

[Mathieu NOE]

My pleasure

Glad it helped