Main Content

envspectrum

Envelope spectrum for machinery diagnosis

Description

example

es = envspectrum(x,fs) returns the envelope spectrum of a signal x sampled at a rate fs. If x is a matrix, then the function computes the envelope spectrum independently for each column and returns the result in the corresponding column of es.

example

es = envspectrum(xt) returns the envelope spectrum of a signal stored in the MATLAB® timetable xt.

example

es = envspectrum(___,Name=Value) specifies additional options for any of the previous syntaxes using name-value arguments. Options include the algorithm used to compute the envelope signal and the frequency band over which to estimate the spectrum.

[es,f,env,t] = envspectrum(___) returns f, a vector of frequencies at which es is computed; env, the envelope signal; and t, the times at which env is computed.

envspectrum(___) with no output arguments plots the envelope signal and the envelope spectrum in the current figure.

Examples

collapse all

Simulate two vibration signals, one from a healthy bearing and one from a damaged bearing. Compute and compare their envelope spectra.

A bearing with a pitch diameter of 12 cm has eight rolling elements. Each rolling element has a diameter of 2 cm. The outer race remains stationary as the inner race is driven at 25 cycles per second. An accelerometer samples the bearing vibrations at 10 kHz.

fs = 10000;
f0 = 25;
n = 8;
d = 0.02;
p = 0.12;

The vibration signal from the healthy bearing includes several orders of the driving frequency. Plot 0.1 second of data.

t = 0:1/fs:1-1/fs;
z = [1 0.5 0.2 0.1 0.05]*sin(2*pi*f0*[1 2 3 4 5]'.*t)/5;

plot(t,z)
xlim([0.4 0.5])

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

A defect in the outer race of the bearing causes a series of 5 millisecond impacts on the bearing. Eventually, those impacts result in bearing wear. The impacts occur at the ball pass frequency outer race (BPFO) of the bearing,

BPFO=12nf0[1-dpcosθ],

where f0 is the driving rate, n is the number of rolling elements, d is the diameter of the rolling elements, p is the pitch diameter of the bearing, and θ is the bearing contact angle. Assume a contact angle of zero and compute the BPFO.

ca = 0;
bpfo = n*f0/2*(1-d/p*cos(ca))
bpfo = 83.3333

Model each impact as a 3 kHz sinusoid windowed by a flat top window. Make the impact periodic by convolving it with a comb function. Plot 0.1 second of data.

fImpact = 3000;
tImpact = 0:1/fs:5e-3-1/fs;
xImpact = sin(2*pi*fImpact*tImpact).*flattopwin(length(tImpact))'/10;

xComb = zeros(size(t));
xComb(1:fs/bpfo:end) = 1;

x = conv(xComb,xImpact,'same')/3;

plot(t,x+z)
xlim([0.4 0.5])

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

Add white Gaussian noise to the signals. Specify a noise variance of 1/30². Plot 0.1 second of data.

yGood = z + randn(size(z))/30;
yBad = x+z + randn(size(z))/30;
plot(t,yGood,t,yBad)
xlim([0.4 0.5])
legend('Healthy','Damaged')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Healthy, Damaged.

Compute and plot the envelope signals and spectra.

envspectrum([yGood' yBad'],fs)
xlim([0 10*bpfo]/1000)

Figure contains 2 axes objects. Axes object 1 with title Envelope Signal, xlabel Time (ms), ylabel Amplitude contains 2 objects of type line. Axes object 2 with title Envelope Spectrum, xlabel Frequency (kHz), ylabel Peak Amplitude contains 2 objects of type line.

Compare the peak locations to the frequencies of harmonics of the BPFO. The BPFO harmonics in the envelope spectrum are a sign of bearing wear.

harmImpact = (1:10)*bpfo;
[X,Y] = meshgrid(harmImpact,ylim);

hold on
plot(X/1000,Y,':k')
legend('Healthy','Damaged','BPFO harmonics')
hold off

Figure contains 2 axes objects. Axes object 1 with title Envelope Signal, xlabel Time (ms), ylabel Amplitude contains 2 objects of type line. Axes object 2 with title Envelope Spectrum, xlabel Frequency (kHz), ylabel Peak Amplitude contains 12 objects of type line. These objects represent Healthy, Damaged, BPFO harmonics.

Compute the Welch spectra of the signals. Specify a frequency resolution of 5 Hz.

figure
pspectrum([yGood' yBad'],fs,'FrequencyResolution',5)
legend('Healthy','Damaged')

Figure contains an axes object. The axes object with title Fres = 5 Hz, xlabel Frequency (kHz), ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Healthy, Damaged.

At the low end of the spectrum, the driving frequency and its orders obscure other features. The spectrum of the healthy bearing and the spectrum of the damaged bearing are indistinguishable.

xlim([0 10*bpfo]/1000)

Figure contains an axes object. The axes object with title Fres = 5 Hz, xlabel Frequency (kHz), ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Healthy, Damaged.

The spectrum of the faulty bearing shows BPFO harmonics modulated by the impact frequency.

xlim((bpfo*[-10 10]+fImpact)/1000)

Figure contains an axes object. The axes object with title Fres = 5 Hz, xlabel Frequency (kHz), ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Healthy, Damaged.

Generate a two-channel signal that resembles the vibration signals from a bearing that completes a rotation every 10 milliseconds. The signal is sampled at 10 kHz for 0.2 seconds, which corresponds to 20 bearing rotations.

fs = 10000;
tmax = 20;
mlt = 0.01;
t = 0:1/fs:mlt-1/fs;

During each 10-millisecond interval:

  • The first channel is a damped sinusoid with damping constant 700 and sinusoid frequency 600 Hz.

  • The second channel is another damped sinusoid with damping constant 800 and sinusoid frequency 500 Hz. The second channel lags the first channel by 5 milliseconds.

Plot the signal.

y1 = sin(2*pi*600*t).*exp(-700*t);
y2 = sin(2*pi*500*t).*exp(-800*t);
y2 = [y2(51:100) y2(1:50)];

T = (0:1/fs:mlt*tmax-1/fs)';
Y = repmat([y1;y2],1,tmax)';

plot(T,Y)

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

Create a duration array using the time interval T. Construct a timetable with the duration array and the two-channel signal.

dt = seconds(T);
ttb = timetable(dt,Y);

Use envspectrum with no output arguments to display the envelope signal and envelope spectrum of the two channels. Compute the spectrum on the whole Nyquist interval, excluding 100 Hz intervals at the ends.

envspectrum(ttb,'Band',[100 4900])

Figure contains 2 axes objects. Axes object 1 with title Envelope Signal, xlabel Time, ylabel Amplitude contains 2 objects of type line. Axes object 2 with title Envelope Spectrum, xlabel Frequency (kHz), ylabel Peak Amplitude contains 2 objects of type line.

The envelope spectra of the signals have peaks at integer multiples of the repetition rate of 1/0.01 = 0.1 kHz. This is just as expected. envspectrum removes the high-frequency sinusoidal components and focuses on the lower-frequency repetition behavior. This is why the envelope spectrum is a useful tool for the analysis of rotational machinery.

Compute the envelope signal and the times at which it is computed. Check the types of the output variables.

[~,~,ttbenv,ttbt] = envspectrum(ttb,'Band',[100 4900]);
whos ttb*
  Name           Size            Bytes  Class        Attributes

  ttb         2000x1             48977  timetable              
  ttbenv      2000x1             48985  timetable              
  ttbt        2000x1             16002  duration               

The time vector is of duration type, like the time values of the input timetable. The output timetable has the same size as the input timetable.

Store each channel of the input timetable as a separate variable. Compute the envelope signal and the time vector. Check the output types.

btb = timetable(dt,Y(:,1),Y(:,2));

[~,~,btbenv,btbt] = envspectrum(btb,'Band',[100 4900]);
whos btb*
  Name           Size            Bytes  Class        Attributes

  btb         2000x2             49199  timetable              
  btbenv      2000x2             49219  timetable              
  btbt        2000x1             16002  duration               

The output timetable has the same size as the input timetable.

Generate a signal sampled at 1 kHz for 5 seconds. The signal consists of 0.01-second rectangular pulses that repeat every T = 0.25 second. Amplitude modulate the signal onto a sinusoid of carrier frequency 150 Hz.

fs = 1e3;
tmax = 5;

t = 0:1/fs:tmax;
y = pulstran(t,0:0.25:tmax,rectpuls=0.01);

fc = 150;
z = modulate(y,fc,fs);

Plot the original and modulated signals. Show only the first few cycles.

plot(t,y,t,z,"-")
grid on
axis([0 1 -1.1 1.1])

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

Compute the envelope and envelope spectrum of the signal. Determine the signal envelope using complex demodulation. Compute the envelope spectrum on a 20 Hz interval centered at the carrier frequency.

[q,f,e,te] = envspectrum(z,fs,Method="demod",Band=[fc-10 fc+10]);

Plot the envelope signal and the envelope spectrum. Zoom in on the interval from 0 to 50 Hz.

subplot(2,1,1)
plot(te,e)
xlabel("Time")
title("Envelope")

subplot(2,1,2)
plot(f,q)
xlim([0 50])
xlabel("Frequency")
title("Envelope Spectrum")

Figure contains 2 axes objects. Axes object 1 with title Envelope, xlabel Time contains an object of type line. Axes object 2 with title Envelope Spectrum, xlabel Frequency contains an object of type line.

The envelope signal has the same period in time, T = 0.25 second, as the original signal. The envelope spectrum has pulses at 1 / T = 4 Hz.

Repeat the computation, but now use the hilbert function to compute the envelope. Bandpass-filter the signal using a 10th-order finite impulse response (FIR) filter. Plot the envelope signal and envelope spectrum using the built-in functionality of envspectrum.

envspectrum(z,fs,Method="hilbert",FilterOrder=10)

Figure contains 2 axes objects. Axes object 1 with title Envelope Signal, xlabel Time (secs), ylabel Amplitude contains an object of type line. Axes object 2 with title Envelope Spectrum, xlabel Frequency (Hz), ylabel Peak Amplitude contains an object of type line.

Embed the signal in white Gaussian noise of variance 1/3. Plot the result.

zn = z + randn(size(z))/3;

plot(t,zn,"-")
grid on
axis([0 1 -1.1 1.1])

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

Compute and display the envelope signal and envelope spectrum. Compute the envelope spectrum using complex demodulation on a 10 Hz interval centered at the carrier frequency. Zoom in on the interval from 0 to 50 Hz.

envspectrum(zn,fs,Band=[fc-5 fc+5])
xlim([0 50])

Figure contains 2 axes objects. Axes object 1 with title Envelope Signal, xlabel Time (secs), ylabel Amplitude contains an object of type line. Axes object 2 with title Envelope Spectrum, xlabel Frequency (Hz), ylabel Peak Amplitude contains an object of type line.

Input Arguments

collapse all

Input signal, specified as a vector or a matrix. If x is a vector, it is treated as a single channel. If x is a matrix, then envspectrum computes the envelope spectrum independently for each column and returns the result in the corresponding column of es.

Example: cos(pi/4*(0:159))+randn(1,160) is a single-channel row-vector signal.

Example: cos(pi./[4;2]*(0:159))'+randn(160,2) is a two-channel signal.

Data Types: single | double
Complex Number Support: Yes

Sample rate, specified as a positive real scalar.

Data Types: single | double

Input timetable. xt must contain increasing finite row times. If xt represents a multichannel signal, then it must have either a single variable containing a matrix or multiple variables consisting of vectors.

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,2)) specifies a two-channel, 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','hilbert','FilterOrder',30,'Band',[0 fs/4] computes the envelope spectrum between 0 and one-half the Nyquist frequency using a 30th-order bandpass filter and computing the envelope of the analytic signal.

Algorithm for computing the envelope signal, specified as either "hilbert" or "demod". For more information, see Algorithms.

Frequency band to compute envelope spectrum, specified as a two-element vector of strictly increasing values between 0 and the Nyquist frequency.

Data Types: single | double
Complex Number Support: Yes

FIR filter order, specified as a positive integer scalar.

  • If Method is "hilbert", then this argument specifies the order of an FIR bandpass filter.

  • If Method is "demod", then this argument specifies the order of an FIR lowpass filter.

Data Types: single | double

Output Arguments

collapse all

Envelope spectrum, returned as a vector or matrix.

Frequencies at which the envelope spectrum is computed, returned as a vector.

Envelope signal, returned as a vector, matrix, or timetable.

If the input to envspectrum is a timetable, then env is also a timetable. The time values of env have the same format as the time values of the input timetable.

  • If the input is a timetable with a single variable containing a matrix, then env has a single variable containing a matrix.

  • If the input is a timetable with multiple variables consisting of vectors, then env has multiple variables consisting of vectors.

Time values at which the envelope signal is computed, returned as a vector.

If the input to envspectrum is a timetable, then t has the same format as the time values of the input timetable.

Algorithms

envspectrum initially removes the DC bias from the input signal, x, and then computes the envelope signal.

  • If Method is set to "hilbert", the function:

    1. Bandpass-filters the signal. The FIR filter has an order specified by FilterOrder and cutoff frequencies at ba(1) and ba(2), where ba is a frequency band specified using Band.

    2. Computes the analytic signal using the hilbert function.

    3. Computes the envelope signal as the absolute value of the analytic signal.

  • If Method is set to "demod", the function:

    1. Performs complex demodulation of the signal. The signal is multiplied by exp(j2πf0t), where f0 = (ba(1) + ba(2))/2.

    2. Lowpass-filters the demodulated signal to compute the analytic signal. The FIR filter has an order specified by FilterOrder and a cutoff frequency of (ba(2)ba(1))/2.

    3. Computes the envelope signal as twice the absolute value of the analytic signal.

After computing the envelope signal, the function removes the DC bias from the envelope and computes the envelope spectrum using the FFT.

References

[1] Randall, Robert Bond. Vibration-Based Condition Monitoring. Chichester, UK: John Wiley & Sons, 2011.

Extended Capabilities

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

Version History

Introduced in R2017b