rpmtrack

Track and extract RPM profile from vibration signal

collapse all in page

Syntax

rpm = rpmtrack(x,fs,order,p)

rpm = rpmtrack(xt,order,p)

rpm = rpmtrack(___,Name=Value)

[rpm,tout] = rpmtrack(___)

rpmtrack(___)

Description

rpm = rpmtrack(x,fs,order,p) returns a time-dependent estimate of the rotational speed rpm from a vibration signal x sampled at a rate fs.

The two-column matrix p contains a set of points that lie on a time-frequency ridge corresponding to a given order. Each row of p specifies one coordinate pair. If you call rpmtrack without specifying both order and p, the function opens an interactive plot that displays the time-frequency map and enables you to select the points.

If you have a tachometer pulse signal, use tachorpm to extract rpm directly.

example

rpm = rpmtrack(xt,order,p) returns a time-dependent estimate of the rotational speed from a signal stored in the MATLAB^® timetable xt.

example

rpm = rpmtrack(___,Name=Value) specifies additional options for any of the previous syntaxes using name-value arguments. Options include the method used to estimate the time-frequency map and the starting time for the RPM profile.

[rpm,tout] = rpmtrack(___) also returns the time vector at which the RPM profile is computed.

example

rpmtrack(___) with no output arguments plots the power time-frequency map and the estimated RPM profile on an interactive figure.

example

Examples

collapse all

RPM Profile of Vibration Signal

Open Script

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

Open Live Script

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)

Figure contains an axes object. The axes object contains an object of type line.

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)

rpmest: 30001x1 timetable

Row Times:

    tout: duration

Variables:

    rpm: double

Statistics for applicable variables and row times:

            NumMissing         Min              Median             Max               Mean              Std      

    tout        0           0 sec             15 sec            30 sec            15 sec            8.6607 sec  
    rpm         0           2.2204e-16        4.3272e+03        8.8798e+03        4.2732e+03        2.5558e+03

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')

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent Original, Reconstructed, Ridge points.

Use the extracted RPM profile to generate the order-RPM map of the signal.

rpmordermap(vib,fs,rpmest.rpm)

Figure Order Map contains 5 axes objects and other objects of type uimenu, uitoolbar, uiflowcontainer. Hidden axes object 1 contains 2 objects of type line. Hidden axes object 2 contains 2 objects of type patch, text. Hidden axes object 3 contains 2 objects of type patch, text. Hidden axes object 4 contains 2 objects of type patch, text. Hidden axes object 5 contains 2 objects of type patch, text.

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])

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Order = 1, Order = 2.4, Order = 3.

Fan Switchoff RPM Profile

Open Live Script

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

$ϕ = 2 π f_{0} T (1 - e^{- t / T}),$

where $f_{0}$ is the initial frequency and $T = 0.75$ 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 $0 . 1^{2}$ .

x = a*cos(phi) + randn(size(phi))/10;
y = a*sin(phi) + randn(size(phi))/10;

plot(t,x,t,y)

Figure contains an axes object. The axes object contains 2 objects of type line.

Use the rpmtrack function to determine the RPM profile. Type

rpmtrack(x,fs)

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])

Figure contains an axes object. The axes object contains an object of type line.

Input Arguments

collapse all

`x` — Input signal
vector

Input signal, specified as a vector.

Example: cos(pi/4*(0:159))+randn(1,160) specifies a noisy sinusoid sampled at 2π Hz.

Data Types: single | double

`fs` — Sample rate
positive real scalar

Sample rate, specified as a positive real scalar.

Data Types: single | double

`order` — Ridge order
positive real scalar

Ridge order, specified as a positive real scalar.

Data Types: single | double

`p` — Ridge points
two-column matrix

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.

Data Types: single | double

`xt` — Input timetable
timetable

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.

Example: timetable(seconds(0:4)',randn(5,1)) specifies a random variable sampled at 1 Hz for 4 seconds.

Data Types: single | double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is 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.

Example: "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 ridge.

`Method` — Type of time-frequency map
`"stft"` (default) | `"fsst"`

Type of time-frequency map used in the estimation process, specified as either "stft" or "fsst".

"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 pspectrum.
"fsst" — Use the synchrosqueezed Fourier transform to compute a time-frequency map. For more details about the synchrosqueezed Fourier transform, see fsst.

`FrequencyResolution` — Frequency resolution bandwidth
numeric scalar

Frequency resolution bandwidth used to compute the time-frequency map, specified as a numeric scalar expressed in hertz.

Data Types: single | double

`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 corresponding order. 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.

Data Types: single | double

`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 frequency.

Data Types: single | double

`StartTime` — Start time for RPM profile estimation
input signal start time (default) | scalar in seconds | `duration` scalar

Start time for RPM profile estimation, specified as a numeric or duration scalar.

Data Types: single | double | duration

`EndTime` — End time for RPM profile estimation
input signal end time (default) | scalar in seconds | `duration` scalar

End time for RPM profile estimation, specified as a numeric or duration scalar.

Data Types: single | double | duration

Output Arguments

collapse all

`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 labeled rpm. The row times of the timetable are labeled tout and are of type duration.

`tout` — Time values
vector

Time values at which the RPM profile is estimated, returned as a vector.

Algorithms

rpmtrack uses a two-step (coarse-fine) estimation method:

Compute a time-frequency map of x and extract a time-frequency ridge based on a specified set of points on the ridge, p, the order corresponding to that ridge, and the optional penalty parameters PowerPenalty and 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.

References

[1] 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.

Extended Capabilities

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.

Version History

Introduced in R2018a

rpmtrack

Syntax

Description

Examples

RPM Profile of Vibration Signal

RPM Profile of Revving Engine

Fan Switchoff RPM Profile

Input Arguments

`x` — Input signal
vector

`fs` — Sample rate
positive real scalar

`order` — Ridge order
positive real scalar

`p` — Ridge points
two-column matrix

`xt` — Input timetable
timetable

Name-Value Arguments

`Method` — Type of time-frequency map
`"stft"` (default) | `"fsst"`

`FrequencyResolution` — Frequency resolution bandwidth
numeric scalar

`PowerPenalty` — Maximum difference in power between adjacent ridge points
`Inf` (default) | numeric scalar in dB

`FrequencyPenalty` — Penalty in coarse ridge extraction
`0` (default) | nonnegative scalar

`StartTime` — Start time for RPM profile estimation
input signal start time (default) | scalar in seconds | `duration` scalar

`EndTime` — End time for RPM profile estimation
input signal end time (default) | scalar in seconds | `duration` scalar

Output Arguments

`rpm` — Rotational speed estimate
vector | timetable

`tout` — Time values
vector

Algorithms

References

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

See Also

Functions

rpmtrack

Syntax

Description

Examples

RPM Profile of Vibration Signal

RPM Profile of Revving Engine

Fan Switchoff RPM Profile

Input Arguments

x — Input signal vector

fs — Sample rate positive real scalar

order — Ridge order positive real scalar

p — Ridge points two-column matrix

xt — Input timetable timetable

Name-Value Arguments

Method — Type of time-frequency map "stft" (default) | "fsst"

FrequencyResolution — Frequency resolution bandwidth numeric scalar

PowerPenalty — Maximum difference in power between adjacent ridge points Inf (default) | numeric scalar in dB

FrequencyPenalty — Penalty in coarse ridge extraction 0 (default) | nonnegative scalar

StartTime — Start time for RPM profile estimation input signal start time (default) | scalar in seconds | duration scalar

EndTime — End time for RPM profile estimation input signal end time (default) | scalar in seconds | duration scalar

Output Arguments

rpm — Rotational speed estimate vector | timetable

tout — Time values vector

Algorithms

References

Extended Capabilities

C/C++ Code Generation Generate C and C++ code using MATLAB® Coder™.

Version History

See Also

Functions

`x` — Input signal
vector

`fs` — Sample rate
positive real scalar

`order` — Ridge order
positive real scalar

`p` — Ridge points
two-column matrix

`xt` — Input timetable
timetable

`Method` — Type of time-frequency map
`"stft"` (default) | `"fsst"`

`FrequencyResolution` — Frequency resolution bandwidth
numeric scalar

`PowerPenalty` — Maximum difference in power between adjacent ridge points
`Inf` (default) | numeric scalar in dB

`FrequencyPenalty` — Penalty in coarse ridge extraction
`0` (default) | nonnegative scalar

`StartTime` — Start time for RPM profile estimation
input signal start time (default) | scalar in seconds | `duration` scalar

`EndTime` — End time for RPM profile estimation
input signal end time (default) | scalar in seconds | `duration` scalar

`rpm` — Rotational speed estimate
vector | timetable

`tout` — Time values
vector

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.