sigqual

Signal Quality Assessment Demo
Introduction
Calculate Quality Matrix
Plot 10 Signals With Highest Kurtosis
Check Problem Data
Correlate With Found Bad Data
Look for Spikes
Check High Skewness

Signal Quality Assessment Demo

This demonstration script, which is part of DuraTools, is utilizing a vibration data base consisting of 102 vibration signals, test1 - test102. Data were recorded on different measurement points on a truck, driving by different speeds on different roads.

The task is to analyze the signals and find potentially bad measurements. We will show a few good tricks that can be made for assessing datac quality.

% Copyright (c) 2003-2006, Axiom EduTech AB, Ljusterö, Sweden. All rights reserved.
% URL: http://www.vibratools.com Email: info@vibratools.com

Introduction

With many recorded signals, it is usually necessary to perform a signal data quality check. As a first test, we make a quality matrix for data using skewness and kurtosis. We sort the matrix after kurtosis size and plot the ten "worst" signals:

CurrDir=pwd;
cd DuraData
FontSize=10;

Calculate Quality Matrix

disp('Calculating Quality Matrix with skewness and kurtosis!')
Q = qualitym('test',1,102,'wk');
% sort the matrix after kurtosis size
[y,I] = sort(Q(:,3));
Q = Q(flipud(I),:);
clc
disp('The ten signals with highest kurtosis:')
disp('  signal #   skew     kurtosis')
disp(Q(1:10,:))
disp(' ')
disp('Press any key to continue!')

Calculating Quality Matrix with skewness and kurtosis!
The ten signals with highest kurtosis:
  signal #   skew     kurtosis
   99.0000  -15.8565  942.0135
   93.0000    5.4500  225.3094
   87.0000    2.4807   56.1924
   95.0000   -0.9698   54.3978
   23.0000    0.5468   46.0364
   27.0000    0.7321   26.1947
   25.0000    0.6891   23.8230
   37.0000    0.4559   23.5546
   29.0000    0.6669   22.8375
   21.0000    0.3464   21.6435

 
Press any key to continue!

Plot 10 Signals With Highest Kurtosis

load test1
fs = headfs(Header);
for n = 1:10
    figure(n)
    eval(['load test',int2str(Q(n,1))]);
    t = maketime(Data,fs);
    plot(t,Data,'k')
    set(gca,'FontSize',FontSize)
    xlabel(headxl(Header),'FontSize',FontSize)
    ylabel(headyu(Header),'FontSize',FontSize)
    title(['Time history     ',headyl(Header)],'FontSize',FontSize)
end

Check Problem Data

There seem to be (at least) two kinds of problems; spikes, as for test99 (and more) and bad contact as for test93 and test87. The typical shape is the result of a high-pass filtered step

Correlate With Found Bad Data

First, we try to find all signals behaving as test87.

disp('Testing for test87 behaviour')
close all
load test87
t = maketime(Data,fs);
%
% Pick Out a Typical Failure Part in the Signal
x  = sigtrunc(Data,fs,113,126);
% correlate with test87 to get reference
xc = xcorr(x,Data);
scale = maximax(xc);
% correlate this with all signals, save maximax
% (maximax is max of absolute value)
for n = 1:102
    eval(['load test',int2str(n)]);
    xc = xcorr(x,Data);
    result(n,:) = [n maximax(xc)/scale];
    home
end
% sort result
%
[y,I] = sort(result(:,2));
result = result(flipud(I),:);
disp('Result of correlation')
result(1:10,:)
disp('It seems like we only have to cancel test87 and test 93!')

% we start a black list
%
blacklist = [87 93];
%
% For the spikes we look at test99 and test91.
% Look at figure 1 and 2, the signals both looks spiky!
% BUT, there is a huge difference, look at one second of
% data in figures 3 and 4
%
load test99
figure(1)
t = maketime(Data,fs);
plot(t,Data,'k')
set(gca,'FontSize',FontSize)
xlabel(headxl(Header),'FontSize',FontSize)
ylabel(headyu(Header),'FontSize',FontSize)
title(['Full Time history     ',headyl(Header)],'FontSize',FontSize)
figure(3)
plot(t,Data,'k')
axis([56.5 57.5 -250 50])
set(gca,'FontSize',FontSize)
xlabel(headxl(Header),'FontSize',FontSize)
ylabel(headyu(Header),'FontSize',FontSize)
title(['One second Time history     ',headyl(Header)],'FontSize',FontSize)
drawnow
%
load test91
figure(2)
t = maketime(Data,fs);
plot(t,Data,'k')
set(gca,'FontSize',FontSize)
xlabel(headxl(Header),'FontSize',FontSize)
ylabel(headyu(Header),'FontSize',FontSize)
title(['Full Time history     ',headyl(Header)],'FontSize',FontSize)
figure(4)
plot(t,Data,'k')
axis([100 101 -60 60])
set(gca,'FontSize',FontSize)
xlabel(headxl(Header),'FontSize',FontSize)
ylabel(headyu(Header),'FontSize',FontSize)
title(['One second Time history     ',headyl(Header)],'FontSize',FontSize)
drawnow

%

Testing for test87 behaviour
Result of correlation

ans =

   87.0000    1.0000
   93.0000    0.9127
   75.0000    0.0253
   81.0000    0.0086
   33.0000    0.0066
   35.0000    0.0064
   31.0000    0.0061
   43.0000    0.0059
   47.0000    0.0056
   45.0000    0.0055

It seems like we only have to cancel test87 and test 93!

Look for Spikes

So, the "spikes" in test91 seems not to be in error, because they are symmetric and oscillating. But the spikes in test99 look like errors! We perform a test, looking for high frequency spikes. Let's filter in a 100 Hz 1/3 octave band filter and look at the kurtosis. We then arrange the result in order of kurtosis size

disp('testing for high frequency spikes')
clear rrr
[b,a] = tersfilt(fs,100);
for n = 1:102
    if isempty(find(n == blacklist)) == 1
       eval(['load test',int2str(n)]);
       x = filter(b,a,Data);
       t = maketime(x,fs);
       [m,I] = max(abs(x));
       rrr(n,:) = [n kurtosis(x) t(I)];
       home;
   end
end
[m,I] = sort(rrr(:,2));
rrr = rrr(flipud(I),:);
disp('  signal #  kurtosis time for max')
disp(rrr(1:10,:))
disp('By inspection we find that the method points out signals with spikes')
%
% We can inspect the worst signals, and we find that the method points out
% high frequency spikes.
% We add the first 8 in the list to our black list
%
blacklist = [blacklist rrr(1:8,1)'];
blacklist = sort(blacklist);
disp('We add the bad signals to our black list')
disp(blacklist)

testing for high frequency spikes
  signal #  kurtosis time for max
   99.0000   70.8247  259.2570
   95.0000   18.4962  305.8600
   34.0000   12.5338  162.2890
   58.0000   10.5932   33.3710
   36.0000    8.7628  164.8280
   37.0000    7.9091  252.9120
   79.0000    7.9026  218.2230
   73.0000    7.6687  216.0180
   75.0000    7.6358   63.2710
   47.0000    7.0972  348.7210

By inspection we find that the method points out signals with spikes
We add the bad signals to our black list
    34    36    37    58    73    79    87    93    95    99

Check High Skewness

Some signals have a very spectacular skewness, so we resort the Q matrix for skewness

disp('we resort the quality matrix after large skewness')
[y,I] = sort(abs(Q(:,2)));
Q = Q(flipud(I),:);
disp('  signal #   skew     kurtosis')
disp(Q(1:30,:))
disp('Add the worst ones to the black list')

% Quality control can go on for ever, but here we stop, removing (rather
% arbitrarely) all signals having an absolute skewness greater than 0.3
%
% So we add to the black list
%
blacklist = [blacklist Q(find(abs(Q(:,2))>0.3),1)'];
%
% Finally, remove repetitions
%
blacklist = unique(blacklist);
disp(blacklist)
disp('QUALITY CHECK FINISHED!')
cd(CurrDir)
save blacklist blacklist

we resort the quality matrix after large skewness
  signal #   skew     kurtosis
   99.0000  -15.8565  942.0135
   93.0000    5.4500  225.3094
   87.0000    2.4807   56.1924
   81.0000    1.8802   16.3026
   75.0000    1.3780   18.1154
   95.0000   -0.9698   54.3978
   27.0000    0.7321   26.1947
   25.0000    0.6891   23.8230
   29.0000    0.6669   22.8375
   39.0000    0.6184   19.9529
   61.0000    0.5938   15.1037
   23.0000    0.5468   46.0364
   41.0000    0.5197   17.9191
   49.0000    0.5169   15.8080
   79.0000    0.5003   11.9982
   73.0000    0.4715   12.6306
   55.0000    0.4656   13.4765
   19.0000    0.4590   21.1061
   37.0000    0.4559   23.5546
   67.0000    0.4259   12.5098
  100.0000    0.3468    4.9242
   15.0000   -0.3466   14.5037
   21.0000    0.3464   21.6435
   94.0000    0.3398    4.8794
   45.0000    0.3116   17.1910
  101.0000   -0.3075   15.7117
   88.0000    0.2982    4.3745
   47.0000    0.2922   18.8049
   43.0000    0.2848   17.7868
   32.0000    0.2846    6.3054

Add the worst ones to the black list
  Columns 1 through 15 

    15    19    21    23    25    27    29    34    36    37    39    41    45    49    55

  Columns 16 through 29 

    58    61    67    73    75    79    81    87    93    94    95    99   100   101

QUALITY CHECK FINISHED!

Contents