MATLAB Answers

calculate speed and distance from gear wheel pulses sensor

조회 수: 1(최근 30일)
Frederico Pereira
Frederico Pereira 2021년 7월 12일
댓글: Mathieu NOE 2021년 7월 15일
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

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 7월 13일
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
% 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
  댓글 수: 2
Mathieu NOE
Mathieu NOE 2021년 7월 15일
My pleasure
Glad it helped

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by