envspectrum

Envelope spectrum for machinery diagnosis

Syntax

es = envspectrum(x,fs)

es = envspectrum(xt)

es = envspectrum(___,Name=Value)

[es,f,env,t] = envspectrum(___)

envspectrum(___)

Description

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.

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

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

Envelope Spectrum of Vibration Signals

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

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=\frac{1}{2}n{f}_0[1-\frac{d}{p}\cos \theta ],$$

where $$f_0$$ 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 the impacts as a periodic train of 3 kHz sinusoids windowed by a flat top window. 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;

x = pulstran(t,0:1/bpfo:t(end),xImpact,fs)/3;

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

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

Compute and plot the envelope signals and spectra.

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

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;
xline(harmImpact/1000,":")
legend("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")

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)

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

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

Envelope Spectrum of Timetable

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)

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

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             49231  timetable              
  btbenv      2000x2             49251  timetable              
  btbt        2000x1             16002  duration               

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

Envelope Spectrum of Modulated Pulses

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

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

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)

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

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

Input Arguments

x — Input signal

vector | matrix

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

fs — Sample rate

positive real scalar

Sample rate, specified as a positive real scalar.

Data Types: single | double

xt — Input timetable

timetable

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.

Method — Algorithm for computing envelope signal

"demod" (default) | "hilbert"

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

Band — Frequency band to compute envelope spectrum

[fs/4 fs*3/8] (default) | two-element vector

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

FilterOrder — FIR filter order

50 (default) | positive integer scalar

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

es — Envelope spectrum

vector | matrix

Envelope spectrum, returned as a vector or matrix.

f — Frequencies

vector

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

env — Envelope signal

vector | matrix | timetable

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.

t — Time values

vector

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πf₀t), where f₀ = (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