Track and extract RPM profile from vibration signal
The two-column matrix
p contains a set of points that lie
on a time-frequency ridge corresponding to a given
Each row of
p specifies one coordinate pair. If you call
rpmtrack without specifying both
p, the function opens an
interactive plot that displays the time-frequency map and enables you to select
If you have a tachometer pulse signal, use
tachorpm to extract
rpmtrack(___) with no output arguments plots
the power time-frequency map and the estimated RPM profile on an interactive
RPM Profile of Vibration Signal
Generate a vibration signal with three harmonic components. The signal is sampled at 1 kHz for 16 seconds. The signal's instantaneous frequency resembles the runup and coastdown of an engine. Compute the instantaneous phase by integrating the frequency using the trapezoidal rule.
fs = 1000; t = 0:1/fs:16; ifq = 20 + t.^6.*exp(-t); phi = 2*pi*cumtrapz(t,ifq);
The harmonic components of the signal correspond to orders 1, 2, and 3. The order-2 sinusoid has twice the amplitude of the others.
ol = [1 2 3]; amp = [5 10 5]; vib = amp*cos(ol'.*phi);
Extract and visualize the RPM profile of the signal using a point on the order-2 ridge.
time = 3; order = 2; p = [time order*ifq(t==time)]; rpmtrack(vib,fs,order,p)
RPM Profile of Revving Engine
Generate a signal that resembles the vibrations caused by revving a car engine. The signal is sampled at 1 kHz for 30 seconds and contains three harmonic components of orders 1, 2.4, and 3, with amplitudes 5, 4, and 0.5, respectively. Embed the signal in unit-variance white Gaussian noise and store it in a MATLAB® timetable. Multiply the instantaneous frequency by 60 to obtain an RPM profile. Plot the RPM profile.
fs = 1000; t = (0:1/fs:30)'; fit = @(a,x) (t-x).^6.*exp(-(t-x)).*((t-x)>=0)*a'; fis = fit([0.4 1 0.6 1],[0 6 13 17]); phi = 2*pi*cumtrapz(t,fis); ol = [1 2.4 3]; amp = [5 4 0.5]'; vib = cos(phi.*ol)*amp + randn(size(t)); xt = timetable(seconds(t),vib); plot(t,fis*60)
Derive the RPM profile from the vibration signal. Use four points at 5 second intervals to specify the ridge corresponding to order 2.4. Display a summary of the output timetable.
ndx = (5:5:20)*fs; order = ol(2); p = [t(ndx) order*fis(ndx)]; rpmest = rpmtrack(xt,order,p); summary(rpmest)
RowTimes: tout: 30001x1 duration Values: Min 0 sec Median 15 sec Max 30 sec TimeStep 0.001 sec Variables: rpm: 30001x1 double Values: Min 2.2204e-16 Median 4327.2 Max 8879.8
Plot the reconstructed RPM profile and the points used in the reconstruction.
hold on plot(seconds(rpmest.tout),rpmest.rpm,'.-') plot(t(ndx),fis(ndx)*60,'ok') hold off legend('Original','Reconstructed','Ridge points','Location','northwest')
Use the extracted RPM profile to generate the order-RPM map of the signal.
Reconstruct and plot the time-domain waveforms that compose the signal. Zoom in on a time interval occurring after the transients have decayed.
xrc = orderwaveform(vib,fs,rpmest.rpm,ol); figure plot(t,xrc) legend([repmat('Order = ',[3 1]) num2str(ol')]) xlim([5 20])
Fan Switchoff RPM Profile
Estimate the RPM profile of a fan blade as it slows down after switchoff.
An industrial roof fan spinning at 20,000 rpm is turned off. Air resistance (with a negligible contribution from bearing friction) causes the fan rotor to stop in approximately 6 seconds. A high-speed camera measures the x-coordinate of one of the fan blades at a rate of 1 kHz.
fs = 1000; t = 0:1/fs:6-1/fs; rpm0 = 20000;
Idealize the fan blade as a point mass circling the rotor center at a radius of 50 cm. The blade experiences a drag force proportional to speed, resulting in the following expression for the phase angle
where is the initial frequency and second is the decay time.
a = 0.5; f0 = rpm0/60; T = 0.75; phi = 2*pi*f0*T*(1-exp(-t/T));
Compute and plot the x- and y-coordinates of the blade. Add white Gaussian noise of variance .
x = a*cos(phi) + randn(size(phi))/10; y = a*sin(phi) + randn(size(phi))/10; plot(t,x,t,y)
rpmtrack function to determine the RPM profile. Type
at the command line to open the interactive figure.
Use the slider to adjust the frequency resolution of the time-frequency map to about 11 Hz. Assume that the signal component corresponds to order 1 and set the end time for ridge extraction to 3.0 seconds. Use the crosshair cursor in the time-frequency map and the Add button to add three points lying on the ridge. Alternatively, double-click the cursor to add the points at the locations you choose. Click Estimate to track and extract the RPM profile.
Verify that the RPM profile decays exponentially. On the Export tab, click Export and select
Generate MATLAB Script. The script appears in the Editor.
% MATLAB Code from rpmtrack GUI % Generated by MATLAB 9.12 and Signal Processing Toolbox 8.7 % Generated on 12-Oct-2021 09:36:49 % Set sample rate fs = 1000.0000; % Set order of ridge of interest order = 1.0000; % Set ridge points on ridge of interest ridgePoints = [... 0.4077 194.6023;... 0.9781 89.4886;... 2.3678 15.6250]; % Estimate RPM [rpmOut,tOut] = rpmtrack(x,fs,order,ridgePoints,... 'Method','stft',... 'FrequencyResolution',11.1612,... 'PowerPenalty',Inf,... 'FrequencyPenalty',0.0000,... 'StartTime',0.0000,... 'EndTime',3.0000);
Run the script. Display the RPM profile in a semilogarithmic plot.
semilogy(tOut,rpmOut) ylim([500 20000])
x — Input signal
Input signal, specified as a vector.
cos(pi/4*(0:159))+randn(1,160) specifies a
noisy sinusoid sampled at 2π Hz.
fs — Sample rate
positive real scalar
Sample rate, specified as a positive real scalar.
order — Ridge order
positive real scalar
Ridge order, specified as a positive real scalar.
p — Ridge points
Ridge points, specified as a two-column matrix containing one time-frequency coordinate on each row. The coordinates describe points on the time-frequency map belonging to the order ridge of interest.
xt — Input timetable
xt must contain increasing, finite,
and equally spaced row times of type
duration. The timetable
must contain only one numeric data vector with signal values.
If a timetable has missing or duplicate time points, you can fix it using the tips in Clean Timetable with Missing, Duplicate, or Nonuniform Times.
timetable(seconds(0:4)',randn(5,1)) specifies a
random variable sampled at 1 Hz for 4 seconds.
Specify optional pairs of arguments as
the argument name and
Value is the corresponding value.
Name-value arguments must appear after other arguments, but the order of the
pairs does not matter.
Before R2021a, use commas to separate each name and value, and enclose
Name in quotes.
"Method","fsst","PowerPenalty",10 specifies that the
time-frequency map is estimated using the synchrosqueezed Fourier transform,
allowing up to 10 decibels of power difference between adjacent points on a
Method — Type of time-frequency map
"stft" (default) |
Type of time-frequency map used in the estimation process, specified
"stft"— Use the short-time Fourier transform to compute a power spectrogram time-frequency map. For more details about the short-time Fourier transform, see
"fsst"— Use the synchrosqueezed Fourier transform to compute a time-frequency map. For more details about the synchrosqueezed Fourier transform, see
FrequencyResolution — Frequency resolution bandwidth
Frequency resolution bandwidth used to compute the time-frequency map, specified as a numeric scalar expressed in hertz.
PowerPenalty — Maximum difference in power between adjacent ridge points
Inf (default) | numeric scalar in dB
Maximum difference in power between adjacent ridge points, specified as a numeric scalar expressed in decibels.
Use this parameter to ensure that the ridge-extraction algorithm of
rpmtrack finds the correct ridge for the
PowerPenalty is useful when the
order ridge of interest crosses other ridges or is very close in
frequency to other ridges, but has a different power level.
FrequencyPenalty — Penalty in coarse ridge extraction
0 (default) | nonnegative scalar
Penalty in coarse ridge extraction, specified as a nonnegative scalar.
Use this parameter to ensure that the ridge-extraction algorithm of
rpmtrack avoids large jumps that could make the
ridge estimate move to an incorrect time-frequency location.
FrequencyPenalty is useful when you want to
differentiate order ridges that cross or are closely spaced in
StartTime — Start time for RPM profile estimation
input signal start time (default) | scalar in seconds |
Start time for RPM profile estimation, specified as a numeric or
EndTime — End time for RPM profile estimation
input signal end time (default) | scalar in seconds |
End time for RPM profile estimation, specified as a numeric or
rpm — Rotational speed estimate
vector | timetable
Rotational speed estimate, returned as a vector expressed in revolutions per minute.
If the input to
rpmtrack is a timetable, then
rpm is also a timetable with a single variable
rpm. The row times of the timetable are labeled
tout and are of type
tout — Time values
Time values at which the RPM profile is estimated, returned as a vector.
rpmtrack uses a two-step (coarse-fine) estimation method:
Compute a time-frequency map of
xand extract a time-frequency ridge based on a specified set of points on the ridge,
ordercorresponding to that ridge, and the optional penalty parameters
FrequencyPenalty. The extracted ridge provides a coarse estimate of the RPM profile.
Compute the order waveform corresponding to the extracted ridge using a Vold-Kalman filter and calculate a new time-frequency map from this waveform. The isolated order ridge from the new time-frequency map provides a fine estimate of the RPM profile.
 Urbanek, Jacek, Tomasz Barszcz, and Jerome Antoni. "A Two-Step Procedure for Estimation of Instantaneous Rotational Speed with Large Fluctuations." Mechanical Systems and Signal Processing. Vol. 38, 2013, pp. 96–102.
C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.
Usage notes and limitations:
Any name-value argument character or string input must be a constant at compile time.