function [fft_stats, T_filt] ...
             = max8_spectrum_properties(s_in, T_filt, FFT_STATS,             ...
                                      FFT_WIN_FUN, FFT_WIN_FUN_NOISE_BW,     ...
                                      FFT_SMOOTHING_COEF, FFT_IDX_OF_RF_SAW, ...
                                      NOISE_POWER_FINAGLING_FACTOR)
%-------------------------------------------------------------------------------
% $Id: max8_spectrum_properties.m,v 1.3 2023/07/26 21:49:41 ddelore Exp $
%
%  David S. De Lorenzo
%  12 December 2019
%
%  Estimate some spectrum properties
%
%-------------------------------------------------------------------------------

% Initialize the output struct
fft_stats = FFT_STATS;

% Copy variables to local storage
nsamp = fft_stats.nsamp;
nfft  = fft_stats.nfft;

% Check whether we're processing the time-domain signal
if ~isempty(s_in)

    % Consistency check
    if ( (nsamp ~= length(s_in))   ...  % the number of input samples
      || (nfft  ~= length(T_filt)) ...  % the number of FFT bins
      || (size(s_in,1) ~= 1) )          % the input must be a row vector
        error('You have a problem with your FFT data!!');
    end

    % Multiply the input signal by the time-domain anti-aliasing
    % window function and zero-pad to the desired FFT length.
    %   's' is complex {1 x nfft}
    s = [(s_in .* FFT_WIN_FUN), complex(zeros(1,nfft-nsamp))];

    % Perform the FFT
    %   'S' is complex {1 x nfft}
    %
    %   Recall, the input already has been multiplied by the time-
    %   domain window function and zero-padded to length 'nfft',
    %   which is why you only see the FFT() call here with no
    %   additional (optional) setting parameters.
    %
    %   We divide by sqrt(NSAMP) in order to normalize by the
    %   number of input samples that are getting sent to the FFT.
    %
    S = fft(s) / sqrt(nsamp);

    % Calculate the FFT magnitude
    %   'T' is real {1 x nfft}
    T = abs(S);

    % Filtered spectrum magnitude
    %   'T_filt' is real {1 x nfft}
    if ~all(isfinite(T_filt))
        % Initialization
        T_filt(:) = T;
    else
        % Apply a low-pass filter
        T_filt(:) = (1-FFT_SMOOTHING_COEF) * T_filt ...
                  +   FFT_SMOOTHING_COEF   *   T;
    end

end % if ~isempty(s_in)

% Shift the zero-frequency component to center of spectrum
%   shift_T_filt is real {1 x nfft}
%
%   This is done so that the index vector defining the FFT
%   frequencies that are within the bandwidth of the front-end
%   band-pass filter is monotonic (i.e., the frequencies do
%   not wrap-around).
%
shift_T_filt = fftshift(T_filt);

% Estimate the median value
%   recall, 'shift_T_filt' already is real-valued & non-negative
%
%   Here is where we limit the frequency range over which the
%   median value is determined to be within the bandwidth of
%   the front-end RF SAW band-pass filter.
%
% median_value = max8_quick_median(shift_T_filt(FFT_IDX_OF_RF_SAW)); % <-- SLOW
% median_value = max8_3pt_median(shift_T_filt(FFT_IDX_OF_RF_SAW));
median_value = median(shift_T_filt(FFT_IDX_OF_RF_SAW));

% Apply a frequency-domain noise power "finagling factor"
%
%   When estimating the median of a filtered spectrum, your
%   result comes out too large by a factor of ~sqrt(17/15).
%
median_value = median_value * sqrt(NOISE_POWER_FINAGLING_FACTOR);

% Calculate the standard deviation
std_dev = std(shift_T_filt(FFT_IDX_OF_RF_SAW));

% Calculate the mean absolute deviation about the median
%
%   This is the deviation of the filtered spectrum about the
%   median value --> where the median has been estimated from a
%   filtered spectrum in order to reduce noise on that estimate.
%
%   We also limit the frequency range over which the mean absolute
%   deviation is calculated to be within the filter bandwidth.
%
%   https://en.wikipedia.org/wiki/Average_absolute_deviation
%         #Mean_absolute_deviation_around_the_median
%   https://www.leeds.ac.uk/educol/documents/00003759.htm Note #4
%
mean_abs_deviation = mean(abs(shift_T_filt(FFT_IDX_OF_RF_SAW) - median_value));

% Bessel's modification, approximates the standard deviation
%
%   Note that the standard deviation estimate is for that
%   frequency range which is within the bandwidth of the
%   front-end RF SAW band-pass filter.
%
%   See S. Gorard (2005), Revisiting a 90-year-old debate: the
%   advantages of the mean deviation - particularly the comment
%   under Note #4
%   https://www.leeds.ac.uk/educol/documents/00003759.htm
%   http://bpi.enschool.org/ourpages/auto/2013/3/15/34398399/
%         Standard%20Deviation%20vs_%20Average%20Deviation.docx
%
%   It should be noted, since it's not covered elsewhere, that
%   this method of approximating the standard deviation breaks
%   down when the dynamic range of the sample distribution is too
%   large, meaning in those cases where the distribution is highly
%   non-Gaussian (i.e., when the ratio of the largest value(s) to
%   the mean or median is/are much greater than the square-root of
%   the second moment).
%
sigma_est = sqrt(pi/2) * mean_abs_deviation;

% Quick estimate of the noise power
%
%   We have already scaled the spectrum by sqrt(NSAMP), so
%   the noise power is effectively the square of the median
%   value, times a sqrt(2) for independent real & imaginary
%   samples.
%
noise_pwr = sqrt(2) * median_value^2;

% Correct for the noise bandwidth of the time-domain window
%
%   The spectrum median value actually depends on the choice
%   of time-domain anti-aliasing window.  Therefore, when
%   estimating the noise power of the signal based on the
%   square of the median, you have to correct by dividing
%   by the noise bandwidth of the time-domain window.
%
noise_pwr = noise_pwr * (1/FFT_WIN_FUN_NOISE_BW);

% Scale the noise power by the fraction of the spectrum which
% is within the bandwidth of the front-end band-pass filter.
%
noise_pwr = noise_pwr * (sum(FFT_IDX_OF_RF_SAW)/nfft);

% Quick estimate of the total in-band signal power
%
%   This value can be compared to the signal power
%   determined according to the time-domain variance
%   (assumes zero-mean signal -- i.e., no DC bias):
%
%     rssi == time-domain signal variance
%           = 1/length(s_in) * sum(s_in .* conj(s_in));
%           = mean((abs(s_in)).^2);
%
rssi = sum(shift_T_filt.^2) / nfft;

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

% At low J/N, this method may under-estimate the
% total in-band signal power.  So, make sure that
% RSSI is at least equal to the noise power.
if rssi < noise_pwr
    rssi = noise_pwr;
end

% Load signal parameters into return struct
fft_stats.median_value = median_value;
fft_stats.sigma_est    = sigma_est;
fft_stats.std_dev      = std_dev;
fft_stats.noise_pwr    = noise_pwr;
fft_stats.rssi         = rssi;

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