Solves the exact solution to the impulse train deconvolution problem. Applies to rotating machines.
MULTIPOINT OPTIMAL MINUMUM ENTROPY DECONVOLUTION ADJUSTED
code by Geoff McDonald (glmcdona@gmail.com), 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): 461477. For a PDF version of the paper, see my blog: http://www.splitcode.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 bestsuited in application to rotating machine faults from vibration signals to deconvolve
the impulselike 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 xaxis. 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 5sample recangular window would be ones(5,1).
window = ones(1,1);
% 1000sample 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 xaxis. 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 xaxis of the resulting spectrum, representing the period in samples.
MKurt:
The yaxis 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
Mostfaulty 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 5sample recangular window would be ones(5,1)
window = ones(1,1);
% 1000sample 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), ')']))
1.0   Added reference to the published paper. 

1.0   Added a picture 

1.0   Updated the scaling of MKurt definition. Now MKurt is generally a number between 0 and 1 as the fault indicator. 
Inspired: Minimum Entropy Deconvolution Multipack (MED, MEDA, OMEDA, MOMEDA, MCKD)