%% 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*

%
% 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*

% The output timetable has the same size as the input timetable.
%
% Copyright 2012 The MathWorks, Inc.