Multipoint Optimal Minimum Entropy Deconvolution with Convolution Adjustment (MOMEDA)

Solves the exact solution to the impulse train deconvolution problem. Applies to rotating machines.

MULTIPOINT OPTIMAL MINUMUM ENTROPY DECONVOLUTION ADJUSTED
code by Geoff McDonald ([e-mail]), 2015
McDonald, Geoff L., and Qing Zhao. "Multipoint Optimal Minimum Entropy Deconvolution and Convolution Fix: Application to vibration fault detection." Mechanical Systems and Signal Processing 82 (2017): 461-477. For a PDF version of the paper, see my blog: http://www.split-code.com/rotation.html

Multipoint Optimal Minimum Entropy Deconvolution (MOMEDA) computation algorithm. This proposed
method solves the optimal solution for deconvolving a periodic train of impulses from a signal.
It is best-suited in application to rotating machine faults from vibration signals to deconvolve
the impulse-like vibration associated with many gear and bear faults.
This upload includes 'momeda' for the base deconvolution and 'momeda_spectrum' for plotting
spectrums. Generally, this spectrum will produce peaks at periods corresponding to critical
frequencies, and the resulting magnitudes may be tracked to monitor machine component health.

MULTIPOINT OPTIMAL MINUMUM ENTROPY DECONVOLUTION ADJUSTED

momeda(x,filterSize,window,period,plotMode)

Inputs:
x:
Signal to generate apply MOMEDA on. Generally this should be around the range
of 1000 to 10,000 samples covering at least 5 rotations of the elements in the machine.

filterSize:
This is the length of the finite inpulse filter filter to
design. Generally a number on the order of 500 or 1000 is good, but may
depend on the dataset length.

window:
This is the window that be convolved with the impulse train target. Generally, a
rectangular window works well, eg [1 1 1 1 1].

period:
This is the periods to test as the spectrum x-axis. It should be a decimal range, like:
range = 5:0.1:300;

plotMode:
If this value is > 0, plots will be generated of the iterative
performance and of the resulting signal.

Outputs:
MKurt:
The Multipoint Kurtosis of the Deconvolution result.

f:
Optimal FIR filter designed.

y:
Filtered output signal.

Example:

% Simple vibration fault model
close all
n = 0:4999;
h = [-0.05 0.1 -0.4 -0.8 1 -0.8 -0.4 0.1 -0.05];
faultn = 0.05*(mod(n,50)==0);
fault = filter(h,1,faultn);
noise = wgn(1,length(n),-25);
x = sin(2*pi*n/30) + 0.2*sin(2*pi*n/60) + 0.1*sin(2*pi*n/15) + fault;
xn = x + noise;

% No window. A 5-sample recangular window would be ones(5,1).
window = ones(1,1);

% 1000-sample FIR filters will be designed
L = 1000;

% Recover the fault signal of period 50
[MKurt f y] = momeda(xn,L,window,50,1);

MULTIPOINT OPTIMAL MINUMUM ENTROPY DECONVOLUTION ADJUSTED SPECTRUM

momeda_spectrum(x,filterSize,window,range,plotMode)

Inputs:
x:
Input signal.

filterSize:
FIR filter length to design.

window:
Window function.

range:
This is the periods to test as the spectrum x-axis. It should be a decimal range, like:
range = 5:0.1:300;

plotMode:
If this value is > 0 plots will be generated

Outputs:
T:
The x-axis of the resulting spectrum, representing the period in samples.

MKurt:
The y-axis of the spectrum, representing the Multipoint Kurtosis of the Deconvolution
result at the corresponding period in T.

f:
Optimal filters designed for each period in T.

y:
Outputs for each filter designe at each period in T.

T_best:
Period T corresponding to the highest MKurt.

MKurt_best:
max(MKurt)

f_best
Filter corresponding to maximum MKurt in the range provided.

y_best
Most-faulty output signal, max(MKurt).

Example:
% Simple vibration fault model
close all
n = 0:9999;
h = [-0.05 0.1 -0.4 -0.8 1 -0.8 -0.4 0.1 -0.05];
faultn = 0.05*(mod(n,50)==0);
fault = filter(h,1,faultn);
noise = wgn(1,length(n),-25);
x = sin(2*pi*n/30) + 0.2*sin(2*pi*n/60) + 0.1*sin(2*pi*n/15) + fault;
xn = x + noise;

% No window. A 5-sample recangular window would be ones(5,1)
window = ones(1,1);

% 1000-sample FIR filters will be designed
L = 1000;

% Plot a spectrum from a period of 10 to 300 (actual fault is at 50)
range = [10:0.1:300];

% Plot the spectrum
[T MKurt f y T_best MKurt_best f_best y_best] = momeda_spectrum(xn,L,window,range,1);

% Now lets extract the fault signal, assuming we know it has a period between 45 and 55
window = ones(1,1); % this is no window
range = [45:0.1:55];
[T MKurt f y T_best MKurt_best f_best y_best] = momeda_spectrum(xn,L,window,range,0);

% Plot the resulting fault signal
figure;
plot( y_best(1:1000) );
title(strcat(['Extracted fault signal (period=', num2str(T_best), ')']))