Aims
ACS6124 Laboratory Part I: Multisensor Systems
Spring 2021
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 purpose.
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
youwanttoknowhowtouseafunctionmean,typehelp 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 feature extraction for health monitoring decision sys- tem.
2. Analyse sensor signal feature extraction design choices.
3. Use MATLAB commands to implement decision system functions for health monitoring.
How to use MATLAB files
Download the files from ACS6124 Course Content at Blackboard at Labs/Multisensor Systems Lab/Data Acquisition and Frequency Analysis Lab files. Put those files in your working direc- tory and add the corresponding working path to MATLAB.
2
Spectral Analysis & Fault Diagnosis
Features are parameters derived from the measured data that robustly indicate the presence of faults. Objective of this lab is to obtain a feature extractor from a set of vibration 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 signals representing five different fault types and extract features from the signals. The extracted features will be used in diagnosing of machine conditions using a pattern classification technique.
• Test Rigs
The vibration 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 acqui- sition card (PCI 1711) and a PC, as shown in Figure 1. The vibration data measured by the 1-axis
Figure 1: Machinery Vibration Monitoring System
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 Signal Analysis
The purpose of this task is to perform frequency analysis and feature extraction of the measure- ment 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 Welch’s method will be applied using pwelch command. Then, plot the estimated PSDs in sepa- rate figures.
In the frequency domain, decide on features, which you believe distinguish different fault types.
Task II: 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 measurement 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:
% Estimate the PSD of y
>>[P,f]=pwelch(y,[],[],[],1000);%1000 is the sampling frequency (Hz)
>>plot(f,p);
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
(2)
f1 = √
N
∥Pi ∥2
where Pi is the output of the Welch’s method and N is the length of that vector. ∥Pi∥2 is the 2-norm or the Euclidean length of Pi. 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.
• Apply the filter to xi using filter MATLAB function:
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:
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 = 50 . 500
[B,A]=butter(11,0.1);
y1= filter(B,A,x1); % for x1 data vector;
[B,A] = butter(13,[0.1 0.4]);
5
• Apply the filter to xi using filter command:
• 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:
• Apply the filter to xi using filter command2:
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:
Fault1 Fault2
G= . (3) .
Fault5
The first 2 principal components of G are then calculated and used in the visualisation by the following MATLAB commands:
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.
y1= filter(B,A,x1); % for x1 data vector
[B,A] = butter(18,0.4,’high’);
y1= filter(B,A,x1); % for x1 data vector
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
6
Lab A.II — Pattern Classification Aims
1. To relate theory and practice of decision systems 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. Understand the design procedure fault diagnosis scheme based on pattern classification approaches.
2. Evaluate decision system performance for design choices on health monitoring applica- tion.
How to use MATLAB files
Download the files from ACS6124 Course Content at Blackboard at Labs/Multisensor Systems Lab/Pattern Classification Lab files. Put those files in your working directory and add the cor- responding working path to MATLAB.
7
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 classifier’s performance. In this case, the training target (traintar) and the test target vectors (testtar) can be created by the following MATLAB commands
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 (4) Na
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…]andtesttaristhetargetvectorofthetestdata.
>> 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];
8
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: 222
d= (a1 −b1) +(a2 −b2) +···+(an −bn) Which can be simplified to just:
n
d = (ai − bi)2 i=1
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);
1-Nearest Neighbours Algorithm
Firstly, the variables need to be initialised:
9
% 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];
10
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
11
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.
12
References
Figure 3: Angular and parallel misalignments (reproduced from [3])
[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/.
13