Code covered by the BSD License

# AR filter + Minimum Entropy Deconvolution for Bearing Fault Diagnosis

### Santhana Raj (view profile)

AR Filter by YuleWalker method combined with Minimum Entropy Deconvolution for bearing fault diagnos

med2d.m
```function [y_final f_final kurtIter] = med2d(x,filterSize,termIter,termDelta,plotMode)
%2D MINIMUM ENTROPY DECONVOLUTION
%  code by Geoff McDonald (glmcdona@gmail.com), February 2011
%  Used in my MSc research at the University of Alberta, Advanced
%  Control Systems Laboratory.
%
% med2d(x,filterSize,termIter,termDelta,plotMode)
%
% Algorithm Reference:
%    R.A. Wiggins, Minimum Entropy Deconvolution, Geoexploration, vol.
%    16, Elsevier Scientific Publishing, Amsterdam, 1978. pp. 2135.
%
% Inputs:
%    x:
%       Signal to perform Minimum Entropy Deconvolution on. If a single
%       column/row of data is specified, a 1d filter is designed to
%       minimize the entropy of the resulting signale. If a 2d data
%       matrix is specified, a single 1d filter will be designed to
%       minimize the averaged entropy of each column of the filtered
%       data.
%
%    filterSize:
%       This is the length of the finite inpulse filter filter to
%       design. Using a value of around 30 is appropriate depending on
%       the data. Investigate the performance difference using
%       different values.
%
%    termIter: (OPTIONAL)
%       This is the termination number of iterations. If the
%       the number of iterations exceeds this number, the MED process
%       will complete. Specify [] to use default value of 30.
%
%    termDelta: (OPTIONAL)
%       This is the termination condition. If the change in kurtosis
%       between iterations is below this threshold, the iterative
%       process will terminate. Specify [] to use the default value
%       of 0.01. You can specify a value of 0 to only terminate on
%       the termIter condition, ie. execute an exact number of
%       iterations.
%
%    plotMode:
%       If this value is > 0, plots will be generated of the iterative
%       performance and of the resulting signal.
%
% Outputs:
%    y_final:
%       The input signal(s) x, filtered by the resulting MED filter.
%       This is obtained simply as: y_final = filter(f_final,1,x);
%
%    f_final:
%       The final 1d MED filter in finite impulse response format.
%
%    kurtIter:
%       Kurtosis according to MED iteration. kurtIter(end) is the
%       final kurtosis, ie. the summed kurtosis of each y_final
%       column of y_final. sum(kurtosis(each column of y_final))
%
% Example:
%       % This will mostly extract the impulse-like
%       % disturbances caused by 0.2*(mod(n,21)==0)
%       % and plot the result.
%       n = 0:999;
%       x = [sin(n/30) + 0.2*(mod(n,21)==0);
%           sin(n/13) + 0.2*(mod(n,21)==0)];
%       [y_final f_final kurt] = med2d(x',30,[],0.01,1);
%
%
% Note:
%   The solution is not guaranteed to be the optimal solution to the
%   entropy minimizataion problem, the solution is just a local
%   minimum of the entropy and therefore a good pick.

% Assign default values for inputs
if( isempty(filterSize) )
filterSize = 30;
end
if( isempty(termIter) )
termIter = 30;
end
if( isempty(termDelta) )
termDelta = 0.01;
end
if( isempty(plotMode) )
plotMode = 0;
end

% Validate the inputs
if( sum( size(x) > 1 ) > 2 )
error('MED:InvalidInput', 'Input signal x must be of either 2d or 1d.')
elseif( sum(size(termDelta) > 1) ~= 0 || termDelta < 0 )
error('MED:InvalidInput', 'Input argument termDelta must be a positive scalar, or zero.')
elseif( sum(size(termIter) > 1) ~= 0 || mod(termIter, 1) ~= 0 || termIter <= 0 )
error('MED:InvalidInput', 'Input argument termIter must be a positive integer scalar.')
elseif(  sum(size(plotMode) > 1) ~= 0 )
error('MED:InvalidInput', 'Input argument plotMode must be a scalar.')
elseif( sum(size(filterSize) > 1) ~= 0 || filterSize <= 0 || mod(filterSize, 1) ~= 0 )
error('MED:InvalidInput', 'Input argument filterSize must be a positive integer scalar.')
end

% If the data is 1d, lets make it a column vector
if( sum(size(x)>1) == 1 )
x = x(:); % A column vector
end
L = filterSize;

% Calculate the weighted toeplitz autocorrelation matrix
% as the average autocorrelation matrix of the rows.
autoCorr = zeros(1,L);
for column = 1:size(x,2);
for k = 0:L-1
% Create the shifted x
x2 = zeros(size(x,1),1);
x2(k+1:end) = x(1:end-k,column);

% Calculate the autocorrelation at this shift
autoCorr(k+1) = autoCorr(k+1) + sum(x(:,column).*x2);
end
end
autoCorr = autoCorr / size(x,2); % Average normalization
A = toeplitz(autoCorr);
A_inv = inv(A);

% Initialize matrix sizes
f = zeros(L,1);
y = zeros(size(x,1),size(x,2));
b = zeros(L,1);
kurtIter = [];

% Assume initial filter as a delayed impulse. This decision
% was made by paper:
% H. Endo and R. Randall, Enhancement of autoregressive model based
% gear tooth fault detection technique by the use of minimum entropy
% deconvolution filter, Mechanical Systems and Signal Processing vol.21,
% no.2, February 2007
f(2) = 1;

% Iteratively adjust the filter to minimize entropy
n = 1;
while n == 1 || ( n < termIter && ( (kurt(filter(f,1,x)) - kurtIter(n-1)) > termDelta ) )

% Compute output signal
y = filter(f,1,x);

% Calculate the kurtosis
kurtIter(n) = kurt(y); %#ok<AGROW>

% Calculate the matrix g = weighted av{ crosscorr( y.^3, x) }
yc = y.^3;
weightedCrossCorr = zeros(L,1);
for column = 1:size(x,2);
for k = 0:L-1
% Create the shifted x
x2 = zeros(size(x,1),1);
x2(k+1:end) = x(1:end-k,column);

% Calculate the crosscorrelation at this shift
weightedCrossCorr(k+1) = weightedCrossCorr(k+1) + sum((y(:,column).^3).*x2);
end
end
weightedCrossCorr = weightedCrossCorr / size(x,2);

% Now we have new filter coefficients calculted as:
% f = A^-1 * g
f = A_inv*weightedCrossCorr;
f = f/sqrt(sum(f.^2)); % Normalize the filter result

% Next iteration
n = n + 1;
end

% Update the final result
f_final = f;
y_final = filter(f_final,1,x);
kurtIter(n) = kurt(y_final);

% Plot the results
if( plotMode > 0 )

figure;
subplot(2,1,1)
plot(x)
title('Input Signal(s)')
ylabel('Value')
xlabel('Sample Number')
axis tight

subplot(2,1,2)
plot(y_final)
title('Signal(s) Filtered by MED')
ylabel('Value')
xlabel('Sample Number')
axis tight

figure;
stem(f_final)
xlabel('Sample Number')
ylabel('Value')
title('Final Filter, Finite Impulse Response')

figure;
plot(kurtIter);
xlabel('MED Algorithm Iteration')
ylabel('Sum of Kurtosis for Filtered Signal(s)')

if( n == termIter )
display('Terminated for iteration condition.')
else
display('Terminated for minimum change in kurtosis condition.')
end
end
end

function [result] = kurt(x)
% This function simply calculates the summed kurtosis of the input
% signal, x.
result = mean( (sum((x-ones(size(x,1),1)*mean(x)).^4)/(size(x,1)-2))./(std(x).^4) );
end
```