```function [y f d_norm] = omeda(x,filterSize,plotMode)
% MINIMUM ENTROPY DECONVOLUTION D-NORM EXACT SOLUTION APPROACH
%       code by Geoff McDonald (glmcdona@gmail.com), 2015
%
% omeda(x,filterSize,overlapMode,plotMode)
%  A noniterative optimal MED (OMED) computation algorithm. This method
%  solves a similar deconvolution problem to MED, and is able to directly solve for
%  the optimal filter solution as proposed by Carlos Cabrelli. This can be used is rotating
%  machine fault detection from vibration signals to detect gear and bearing faults, and OMED
%  is applied to extract the impulse-like features in the vibration.
%
%  This implementation uses the convolution adjustment proposed by myself in the second paper
%  reference, which is important to prevent this method from reaching the trivial solution of
%  deconvolving the convolution discontinuity.
%
%  Note to readers, you may want to refer to some of my other MED-based submissions:
%      MED:
%             The iterative non-optimal solution is often better for vibration fault detections. It
%             is often better than this optimal solution, since the MED problem aims to deconvolve
%             only a single-impulse as the result. As a result, OMED is able to better-extract the
%             solution of deconvolving only the single-impulse, whereas MED more commonly
%             reaches the solution of deconvolving the desired impulse train.
%      MOMEDA:
%             This is the optimal solution to the periodic impulses deconvolution problem and is recommended
%             for rotating machine faults instead of MED or OMED. Since it is non-iterative, it is able to
%             quickly generate spectrum's to diagnose machine health.
%
% Algorithm Reference:
%   Original derivation:
%    C A. Cabrelli, Minimum entropy deconvolution and simplicity: A
%    non-iterative algorithm, Geophysics, Vol. 50. No. 3. March 1984.
%
%    G.L. McDonald, Qing Zhao, Multipoint Optimal Minimum Entropy Deconvolution and Convolution
%    Fix: Application to Vibration Fault Detection, unpublished
%
% Inputs:
%    x:
%       Signal to perform Minimum Entropy Deconvolution on.
%
%    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.
%
%    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.
%
%    d_norm:
%       Final D-Norm of the filtered signal.
%
% Example:
% % Simple vibration fault model
% close all
% n = 0:1999;
% 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),-40);
% x = sin(2*pi*n/30) + 0.2*sin(2*pi*n/60) + 0.1*sin(2*pi*n/15) + fault;
% xn = x + noise;
%
% % 20-sample FIR filters will be designed
% L = 20;
%
% % Recover the fault signal
% [y, f, d_norm] = omeda(xn,L,1);
%

% Assign default values for inputs
if( isempty(filterSize) )
filterSize = 50;
end
if( isempty(plotMode) )
plotMode = 0;
end

% Validate the inputs
overlapMode = 'valid'; % This setting forces it to use the convolution adjusted solution and is highly recommended.
% I didn't provide this as an input argument since using the original 'full' mode would be
% erroneous in most cases.

if( strcmp(overlapMode,'full') == 1 )
overlap_full = 1; % not recommended
elseif( strcmp(overlapMode,'valid') == 1 )
overlap_full = 0;
else
error('OMEDA:NoOverlapMode', 'overlapMode argument must be "valid" or "full".')
end

if( sum( size(x) > 1 ) > 1 )
error('OMEDA:InvalidInput', 'Input signal x must be 1d.')
elseif(  sum(size(plotMode) > 1) ~= 0 )
error('OMEDA:InvalidInput', 'Input argument plotMode must be a scalar.')
elseif( sum(size(filterSize) > 1) ~= 0 || filterSize <= 0 || mod(filterSize, 1) ~= 0 )
error('OMEDA:InvalidInput', 'Input argument filterSize must be a positive integer scalar.')
end

L = filterSize;
x = x';

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

% Calculate X0
N = length(x);
X0 = zeros(L,N+L-1); % y = f*x where x is padded
for( l =1:L )
if( l == 1 )
X0(l,1:N) = x(1:N);
else
X0(l,2:end) = X0(l-1, 1:end-1);
end
end
if( ~overlap_full )     % "valid" region only. This fixes the convolution definition.
X0 = X0(:,L:N-1);   % y = f*x where only valid x is used
% y = X0'x to get valid output signal
end

% Calculate the inverse component
autocorr_inv = pinv(X0*X0');

% Now we need to calculate all the possible filters solutions
d_norms = zeros(size(X0,2),1); % Holds the D-Norms
k_norms = zeros(size(X0,2),1); % Holds the Kurtosis Norms

d_norm_best = -1;
f_best = zeros(filterSize,1);
j_best = 0;

% Calculate the filter solution matrix
F = autocorr_inv * X0;
Y = X0' * F;

% Find the best-result filter as the optimal answer
for j = 1:size(X0,2)
norm = dnorm(Y(:,j),j);

if( norm > d_norm_best )
% Remember this best result
j_best = j;
d_norm_best = norm;
kurt_best = kurtosis(Y(:,j));
f_best = F(:,j);
y_best = Y(:,j);
end

% Store this position result
k_norms(j) = kurtosis(Y(:,j));
d_norms(j) = norm;
end

% Select the best result as the optimal solution answer
f = f_best;
y = y_best;
d_norm = d_norm_best;

% Plot the results
if( plotMode > 0 )
figure;
subplot(5,1,1)
plot(d_norms')
hold on
stem(j_best,d_norms(j_best),'black');
hold off
title('D-Norm')

subplot(5,1,2)
plot(k_norms')
hold on
stem(j_best,k_norms(j_best),'black');
hold off
title('Kurtosis')

subplot(5,1,3)
plot(y)
title('Output')

subplot(5,1,4)
plot(x)
title('Input')

subplot(5,1,5)
stem(f)
title('Filter')
end
end

function [result] = dnorm(x,k)
result = abs(x(k))/(sum(sum(x.^2))^(1/2));
end

function [result] = kurtosis(x)
% This function simply calculates the summed kurtosis of the input
% signal, x.
result = sum(x.^4)/(sum(x.^2)^2);
end
```