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.
    % 
    %   Convolution adjustment:
    %    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