function [spectral_peaks] ...
             = max8_find_peaks(S, T, fft_stats, MAXPKS, SPECTRAL_PEAKS, ...
                               FFT_WIN_FUN_NOISE_BW, PEAK_DETECTION_THRESHOLD)
%-------------------------------------------------------------------------------
% $Id: max8_find_peaks.m,v 1.8 2023/07/26 21:49:41 ddelore Exp $
%
%  David S. De Lorenzo
%  12 December 2019
%
%  Find spectral peaks
%
%-------------------------------------------------------------------------------
%
%  Search for spectral peaks, one at a time.  This function finds
%  statistics for the largest peaks in an array.
%
%  Note, to estimate the power of a CW (single-bin) feature:
%
%                spectrum   the number of data samples
%                 height    which were sent to the FFT
%                    |           |
%                    v           v
%       power = (amplitude / n_samples)^2;
%
%  And this estimates the power of a multi-bin feature, correcting for
%  rounding of the spectrum caused by the time-domain window function:
%
%       idx_start = floor(idx + delta_idx - w/2);
%       idx_stop  = ceil( idx + delta_idx + w/2);
%       amplitude_vector = shift_T(idx_start:idx_stop);
%       power = sum(amplitude_vector.^2) ...
%               / (2*sqrt(WIN_FUN_NOISE_BW)*n_samples^2);
%
%  TBD:
%        1. add wraparound processing since the right end of the FFT
%           spectrum actually connects back around to the left end.
%        2. correct for non-flat peaks that have a frequency estimate
%           that is biased away from mid-way between either roll-off.
%
%-------------------------------------------------------------------------------

% Initialize the output struct
spectral_peaks = repmat(SPECTRAL_PEAKS, 1, MAXPKS);

% Copy variables to local storage
fs           = fft_stats.fs;
fi           = fft_stats.fi;
nfft         = fft_stats.nfft;
median_value = fft_stats.median_value;
noise_pwr    = fft_stats.noise_pwr;

% Consistency check
if ((nfft ~= length(T)) || (~isempty(S) && (nfft ~= length(S))))
    error('You have a problem with your FFT data!!');
end

% Shift zero-frequency component to center of spectrum
%   'shift_S' is complex {1 x nfft}
%   'shift_T' is  real   {1 x nfft}
shift_S = fftshift(S);
shift_T = fftshift(T);

% Initialize the number of spectral peaks
peak_ctr = 0;

% Find peaks in the current spectrum
%
%   We'll conduct the search (2 x MAXPKS) times since we want to
%   make sure that short-lived "tall & skinny" peaks which fall
%   below our J/N detection threshold do not over-ride detection
%   of longer-lived peaks.  If/when that happens, it could cause
%   the peaks-to-states association to fail, leading to RFI state(s)
%   being incorrectly deleted (which would have follow-on impacts to
%   notch filter initialization).
%
%   Note, the loop variable 'ff' is only a counter; it is not used
%   for indexing.
%
for ff = 1:2*MAXPKS

    % Look for the highest (remaining) peak using native FINDPEAKS()
    %{
        %
        %   We don't grab the peak width return argument since we'll
        %   calculate that below based on start & stop frequencies.
        %      |
        %      |   /--- We also don't grab the "prominence" since it'll
        %      |  /       not get used in the search algorithm below.
        %      v v
        [~,idx,~,~] = findpeaks(shift_T, ...
                                 'SortStr',           'descend',  ...
                                 'NPeaks',               1,       ...
                                 'MinPeakHeight',      -inf,      ...
                                 'MinPeakProminence',    0,       ...
                                 'WidthReference',    'halfprom', ...
                                 'Annotate',          'extents');

        % Check FINDPEAKS() output
        if isempty(idx)
            % Terminate the loop
            break;
        end
    %}

    % Look for the highest (remaining) peak using a simple call to MAX()
    % {
        [~,idx] = max(shift_T);
    %}

    % The peak is at one end of the vector
    if ((idx == 1) || (idx == nfft))

        % Set the index correction to zero
        delta_idx = 0;

        % This is the amplitude estimate
        ampl = shift_T(idx);

        % This is the phase estimate
        if ~isempty(shift_S)
            cc = shift_S(idx);                  % complex valued!!
            phi = rad2deg(angle(cc));           % units of degrees!!
        else
            phi = NaN;                          % no complex spectrum
        end

    % The peak is somewhere in the middle
    else

        % How far off the index is the actual peak
        %   (parabolic interpolation)
        %
        %    Note, strictly speaking, 'delta_idx' is ~(-0.5,+0.5);
        %    however, due to limited machine precision, it could
        %    be ~[-0.5,+0.5].
        %
        delta_idx = 1/2 * (shift_T(idx+1) - shift_T(idx-1)) ...
                    / (2*shift_T(idx) - shift_T(idx+1) - shift_T(idx-1));

        % This is the amplitude estimate
        %   (quadratic interpolation)
        ampl = shift_T(idx) - (delta_idx/4) * (shift_T(idx-1) - shift_T(idx+1));

        % This is the phase estimate
        %   Linear interpolation
        if ~isempty(shift_S)
            if (delta_idx >= 0)
                aa = shift_S(idx);
                bb = shift_S(idx+1);
                cc = aa + delta_idx*(bb-aa);    % complex valued!!
            else
                aa = shift_S(idx);
                bb = shift_S(idx-1);
                cc = aa - delta_idx*(bb-aa);    % complex valued!!
            end
            phi = rad2deg(angle(cc));           % units of degrees!!
        else
            phi = NaN;                          % no complex spectrum
        end

    end % if/else on peak index

    % This is the frequency estimate
    %
    %   Convert location from bins to frequency
    %
    freq = (-nfft/2 + idx + delta_idx - 1)*(fs/nfft) + fi;

    % Calculate a threshold to define the start & stop of this feature
    %
    %   Threshold is one-quarter of the way above the noise floor (in dB)
    %
    %                 noise floor [dB]    1/4 of the peak above the noise floor
    %                 _______|_______     _________________|__________________
    %                /               \   /                                    \
    %    thresh_dB = 10*log10(med_val) + (10*log10(ampl) - 10*log10(med_val))/4;
    %    thresh = 10^(thresh_dB/10);
    %
    % Simplified calculation in regular units
    %
    thresh = (median_value^(3/4) * ampl^(1/4));

    % Also set a 3dB-down threshold for interpolating the center frequency
    thresh_3dB = (ampl / 2);

    % Make sure the 3dB-down threshold is at least 'thresh' in magnitude
    if (thresh_3dB < thresh)
        thresh_3dB = thresh;
    end

    % Find the start index & stop index of this feature
    %
    %   This is the first index away from the feature (in the left & right
    %   directions) where the signal amplitude falls from its peak value
    %   of 'ht' to the first instance where it's less than the thresholds
    %   calculated above.
    %
    temp = find(shift_T(floor(idx+delta_idx):-1:1) <= thresh, 1, 'first');
    if isempty(temp)
        start_idx = 1;
    else
        start_idx = floor(idx + delta_idx) - temp + 1;
    end
    temp = find(shift_T(ceil(idx+delta_idx):nfft) <= thresh, 1, 'first');
    if isempty(temp)
        stop_idx = nfft;
    else
        stop_idx = ceil(idx + delta_idx) + temp - 1;
    end
    temp = find(shift_T(floor(idx+delta_idx):-1:1) <= thresh_3dB, 1, 'first');
    if isempty(temp)
        start_idx_3dB = 1;
    else
        start_idx_3dB = floor(idx + delta_idx) - temp + 1;
    end
    temp = find(shift_T(ceil(idx+delta_idx):nfft) <= thresh_3dB, 1, 'first');
    if isempty(temp)
        stop_idx_3dB = nfft;
    else
        stop_idx_3dB = ceil(idx + delta_idx) + temp - 1;
    end

    % Calculate the "bandwidth" estimates
    %
    %   Convert index range from bins to frequency
    %
    %   The bandwidth is calculated explicitly rather than converting
    %   the "width" return value from FINDPEAKS(), since the MATLAB
    %   function does not deal properly with peak asymmetry.
    %
    %   Note, the simple start & stop index search above finds the
    %   first indexes that are LESS THAN the threshold, so strictly
    %   speaking they over-estimate the extent of the spectral
    %   feature.  In order to be a bit more precise, without running
    %   a whole interpolation routine here, we'll just exploit the
    %   fact that the errors in the start & stop frequency estimates
    %   are distributed (0,1), so that the sum of their expected
    %   values is 1.0 (in units of bins).  That's where the '-1'
    %   comes from in the formula below.
    %
    %   Note, it is guaranteed by construction that the stop index
    %   will be at least 1 greater than the start index.
    %
    bw     = (  stop_idx   -   start_idx   - 1) * (fs/nfft);
    bw_3dB = (stop_idx_3dB - start_idx_3dB - 1) * (fs/nfft);

    % Calculate the start & stop frequencies of this feature
    %
    %   Convert from bins to frequency
    %
    %   Note, the bandwidth calculated above will not exactly match
    %   the difference of the start & stop frequencies, since we
    %   put in a '-1 bin' correction term into the bandwidth
    %   calculation in order to correct the over-estimate bias
    %   which would occur if we simply differenced the start & stop
    %   indexes that were found in the search above.
    %
    start_freq     = (-nfft/2 +   start_idx   - 1)*(fs/nfft) + fi;
    stop_freq      = (-nfft/2 +   stop_idx    - 1)*(fs/nfft) + fi;
    start_freq_3dB = (-nfft/2 + start_idx_3dB - 1)*(fs/nfft) + fi;
    stop_freq_3dB  = (-nfft/2 + stop_idx_3dB  - 1)*(fs/nfft) + fi;

    % Calculate the mid-point between the 3dB start & stop frequencise
    %
    %   We do this mid-point calculation in order that the frequency
    %   estimate will lie in the center of the feature, rather than
    %   potentially have the frequency estimate be biased by a
    %   highly-asymmetric feature shape.
    %
    freq_3dB = (start_freq_3dB + stop_freq_3dB) / 2;

    % Estimate the power
    %
    %   Note, to estimate the power of a CW (single-bin) feature:
    %
    %                spectrum   number of bins
    %                 height      in the FFT
    %                    |           |
    %                    v           v
    %       power = (amplitude^2 / nfft);
    %
    %   The following calculation estimates the power of a multi-bin
    %   feature, after converting from dB to regular units.
    %
    pwr = sum(shift_T(start_idx:stop_idx).^2) / nfft;

    % Correct for the noise bandwidth of the time-domain window
    pwr = pwr * (1/FFT_WIN_FUN_NOISE_BW);

    % Estimate the J/N power ratio, in dB
    j2n = 10*log10(pwr / noise_pwr);

    % If the J/N is above a threshold, then mark "detected"
    if (j2n >= PEAK_DETECTION_THRESHOLD)

        % Increment the peak counter
        %
        %   We increment the counter here since it's also used
        %   as the index into the 'spectral_peaks' structure.
        %
        peak_ctr = peak_ctr + 1;

        % Mark this peak as "detected"
        spectral_peaks(peak_ctr).is_detected = true;

        % Load feature parameters into return struct
        spectral_peaks(peak_ctr).idx        = idx;
        spectral_peaks(peak_ctr).idx_exact  = idx + delta_idx;
        spectral_peaks(peak_ctr).freq       = freq;
        spectral_peaks(peak_ctr).ampl       = ampl;
        spectral_peaks(peak_ctr).pwr        = pwr;
        spectral_peaks(peak_ctr).j2n        = j2n;
        spectral_peaks(peak_ctr).start_freq = start_freq;
        spectral_peaks(peak_ctr).stop_freq  = stop_freq;
        spectral_peaks(peak_ctr).bw         = bw;
        spectral_peaks(peak_ctr).bw_3dB     = bw_3dB;
        spectral_peaks(peak_ctr).freq_3dB   = freq_3dB;
        spectral_peaks(peak_ctr).phi        = phi;

    end

    % And "remove" this peak from the spectrum
    %
    %   While we only need to do this for the first (MAXPKS-1) peaks,
    %   we're going to do it for all for consistency purposes, in
    %   case we end up needing to continue processing the "de-peaked"
    %   spectrum for some reason.
    %
    %   Note, we explicitly do not check whether the peak was marked
    %   as "detected" [i.e., the J/N ratio was above our threshold],
    %   since we want to remove this peak from the spectrum even if
    %   it were "tall & skinny", thus falling below the J/N threshold.
    %   This way, we'll still have a chance of picking up possibly
    %   "short & wide" features which nonetheless have J/N greater
    %   than the threshold.
    %
    X = [(1:start_idx) (stop_idx:nfft)];
    shift_T(start_idx+1:stop_idx-1) ...
        = interp1(X, shift_T(X), start_idx+1:stop_idx-1);

    % Finally, check whether we've reached our MAXPKS target
    if (peak_ctr >= MAXPKS)
        % Terminate the loop
        break
    end

end % loop over peak number

% Quick bookkeeping check
%
%   Make sure we didn't introduce any new fields
%   into the SPECTRAL_PEAKS return struct.
%
if ~isempty(setxor(fieldnames(spectral_peaks),fieldnames(SPECTRAL_PEAKS)))
    error('YOU HAVE A PROBLEM!!');
end

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