function [median_val] = max8_quick_median(data_in,varargin)
%-------------------------------------------------------------------------------
% $Id: max8_quick_median.m,v 1.16 2023/07/26 21:49:41 ddelore Exp $
%
%  David S. De Lorenzo
%  12 December 2019
%
%  Quick median-value finder
%
%-------------------------------------------------------------------------------
%
%  This function estimates the median value of a 1-D real-valued
%  input array and returns a result in the same units as the input,
%  without relying on the native MATLAB call to MEDIAN().
%
%  Note, it is possible to exactly determine the median value of an
%  array via a sorting type approach, or any number of QuickSort-ing
%  type approaches.  These typically require either sorting in-place,
%  copying the input array to temporary storage and sorting on that,
%  or creating a number of sub-arrays.  Of course, it would be
%  possible to implement a binning (i.e., histogram) approach across
%  a finite number of integer-valued bins.  In any case, the purpose
%  of computing the median value for our application is as a way of
%  estimating the noise floor -- in other words, so long as our
%  estimate is "close enough" to the actual median, then the estimate
%  is satisfactory.
%
%  The very simple algorithm below implements an iterative-based
%  median estimator, basically a "median-seeking" filter operating
%  on the input data samples.  One major caveat to this iterative
%  approach is that our data samples (the FFT values) are highly
%  correlated -- in other words, the spectrum estimate for one
%  frequency bin is correlated _very_strongly_ with the estimates
%  for nearby bins.  Because of this order-dependence pathology, a
%  random permutation of the input samples is required in order to
%  avoid estimation biases.
%
%  From https://stackoverflow.com/a/7485717
%      incremental median estimator, choose eta approx. 0.001:
%        median += eta * sgn(sample - median)
%
%  From https://stackoverflow.com/a/2144754
%      it's clear eta (for the estimate of the median) has to have
%      the same units as the measurements, and so a generic value
%      like 0.001 simply doesn't make any sense. A seemingly-better
%      approach is to set eta from a running estimate of the absolute
%      deviation: for each new value sample, update
%        cumadev += abs(sample - median)  [i.e., cum. abs. deviation]
%      then set
%        eta = 1.5*cumadev/(k*k)
%      where k is the number of samples seen so far.
%
%  From https://stackoverflow.com/a/2144754
%      to modify incremental median estimator for arbitrary quantiles:
%        quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)
%      (note that for p=0.5, this reduces to the median estimator)
%
%  Therefore, I'm going to calculate the CUMULATIVE_ABSOLUTE_DEVIATION
%  for a random permutation of the input samples, then calculate an
%  update weight ETA as
%        ETA = 1.25 * CUMULATIVE_ABSOLUTE_DEVIATION / (N*N)
%
%  One of the 'N' factors in the denominator produces MAD, the mean
%  absolute deviation; then the 1.25 factor is approximately Bessel's
%  modification of sqrt(pi/2) == 1.2533, which also is the approximation
%  for standard deviation when given the mean absolute deviation [see
%  S. Gorard (2005) Revisiting a 90-year-old debate: the advantages of
%  the mean deviation - particularly the comment under Note #4, available
%  at https://www.leeds.ac.uk/educol/documents/00003759.htm].  Then, the
%  other 1/N factor in ETA can be thought of as a decreasing weight for
%  the updated estimate of median value.
%
%  Here's some test code to show median estimation convergence;
%  just expand this section, highlight/select the code (lines
%  77-207), and hit the 'F9' button to run the selection.
%{
    %% -------------------------------------------------------------------------

    fprintf('\nQUICK-MEDIAN TEST SCRIPT\n');

    clear;

    % Make some distributions
    a =  2 + ( sqrt(1)  * randn(1,1400));
    b =  8 + ( sqrt(2)  * randn(1,400));
    c = -4 + (sqrt(0.2) * randn(1,400));

    % Concatenate and randomize the order
    d = [a b c];
    d = d(randperm(length(d)));

    % Optional - just make "d" ~N(0,1)
    %{
        d =  0 + (sqrt(1)  * randn(1,1e5)); %#ok<*NASGU>
    %}

    % Optional - sort your data
    %   This is pathological, and maximizes
    %   sample-to-sample correlation
    %{
        d = sort(d,'ascend');
    %}

    % Optional - start with the maximum value
    %   This will bias the initial median estimate
    %   as large as it could go.
    %{
        [m,idx] = max(d);
        d = [m d(1:idx-1) d(idx+1:end)];
    %}

    % Optional - start with the minimum value
    %   This will bias the initial median estimate
    %   as small as it could go.
    %{
        [m,idx] = min(d);
        d = [m d(1:idx-1) d(idx+1:end)];
    %}

    % Optional - start with the median value
    %   This will reduce the initial bias.
    %{
        m = median(d);
        [~,idx] = min(abs(d-m));
        d = [m d(1:idx-1) d(idx+1:end)];
    %}

    % True mean, variance, and median values
    mean_d   = mean(d);
    var_d    = var(d);
    median_d = median(d);

    % How many points did we make
    npts = length(d);

    % Dimension storage
    median_est = NaN(1,npts);

    % Calculate the median estimate
    %   Looping this way isn't all that efficient, but we're
    %   only doing it for visualization/verification purposes.
    for ii = 1:npts
        median_est(ii) = max8_quick_median(d(1:ii));
    end

    fprintf('\n=== RESULTS ===\n\n');
    fprintf('     distribution is ~N(%.2f,%.3f)\n\n', mean_d, var_d);
    fprintf('     true median value = %6.3f <---\n', median_d);
    fprintf('      median estimate  = %6.3f\n', median_est(end));

    figure(3);
    clf(3);
    subplot(1,2,1);
    box('on');
    grid('on');
    hold('on');
    histogram(d,50,'Normalization','probability');
    line(median_d*[1 1],ylim,'Color','r','LineWidth',4);
    line(median_est(end)*[1 1],ylim,'Color','b','LineWidth',2);
    xlim([-6 14]);
    view(90,-90);
    xlabel('Data Value','FontSize',12);
    ylabel('Probability','FontSize',12);
    title('Signal Histogram','FontSize',14);
    subplot(1,2,2);
    box('on');
    grid('on');
    hold('on');
    line([0 npts],median_d*[1 1],'Color','r','LineWidth',4,...
         'DisplayName','  True Median');
    plot(1:npts,median_est,'b','LineWidth',2,...
         'DisplayName','  Median Filter');
    xlim([0 npts]);
    ylim([-6 14]);
    xlabel('Iteration Number','FontSize',12);
    ylabel('Median Estimate','FontSize',12);
    title('Convergence of Estimate','FontSize',14);
    legend('show','Location','northeast');
    drawnow();

    figure(4);
    clf(4);
    box('on');
    grid('on');
    hold('on');
    histogram(d,50,'Normalization','probability','HandleVisibility','off');
    line(median_d*[1 1],ylim,'Color','r','LineWidth',4,...
         'DisplayName','  True Median');
    line(median_est(end)*[1 1],ylim,'Color','b','LineWidth',2,...
         'DisplayName','  Median Filter');
    xlim([-6 14]);
    xlabel('Data Value','FontSize',12);
    ylabel('Probability','FontSize',12);
    title('Signal Histogram','FontSize',14);
    legend('show','location','northeast');
    drawnow();

    % Optional - plot the raw data
    %{
        figure(3);
        subplot(1,2,2);
        h = plot(1:npts,d,'.','Color',[0.8 0.8 0.8],'HandleVisibility','off');
        uistack(h,'bottom');
        drawnow();
    %}

    %% -------------------------------------------------------------------------
%}
%
%-------------------------------------------------------------------------------

% Constant determining the range of the RANDI() call
%
%   In C/C++, the source of integer randomness is the C library
%   function rand(), which we replace herein with the MATLAB randi().
%   Recall, the C/C++ rand() function returns a pseudo-random number
%   in the range of 0 to RAND_MAX, where RAND_MAX is a constant whose
%   default value may vary between implementations but is guaranteed
%   to be at least 32,767.
%
RAND_MAX = 32767;

% Provide a way for the calling routine to specify a
% particular percentile for this "median-seeking filter".
if (nargin == 1)
    % do nothing -- fall through
elseif (nargin == 2)
    % Which percentile should this filter seek
    %
    %   Note:  if your input is Gaussian, then ~68% of your input
    %   data fall within one standard deviation of the mean value,
    %   and setting your percentile value to (0.50 + 0.68/2) = 0.84
    %   will cause this function to return an estimate of the standard
    %   deviation.  You can try the following command to verify:
    %
    %       max8_quick_median(randn(1,1e4),0.84)
    %
    PERCENTILE = varargin{1};

    % Make sure the input percentile is between 0 & 1
    if (PERCENTILE >= 1) || (PERCENTILE <= 0)
        error('Your percentile must be between 0 & 1!!');
    end
else
    % Make sure your number of inputs is OK
    error('Your number of inputs is wrong!!');
end

% Make a random permutation of the input vector in order to
% eliminate possible sample-to-sample correlation.
%
% This is the Fisher-Yates Shuffle, also known as the Knuth Shuffle.
% In C/C++, the source of randomness for this shuffle is the C library
% function rand(), which we replace here with the MATLAB randi().
% Recall, the rand() function returns a pseudo-random number in the
% range of 0 to RAND_MAX, where RAND_MAX is a constant whose default
% value may vary between implementations but is guaranteed to be at
% least 32,767.
%
% Note, as a programmatic short-cut we don't bother worrying about
% modulo bias in the Fisher-Yates Shuffle.  Recall, there is a
% systematic bias in favor of "small" remainders when calculating
% the destination index for each shuffle iteration; however, this
% bias is on the order of the length of our input array divided by
% 'RAND_MAX', and for purposes of this function can be safely ignored.
%
% Note, we run this loop until i == 2 (for MATLAB 1-based indexing),
% since there's no sense in trying to permute the final entry (i=1)
% with itself.
%
% This loop replaces the following MATLAB command:
%   data_in = data_in(randperm(length(data_in)));
for ii = length(data_in):-1:2
    jj = mod(randi([0 RAND_MAX]), ii) + 1;
    temp        = data_in(ii);
    data_in(ii) = data_in(jj);
    data_in(jj) = temp;
end

% Initialize the sample counter
n = 0;

% Initialize these variables in order to avoid a Simulink
% "variable is not fully defined on some execution paths" error.
cad        = 0;
median_val = 0;

% Loop over the permuted input data.  Note, we go from the highest
% to the lowest index (not usual!!), since the random permutation
% above could have been implemented (in C/C++) with a Fisher-Yates
% Shuffle (aka a Knuth Shuffle), which traditionally increments
% from (len-1) to 0 (or 1).  Thus, for the case where we will only
% be processing part of the input vector for estimating the median
% value, there'll be no need to permute the entire input array, and
% we can terminate this loop early.
for ii = length(data_in):-1:1

    % Increment the sample counter
    n = n + 1;

    % Grab the current sample
    c = data_in(ii);

    if (n == 1)
        % Initialize values on the first iteration
        cad = 0;            % cumulative absolute deviation
        median_val = c;
    else
        % Update values for the rest of the loop
        cad = cad + abs(c - median_val);
        eta = (1.25*cad) / (n*n);
        median_val = median_val + (eta * sign(c - median_val));

        % Note: if we want to calculate a particular percentile of
        % the input, then we modify the computation as follows:
        if exist('PERCENTILE','var')
            median_val = median_val + (eta * (2*PERCENTILE - 1));
        end
    end

end

return
%===============================================================================