ACS6124 Laboratory

I: Multisensor Systems

Visakan Kadirkamanathan and Andrew Hills

2019

Aims

1. To relate theory and practice of multisensor and decision systems in application.

2. To demonstrate a comprehensive overview how different monitoring approaches are im-

plemented and designed in practice.

Objectives

At the end, the student will be able to:

1. Appreciate the design procedure of sensor signal feature extraction.

2. Extract the informative features for health monitoring application.

3. Implement the classification algorithm based decision-making system.

4. Evaluate decision system performance for design choices.

How to use MATLAB files

1. Copy all files for each lab into your working directory.

2. To get help, type help following by the name of the m-file or function. For example, if

you want to know how to use a function mean, type help mean.

1

Lab A.I — Data Acquisition and Frequency Analysis

Aims

1. To relate theory and practice of sensor signal analysis for health monitoring applications.

2. To demonstrate how different monitoring approaches are implemented and designed in

practice.

Objectives

At the end, the student will be able to

1. Appreciate the design procedure of sensor signal feature extraction for health monitoring

decision system.

2. Analyse sensor signal feature extraction design choices.

3. Use MATLAB commands to implement sensor signal analysis functions for health moni-

toring.

How to use MATLAB files

Copy all files in the directory \\uosfstore.shef.ac.uk\Shared\studentacse\public\

vk\ACS6124\FeaturesLab into your working directory in U drive and add the corresponding

working path to MATLAB.

2

Sensor Signal based Fault Diagnosis

Features are parameters derived from the measured sensor data that robustly characterise the

system it is observing, and in this laboratory, indicate the presence of faults. Objective of this

lab is to obtain a feature extractor from a set of vibration sensor signals for the purposes of fault

diagnosis of rotating machinery. The vibration signals of a machine normally carry the dynamic

information of the machine and analysing them will help us in monitoring the condition of

the machine. Here your tasks are to analyse a set of vibration sensor signals representing five

different fault types and extract features from the sensor signals. The extracted features will be

used in diagnosing of machine conditions using a pattern classification technique.

• Test Rigs

The vibration sensor data were collected from five test rigs that are used to demonstrate

effects of the following built-in faults:

Fault 1: Bearing Defect

Fault 2: Gear mesh

Fault 3: Resonance

Fault 4: Imbalance

Fault 5: Misalignment

A brief description of each fault type is given in the Appendix.

Each test rig is composed of: a rotating machinery system, a 1-axis accelerometer, a data acquisi-

tion card (PCI 1711) and a PC, as shown in Figure 1. The vibration sensor data measured by the

Figure 1: Machinery Vibration Monitoring System

1-axis accelerometer is captured by a PC-based data acquisition system via PCI 1711. The Real-

Time Windows Target (RTWT) of MATLAB is used to facilitate the communication between PC

and the A/D card.

3

Task I: Vibration Sensor Signal Analysis

The purpose of this task is to perform frequency analysis and feature extraction of the vibration

sensor measurement signals.

• Files Used

bearing.mat: contains a variable bearing which is a 50000×1 column vector of vibration

measurements of the bearing-defect rig sampled over 50 sec at a frequency of 1 kHz.

gearmesh.mat: contains a variable gearmesh which is a 50000 × 1 column vector of

vibration measurements of the gearmesh rig sampled over 50 sec at a frequency of 1 kHz.

misalignment.mat: contains a variable misalignment which is a 50000×1 column vec-

tor of vibration measurements of the misalignment rig sampled over 50 sec at a frequency

of 1 kHz.

imbalance.mat: contains a variable imbalance which is a 50000 × 1 column vector of

vibration measurements of the imbalance rig sampled over 50 sec at a frequency of 1 kHz.

resonance.mat: contains a variable resonance which is a 50000 × 1 column vector of

vibration measurements of the resonance rig sampled over 50 sec at a frequency of 1 kHz.

To analyse the signals, follow the subsequent steps:

1. Load data files of all fault types and plot them in separate figures. Consider the vibra-

tion signature in the time-domain of each fault and decide on features, which you believe

distinguish different fault types.

2. Analyse signals in the frequency domain by spectral analysis. The energy distribution of

the signal can be made by estimating its power spectral density (PSD). Here the Welchs

method will be applied using pwelch command. Then, plot the estimated PSDs in sepa-

rate figures.

% Estimate the PSD of y

>>[P,f]=pwelch(y,[],[],[],1000);%1000 is the sampling frequency (Hz)

>>plot(f,p);

In the frequency domain, decide on features, which you believe distinguish different fault types.

Task II: Signal Feature Extraction

This task allows you to extract appropriate features from the vibration data indicative of the

fault condition. From the energy distributions, it can be seen that a large variety of features

can be extracted to describe the characteristics of signals. Here we consider the energy levels

which are given by the value of root mean square (RMS) of the signal. The energy levels in three

frequency bands: 0–50 Hz, 50–200 Hz and 200–500 Hz will be included in the feature set. In

summary, four selected features are given in Table 1.

A 50000×1 sensor output vector of each fault case will be divided into (non-overlapping) 50

blocks each of length 1000. Features mentioned above are calculated based on each correspond-

ing a one-second period of the measurements and providing 4 features. Features f1 to f4 of each

block can be calculated by the following steps:

4

Feature Name Detail of the Feature

f1 Root Mean Square (RMS) of the signal

f2 RMS of low frequency: 0–50 Hz

f3 RMS of middle frequency: 50–200 Hz

f4 RMS of high frequency: 200–500 Hz

Table 1: Features to be extracted

• Normalisation: Each 1000× 1 vector is normalised for zero mean by

xnormalised = x− xmean (1)

where x is any 1000 × 1 data block. Let’s call the first 1000 normalised data x1, the next

1000 normalised data x2 and so on. Each xi vector is then used to extract the features

f1 − f4 as follows:

1. Feature f1: calculate the PSD of xi by Welch’s method. The RMS (feature f1) is calculated

by

f1 =

‖xi‖2√

N

(2)

where N is the length of the vector xi which is 1000 in this case. ‖xi‖2 is the 2-norm or the

Euclidean length of xi. The 2-norm can be obtained by using a MATLAB command norm.

2. Feature f2: signal is filtered through a low pass filter with a cut-off frequency of 50 Hz and

the RMS-value is then calculated. Type of the filter here is Butterworth, which provides a

maximally flat response in pass-band as well as in stop-band.

Create an 11th order low pass Butterworth digital filter by butter command 1 which re-

turns the filter coefficients vectors B (numerator) and A (denominator), i.e.

[B,A]=butter(11,0.1);

• Apply the filter to xi using filter MATLAB function:

y1= filter(B,A,x1); % for x1 data vector;

f2 is then obtained by calculating the RMS of the PSD of the filtered xi using (1).

3. Feature f3: Signal is filtered through a band pass filter with a cut-off frequency of 50–200 Hz

and the RMS-value is then calculated.

• Create a 13th order band pass Butterworth digital filer by butter command:

[B,A] = butter(13,[0.1 0.4]);

1The second argument 0.1 (ωn) comes from the scaled cutoff frequency 50 Hz with 1.0 corresponding to half the

sample rate, i.e. ωn = 50500 .

5

• Apply the filter to xi using filter command:

y1= filter(B,A,x1); % for x1 data vector

• f3 is then obtained by calculating the RMS of the PSD of the filtered xi using (1)

4. Feature f4: Signal is filtered through a high pass filter with a cut-off frequency of 200 Hz

and the RMS-value is then calculated.

• Create an 18th order high pass Butterworth digital filter by butter command:

[B,A] = butter(18,0.4,’high’);

• Apply the filter to xi using filter command2:

y1= filter(B,A,x1); % for x1 data vector

Task III: Data Visualisation

From Task II, you should have a set of five feature matrices each of which is composed of four

features of signal energy levels in different frequency bands. As the feature dimension is four,

it is impossible to visualise them and a dimension reduction is required. To this purpose, the

four-dimensional feature vectors will be mapped to two dimensions by principal component

analysis (PCA).

First, use load to load in the feature files of fault types 1 to 5 created in Lab. Then,combine

the feature matrices of fault cases 1–5 into a new matrix G:

G =

Fault1

Fault2

...

Fault5

(3)

The first 2 principal components of G are then calculated and used in the visualisation by the

following MATLAB commands:

c=corrcoef(G); % Calculates a correlation coefficient matrix c of G

[v,d] =eig(c); % Find the eigenvectors v and the eigenvalues d of G

T=[ v(:,end)’;v(:,end-1)’]; % Create the transformation matrix T from

% the first two principal components

z=T*G’; % Create a 2-dimensional feature vector z

plot(z(1,:),z(2,:),’o’) % Scatter plot of the 2-dimensional features

Is each fault type completely separate from the others?

2All filter orders are chosen by using buttord function which returns the lowest order digital Butterworth filter

that loses no more than 3 dB in the passband and has at least 20 dB of attenuation in the stopband.

6

Feature Name Filter Filter Order Wn

f1 Low-pass 25 Hz 7 0.05

f2 Band-pass 25 – 50 Hz 6 [0.05 0.1]

f3 Band-pass 50 – 100 Hz 9 [0.1 0.2]

f4 Band-pass 100 – 200 Hz 8 [0.2 0.4]

f5 Band-pass 200 – 350 Hz 9 [0.4 0.7]

f6 High-pass 350 Hz 16 0.7

Table 2: Features to be extracted

Task IV: Optional Additional Investigation

Now we consider the energy levels in six frequency bands: 0–25 Hz, 25–50 Hz, 50–100 Hz, 100–

200 Hz, 200–350 Hz and 350–500 Hz will be included in the feature set. Filter orders and Wn for

each frequency band can be found in Table 2. Use these values and repeat Task II and Task III.

Create and plot the 2-dimensional feature and compare with the result obtained in Task III.

7

Lab A.II — Pattern Classification

Aims

1. To relate theory and practice of decision systems for health monitoring application.

2. To demonstrate how different monitoring approaches are implemented and designed in

practice.

Objectives

At the end, the student will be able to

1. Understand the design procedure decision systems based on pattern classification ap-

proaches.

2. Evaluate decision system performance for design choices on health monitoring applica-

tion.

How to use MATLAB files

Copy all files in the directory \\uosfstore.shef.ac.uk\Shared\studentacse\public\

vk\ACS6124\ClassificationLab into your working directory in U drive and add the cor-

responding working path to MATLAB.

8

Pattern Classification

The task here is to classify machine conditions based on a simple classifier, namely the 1-nearest

neighbour classifier. The 4-dimensional feature vectors are first transformed into 2-dimensional

data to observe the data groups (this is already done in Lab A.I). The feature data are then

partitioned for training and testing. The nearest neighbour classifier is implemented and its

performance is estimated.

Training data and test data matrices:

The matrices used as the training data and the testing data will be created from the extracted

features in Lab A.I. The training data are used in construction of the classifier and the test data

is the data set used to evaluate the goodness or performance of the classifier. Here, the first 35

feature vectors of each fault case are used as the training data and the last 15 feature vectors are

then used as the test data.

Set two MATLAB variables as the training data and test data as follows:

>> trainingSet = [Fault1(1:35,:); Fault2(1:35,:); Fault3(1:35,:);

Fault4(1:35,:); Fault5(1:35,:)];

>> testingSet = [Fault1(36:end,:); Fault2(36:end,:);Fault3(36:end,:);

Fault4(36:end,:); Fault5(36:end,:)];

By using the training data as a reference set of samples with known fault types (fault classes

1–5), an unknown sample of the test data is assigned to class of its nearest sample of the training

data (see the course handout). You may also find it useful to generate training target and test

target vectors corresponding to the associated fault class labels (1 for Fault 1, 2 for Fault 2 and

so on) of the training data and the test data, respectively. These two target vectors can be used

in the nearest neighbour classifier to assign a fault class to a sample and evaluate the classifiers

performance. In this case, the training target (traintar) and the test target vectors (testtar) can be

created by the following MATLAB commands

>> trainingTarget=[ones(1,35) ones(1,35)*2 ones(1,35)*3 ones(1,35)*4

ones(1,35)*5];

>> testingTarget=[ones(1,15) ones(1,15)*2 ones(1,15)*3 ones(1,15)*4

ones(1,15)*5];

Write a program to classify fault conditions based on the 1-nearest neighbour classifier (1-NN)

by using train as a training data set and test as a test data set and evaluate accuracy of classifi-

cation (ACC) defined by

%ACC = 100× Nc

Na

(4)

where Nc is the number of all correctly classified samples and Na is the number of samples.

Hint: Nc may be calculated by using find function e.g.

Nc = length(find(group==testtar))

where group is a vector of the class labels (1–5) assigned to the test data , e.g. group=[1 1 2

2 2 2 5...] and testtar is the target vector of the test data.

9

Bonus Marks

Extend the 1-nearest neighbours algorithm developed above to create a k-nearest neighbours

solution (where k = 1, 3 and 5) and compare the classification accuracy.

Template for Nearest Neighbours algorithm

To assist you in this exercise, the template code is provided to create the nearest neighbours

algorithm. Please note that the implementation of the algorithm is not unique, so your own

solution may differ to that provided in the template code.

Distance Measure

The Euclidean distanc, d, between points a and b (in n-dimensional space) is defined as:

d =

√

(a1 − b1)2 + (a2 − b2)2 + · · ·+ (an − bn)2

Which can be simplified to just:

d =

√√√√ n∑

i=1

(ai − bi)2

Implemented in MATLAB as the function euc.m, this is:

function distance = euc(a, b)

%Euclidean Distance

% Calculates the Euclidean distance between two cases which an equal

% number of features.

%

% Author: Andrew Hills

% Date: 23/02/07

if nargin ˜= 2

error(’Two input arguments required.’);

return;

end

if ˜all(size(a) == size(b))

error(’Dimensions of inputs are not equal.’);

return;

end

if min(size(a)) ˜= 1

error(’Input is not a vector’);

return;

end

% Calculate the Euclidean Distance using the MATLAB’s norm function

distance = norm(a - b);

10

1-Nearest Neighbours Algorithm

Firstly, the variables need to be initialised:

% Specify the number of training cases:

numberOfTrainingCases = 35;

trainingSet = [Fault1(1:numberOfTrainingCases,:);

Fault2(1:numberOfTrainingCases,:);

Fault3(1:numberOfTrainingCases,:);

Fault4(1:numberOfTrainingCases,:);

Fault5(1:numberOfTrainingCases,:)];

testingSet = [Fault1(numberOfTrainingCases+1:end,:);

Fault2(numberOfTrainingCases+1:end,:);

Fault3(numberOfTrainingCases+1:end,:);

Fault4(numberOfTrainingCases+1:end,:);

Fault5(numberOfTrainingCases+1:end,:)];

Then both the training set and testing set need to be labelled. This can be easily accom-

plished by labelling Fault1 as 1, Fault2 as 2, etc. This is done for both the training and testing

sets.

% Note the below works because all faults are of equal lengths.

numberOfTestingCases = length(Fault1) - numberOfTrainingCases;

trainingTarget = [ones(1,numberOfTrainingCases),

ones(1,numberOfTrainingCases)*2,

ones(1,numberOfTrainingCases)*3,

ones(1,numberOfTrainingCases)*4,

ones(1,numberOfTrainingCases)*5];

testingTarget = [ones(1,numberOfTestingCases),

ones(1,numberOfTestingCases)*2,

ones(1,numberOfTestingCases)*3,

ones(1,numberOfTestingCases)*4,

ones(1,numberOfTestingCases)*5];

11

Using the above initial variables, the algorithm that performs the 1-nearest neighbour search

will look as follows:

% Calculate the total number of test and train classes

totalNumberOfTestingCases = numberOfTestingCases * 5;

totalNumberOfTrainingCases = numberOfTrainingCases * 5;

% Create a vector to store assigned labels

inferredLabels = zeros(1, totalNumberOfTestingCases);

% This loop cycles through each unlabelled item:

for unlabelledCaseIdx = 1:totalNumberOfTestingCases

unlabelledCase = testingSet(unlabelledCaseIdx, :);

% As any distance is shorter than infinity

shortestDistance = inf;

shortestDistanceLabel = 0; % Assign a temporary label

% This loop cycles through each labelled item:

for labelledCaseIdx = 1:totalNumberOfTrainingCases

labelledCase = trainingSet(labelledCaseIdx, :);

% Calculate the Euclidean distance:

currentDist = euc(unlabelledCase, labelledCase);

% Check the distance

if currentDist < shortestDistance

shortestDistance = currentDist;

shortestDistanceLabel = trainingTarget(labelledCaseIdx);

end

end % inner loop

% Assign the found label to the vector of inferred labels:

inferredLabels(unlabelledCaseIdx) = shortestDistanceLabel;

end % outer loop

12

Appendix — Fault Types

Figure 2: Machine components and typical fault types (reproduced from [3] )

• Bearing defect:

Bearing defect is a change in the condition of a bearing from the new state of installation

either as a result of wear of the rollers or races- or as a result of sub surface cracking,

pitting, corrosion or overheating, or cage/roller fracture of indenting due to over loading.

Damage is often associated with lack of lubricating.

• Gear damage:

Typical problems of gear boxes are: eccentric gears, gear mesh wear, backlash, broken or

chipped teeth and gear box distortion [2]. A key indicator of gear tooth wear is excitation

of the Gear Natural Frequency (GNF), which is a potential vibration frequency on any

machine that contains gears; equal to the number of teeth multiplied by the rotational

frequency of the shaft, along with sidebands around it spaced at the running speed of the

bad gear.

• Resonance:

Resonance occurs when a forcing frequency coincides with a system natural frequency,

and can cause dramatic amplitude amplification which can result in premature or even

catastrophic failure. This may be a natural frequency of the rotor but can often originate

from a support frame, foundation, gearbox or drive belts.

• Imbalance:

Imbalance in a rotor denotes that the centre of gravity and the geometric centre of a disk

are not at the same location. These two points can never be same even for a perfectly

made disk, since no material is homogeneous [2]. That is why the imbalance of the rotor

becomes one of the most vexing problems for rotating machine maintenance.

• Misalignment:

Alignment is a condition whereby machine components have the correct angular position

relative to each other; either coincident, parallel, or perpendicular, according to design

requirements. Poor alignment is one of the major causes of premature machine failure.

13

Figure 3: Angular and parallel misalignments (reproduced from [3])

References

[1] Leeds Polytechnic, An Introduction to Plant Condition Monitoring and Its Management,

1991.

[2] J.S. Rao, Vibratory Condition Monitoring of Machines, Narosa Publishing House, 2000.

[3] Machine Vibration Diagnosis, http://www.vibanalysis.co.uk/.

14