```function [y_final f_final kurtIter] = med2d(x,filterSize,termIter,termDelta,overlapMode,plotMode)
%2D MINIMUM ENTROPY DECONVOLUTION WITH ADJUSTED CONVOLUTION FIX
%  code by Geoff McDonald (glmcdona@gmail.com), February 2011
%  Used in my MSc research at the University of Alberta, Advanced
%  Control Systems Laboratory. This is a reference for my paper:
%      McDonald, Geoff L., Qing Zhao, and Ming J. Zuo. "Maximum correlated Kurtosis deconvolution
%      and application on gear tooth chip fault detection." Mechanical Systems and Signal Processing
%      33 (2012): 237-255.
%
%  Updated in 2015 to include the convolution adjustment fix. Proposed
%  in the second reference paper. This is important to fix MED from
%  deconvolving the trivial solution at the convolution discontinuity.
%  Using overlapMode of 'valid' uses the fix, while 'full' uses the
%  original definition used by R.A. Wiggins (not recommended).
%
%  Note to readers: This in general does not produce the optimal solution to
%  the proposed deconvolution problem. You may want to consider the other implementation
%  I've uploaded called 'omeda', which non-iteratively solves for the
%  optimal solution to a similar problem based on the Minimum Entropy Deconvolution problem proposed
%  by Carlos Cabrelli. Or if you are dealing with periodic impulses
%  I recommend looking at 'momeda' which solves the optimal solution
%  for deconvolving periodic impulses as a deconvolution goal.
%
%
% med2d(x,filterSize,termIter,termDelta,overlapMode,plotMode)
%
% Algorithm Reference:
%   Minimum Entropy Deconvolution:
%    R.A. Wiggins, Minimum Entropy Deconvolution, Geoexploration, vol.
%    16, Elsevier Scientific Publishing, Amsterdam, 1978. pp. 2135.
%
%   Convolution Adjustment Derivation:
%    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. 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.
%
%   overlapMode: (OPTIONAL)
%       You should always use 'valid' for this parameter to include the
%       convolution fix that corrects MED erroneously deconvolving
%       spurious impulses. See algorithm reference for Minimum Entropy
%       Deconvolution Adjusted Convolution (MEDA) for details on the convolution
%       adjustment. You can use 'full' if you want to reproduce the original
%       MED method by R.A Wiggins, but it is not recommended for the above reason.
%
%    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:
% % -------- 1d deconvolution example ------
% n = 0:999;
% x = [sin(n/30) + 0.2*(mod(n,21)==0)];
% % Run for 100 iterations, 30 sample FIR filter
% [y_final f_final kurt] = med2d(x',30,100,[],'valid',1);
%
% % -------- 2d deconvolution 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)];
% % Run for 100 iterations, 30 sample FIR filter
% [y_final f_final kurt] = med2d(x',30,100,[],'valid',1);
%
%
% Note:
%   The solution is not 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 = 100;
end
if( isempty(termDelta) )
termDelta = -1000;
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

% Validate the inputs
if( strcmp(overlapMode,'full') == 1 )
overlap_full = 1; % not recommended
warning('MED WARNING: Using the full option often erroneously deconvolves an impulse near the start of the signal. It is recommended to use the valid option instead')
elseif( strcmp(overlapMode,'valid') == 1 )
overlap_full = 0;
else
error('MEDD:NoOverlapMode', 'overlapMode argument must be "valid" or "full".')
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;

%%% First we need to calculate the autocorrelation matrix R
N = length(x);
Xm0 = zeros(L,N+L-1); % y = f*x where x is padded

for( l =1:L )
if( l == 1 )
Xm0(l,1:N) = x(1:N);
else
Xm0(l,2:end) = Xm0(l-1, 1:end-1);
end
end

if( ~overlap_full )     % "valid" region only
Xm0 = Xm0(:,L:N-1); % y = f*x where only valid x is used
% y = Xm0'x to get valid output signal
end

average_autocorr = Xm0*Xm0';

% Generate the toeplitz matrix R and it's inverse
R_inv = pinv(average_autocorr);

% 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 difference filter
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 = Xm0'*f;

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

yc = y.^3;
weightedCrossCorr = Xm0*yc;

% Now we have new filter coefficients calculted as:
% f = A^-1 * g
f = R_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 = Xm0'*f_final;
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) );
%result = sum(x.^4)/(sum(x.^2)^2);
end```