function [median_val] = max8_3pt_median(data_in)
%-------------------------------------------------------------------------------
% $Id: max8_3pt_median.m,v 1.9 2023/07/26 21:49:41 ddelore Exp $
%
%  David S. De Lorenzo
%  5 June 2020
%
%  Quick 3-pt 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,
%  by finding the middle value of sub-arrays which are selected from
%  the input 3 samples at a time, and iterating until all of the data
%  have been used up.  Each time the middle values are found, they are
%  transcribed back into the first (npts/3) entries of the input data
%  array, so this can be thought of as an "in-place" median finder,
%  which does not require allocating new memory on each loop iteration
%  or re-sizing the input array as the data are used up.  Finally, the
%  sub-arrays are not always taken from the first 'npts' entries of
%  the input, as that could (slightly) bias results by under-sampling
%  from the right-hand or "upper indexes" of the array in favor of the
%  left-hand or "lower indexes" of the array.
%
%  Here's some test code to show median estimation; just expand
%  this section, highlight/select the code (lines 30-123), and
%  hit the 'F9' button to run the selection.
%{
    %% -------------------------------------------------------------------------

    fprintf('\n3-PT 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);

    % Calculate the median estimate
    median_est = max8_3pt_median(d);

    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);
    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','k','LineWidth',2);
    xlim([-6 14]);
    view(90,-90);
    xlabel('Data Value','FontSize',12);
    ylabel('Probability','FontSize',12);
    title('Signal Histogram','FontSize',14);
    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','k','LineWidth',2,...
         'DisplayName','3-sample Mid-point');
    xlim([-6 14]);
    legend('show','location','northeast');
    drawnow();

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

% The number of points in the input array
npts = length(data_in);

% Loop until there is only one point left
while (npts > 1)
    if (npts == 2)
        % When there are only 2 points in an array,
        % then the median is defined as the mean value.
        data_in(1) = mean(data_in(1:2));
        % Update the number points for next time
        npts = 1;
    else
        % When there are 3 (or more) points, then we take
        % the median of 3-point long sub-array(s), according
        % to the following indexing rule:
        %
        %      # of samples      index range
        %           3               1 - 3       take all of the data
        %           4               1 - 3       ignore the final point
        %           5               2 - 4       ignore both end points
        %           6               1 - 6       take all of the data
        %           7               2 - 7       ignore the first point
        %           8               2 - 7       ignore both end points
        %           9               1 - 9       take all of the data
        %          10               1 - 9       ignore the final point
        %          11               2 - 10      ignore both end points
        %          12               1 - 12      take all of the data
        %          13               2 - 13      ignore the first point
        %          14               2 - 13      ignore both end points
        %         etc.               ...        continue in same pattern
        %
        % Define the index range based on the above rule
        start_idx = 1 + (mod(npts,3) == 2) + (mod(npts,6) == 1);
        stop_idx  = start_idx + 3*floor(npts/3) - 1;
        % Take the mid-points of the sub-arrays, and write
        % them back into the first part of the data array.
        data_in(1:floor(npts/3)) ...
            = median(reshape(data_in(start_idx:stop_idx),3,[]),1);
        % Update the number points for next time
        npts = floor(npts/3);
    end
end

% Transcribe to the output variable
if (npts < 1)
    % Protect against empty input
    median_val = NaN;
else
    % Return entry #1 as the median value
    median_val = data_in(1);
end

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