%% Compute Hilbert Spectrum of Vibration Signal
% Simulate a vibration signal from a damaged bearing. Compute the Hilbert spectrum of this signal and look for defects.
% 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.

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

% A resonance is excited in the bearing vibration halfway through the measurement process.

yHealthy = (1+1./(1+linspace(-10,10,length(yHealthy)).^4)).*yHealthy;

%
% The resonance introduces a defect in the outer race of the bearing that results in progressive wear. The defect causes a series of impacts that recur at the ball pass frequency outer race (BPFO) of the bearing:
%  {\rm BPFO} = {1\over2}nf_0\left[1-{d\over p}\cos\theta\right],
% 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 \theta is the bearing contact angle. Assume a contact angle of 15° and compute the BPFO.

ca = 15;
bpfo = n*f0/2*(1-d/p*cosd(ca));

%
% Use the pulstran function to model the impacts as a periodic train of 5-millisecond sinusoids. Each 3 kHz sinusoid is windowed by a flat top window. Use a power law to introduce progressive wear in the bearing vibration signal.

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

tx = 0:1/bpfo:t(end);
tx = [tx; 1.3.^tx-2];

nWear = 49000;
nSamples = 100000;
yImpact = pulstran(t,tx',xImpact,fs)/5;
yImpact = [zeros(1,nWear) yImpact(1,(nWear+1):nSamples)];

%
% Generate the BPFO vibration signal by adding the impacts to the healthy bearing signal. Plot the signal and select a 0.3-second interval starting at 5.0 seconds.

yBPFO = yImpact + yHealthy;

xLimLeft = 5.0;
xLimRight = 5.3;
yMin = -0.6;
yMax = 0.6;

plot(t,yBPFO)

hold on
[limLeft,limRight] = meshgrid([xLimLeft xLimRight],[yMin yMax]);
plot(limLeft,limRight,'--')
hold off

%
% Zoom in on the selected interval to visualize the effect of the impacts.

xlim([xLimLeft xLimRight])

%
% Add white Gaussian noise to the signals. Specify a noise variance of 1/150^2.

rn = 150;
yGood = yHealthy + randn(size(yHealthy))/rn;
yBad = yBPFO + randn(size(yHealthy))/rn;

plot(t,yGood,t,yBad)
xlim([xLimLeft xLimRight])
legend('Healthy','Damaged')

%
% Use emd to perform an empirical mode decomposition of the healthy bearing signal. Compute the first five intrinsic mode functions (IMFs). Use the 'Display' name-value argument to output a table showing the number of sifting iterations, the relative tolerance, and the sifting stop criterion for each IMF.

imfGood = emd(yGood,'MaxNumIMF',5,'Display',1);

% Use emd without output arguments to visualize the first three IMFs and the residual.

emd(yGood,'MaxNumIMF',5)

%
% Compute and visualize the IMFs of the defective bearing signal. The first empirical mode reveals the high-frequency impacts. This high-frequency mode increases in energy as the wear progresses.

imfBad = emd(yBad,'MaxNumIMF',5,'Display',1);
emd(yBad,'MaxNumIMF',5)

%
% Plot the Hilbert spectrum of the first empirical mode of the defective bearing signal. The first mode captures the effect of high-frequency impacts. The energy of the impacts increases as the bearing wear progresses.

figure
hht(imfBad(:,1),fs)

% The Hilbert spectrum of the third mode shows the resonance in the vibration signal. Restrict the frequency range from 0 Hz to 100 Hz.

hht(imfBad(:,3),fs,'FrequencyLimits',[0 100])

% For comparison, plot the Hilbert spectra of the first and third modes of the healthy bearing signal.

subplot(2,1,1)
hht(imfGood(:,1),fs)
subplot(2,1,2)
hht(imfGood(:,3),fs,'FrequencyLimits',[0 100])

% Copyright 2019 The MathWorks, Inc.