function [turningPointTable,histogramPeaks,turningAmplitude] = ...
    turningPointHistogram(amplitude,includeEndpointsInHistogram,numberOfBins,YS,fig_num)
%TURNINGPOINTHISTOGRAM Reduce a signal to turning points and plot an
%absolute-amplitude histogram.
%
% Syntax:
%
%   turningPointHistogram(amplitude)
%   turningPointHistogram(amplitude,includeEndpointsInHistogram)
%   turningPointHistogram(amplitude,includeEndpointsInHistogram,numberOfBins)
%
%   [turningPointTable,histogramPeaks] = turningPointHistogram(...)
%
% Inputs:
%
%   amplitude
%       Numeric vector containing the signal amplitudes.
%
%   includeEndpointsInHistogram
%       Logical value specifying whether the first and last points of each
%       finite signal segment are included in the histogram.
%       Default: false
%
%   numberOfBins
%       Number of histogram bins.  Enter [] to use the square-root rule.
%       Default: []
%
% Outputs:
%
%   turningPointTable
%       Table of segment number, sample number, signed turning-point
%       amplitude, absolute amplitude, and endpoint flag.
%
%   histogramPeaks
%       Absolute amplitudes used to construct the histogram.
%
% Example:
%
%   data = readmatrix('amplitude_data.xlsx');
%   [turningPointTable,histogramPeaks] = ...
%       turningPointHistogram(data(:,2),false,20);
%
% Notes:
%
%   1. Consecutive samples with identical amplitudes are treated as a
%      plateau.  The midpoint sample of each plateau is retained.
%
%   2. NaN and Inf values are treated as gaps separating signal segments.
%
%   3. The beginning and ending points of each finite segment are retained
%      in the turning-point output, but are excluded from the histogram by
%      default because they are not necessarily local maxima or minima.

    %% Set defaults and validate inputs

    if nargin < 2 || isempty(includeEndpointsInHistogram)
        includeEndpointsInHistogram = false;
    end

    if nargin < 3
        numberOfBins = [];
    end

    validateattributes(amplitude,{'numeric'},{'vector','nonempty'}, ...
        mfilename,'amplitude',1);

    validateattributes(includeEndpointsInHistogram,{'logical','numeric'}, ...
        {'scalar'},mfilename,'includeEndpointsInHistogram',2);

    includeEndpointsInHistogram = logical(includeEndpointsInHistogram);

    if ~isempty(numberOfBins)
        validateattributes(numberOfBins,{'numeric'}, ...
            {'scalar','integer','positive','finite'}, ...
            mfilename,'numberOfBins',3);
    end

    amplitude = amplitude(:);

    %% Extract turning points

    [turningAmplitude,turningIndex,segmentNumber,isEndpoint] = ...
        extractTurningPoints(amplitude);

    if isempty(turningAmplitude)
        error('turningPointHistogram:NoFiniteData', ...
            'The amplitude vector contains no finite data.');
    end

    fprintf('Original number of samples: %d\n',numel(amplitude));
    fprintf('Number of turning points:   %d\n',numel(turningAmplitude));

    turningPointTable = table( ...
        segmentNumber,turningIndex,turningAmplitude, ...
        abs(turningAmplitude),isEndpoint, ...
        'VariableNames',{'SegmentNumber','SampleNumber', ...
        'TurningPointAmplitude','AbsoluteAmplitude','IsEndpoint'});

    %% Plot the original signal and turning points

    % figure('Name','Turning-Point Reduction','NumberTitle','off');

    %% Plot the original signal and turning points
 
    hp=figure('Name','Turning-Point Reduction','NumberTitle','off');
 
    plot(1:numel(amplitude),amplitude,'-');
    hold on
    plot(turningIndex,turningAmplitude,'o-','LineWidth',1.0,'MarkerSize',2);
    hold off
    grid on
 
    xlabel('Sample Number');
    title('Original Signal and Turning Points');
    legend('Original Signal','Turning Points','Location','best');
    ylabel(YS);

    set(hp, 'Position', [0 0 1100 450]);

    %% Select amplitudes for the histogram

    histogramMask = includeEndpointsInHistogram | ~isEndpoint;
    histogramPeaks = abs(turningAmplitude(histogramMask));

    if isempty(histogramPeaks)
        error('turningPointHistogram:NoHistogramPeaks', ...
            ['No interior turning points are available for the ' ...
             'histogram.  Set includeEndpointsInHistogram to true ' ...
             'or use a longer signal.']);
    end

    if isempty(numberOfBins)
        numberOfBins = max(5,ceil(sqrt(numel(histogramPeaks))));
    end

    %% Plot the absolute turning-point histogram

    figure('Name','Absolute Peak Histogram','NumberTitle','off');

    histogram(histogramPeaks,numberOfBins);

    grid on
    xlabel(sprintf('Absolute Turning-Point Amplitude   %s',YS));
    ylabel('Count');
    title('Absolute Peak Histogram');

end


function [turningAmplitude,turningIndex,segmentNumber,isEndpoint] = ...
    extractTurningPoints(amplitude)
%EXTRACTTURNINGPOINTS Extract turning points, one finite segment at a time.

    turningAmplitude = [];
    turningIndex     = [];
    segmentNumber    = [];
    isEndpoint       = [];

    finiteMask = isfinite(amplitude);

    if ~any(finiteMask)
        return
    end

    segmentStart = find(finiteMask & ~[false; finiteMask(1:end-1)]);
    segmentEnd   = find(finiteMask & ~[finiteMask(2:end); false]);

    numberOfSegments = numel(segmentStart);

    amplitudeCell = cell(numberOfSegments,1);
    indexCell     = cell(numberOfSegments,1);
    segmentCell   = cell(numberOfSegments,1);
    endpointCell  = cell(numberOfSegments,1);

    for k = 1:numberOfSegments

        segmentIndex = (segmentStart(k):segmentEnd(k))';

        [amplitudeCell{k},indexCell{k}] = ...
            segmentTurningPoints(amplitude(segmentIndex),segmentIndex);

        numberOfPoints = numel(amplitudeCell{k});

        segmentCell{k} = repmat(k,numberOfPoints,1);

        endpointFlag = false(numberOfPoints,1);
        endpointFlag([1 end]) = true;
        endpointCell{k} = endpointFlag;

    end

    turningAmplitude = vertcat(amplitudeCell{:});
    turningIndex     = vertcat(indexCell{:});
    segmentNumber    = vertcat(segmentCell{:});
    isEndpoint       = vertcat(endpointCell{:});

end


function [turningAmplitude,turningIndex] = segmentTurningPoints(amplitude,index)
%SEGMENTTURNINGPOINTS Turning points of one finite-data segment.
%
% Consecutive equal-amplitude samples are collapsed to the sample nearest
% the center of the plateau.  The segment endpoints are always retained.

    changeLocation = find(diff(amplitude) ~= 0);

    plateauMiddle = round( ...
        ([1; changeLocation+1] + [changeLocation; numel(amplitude)])/2);

    amplitude = amplitude(plateauMiddle);
    index     = index(plateauMiddle);

    keepMask = true(numel(amplitude),1);

    if numel(amplitude) > 2
        slopeSign = sign(diff(amplitude));
        keepMask(2:end-1) = slopeSign(1:end-1).*slopeSign(2:end) < 0;
    end

    turningAmplitude = amplitude(keepMask);
    turningIndex     = index(keepMask);

end