%==============================================================================%
% $Id: adc_quantization_study.m,v 1.8 2023/07/26 21:49:41 ddelore Exp $
%                                                                              %
%  Determine signal scaling for various quantizers                             %
%  David S. De Lorenzo                                                         %
%  February 05, 2021                                                           %
%                                                                              %
%------------------------------------------------------------------------------%
%                                                                              %
%  A signal of arbitrary power comes in.  How ought that signal be scaled      %
%  such that the post-A/D distribution follows that of the desired quantizer.  %
%  Specifically, when we reserve "M0" bits for WGN quantization yet have a     %
%  total of "M" bits of total dynamic range, what is the post-A/D signal       %
%  variance, in units of the A/D output, that the WGN should take.             %
%                                                                              %
%% -- NOTES ---------------------------------------------------------------%
%
%  This study will determine:
%
%    1. The desired post-A/D noise variance
%
%       This is in digitized "post-quantization" units, and
%       represents the desired signal variance of the WGN-like
%       component in the sampled data, as estimated (for example)
%       by the noise floor of an averaged spectrum.
%
%       Note:  the noise variance has a dependence not only on
%       the value of M0 (the number of A/D bits allocated to WGN)
%       but also on M (the total dynamic range).  Thus, when running
%       this study to generate the desired post-A/D noise varniance,
%       you'll want to set the simulation parameters as follows:
%
%         (1) set "M" to a high value --> I used M = 12
%         (2) set "npts" to a high value --> I used npts = 1e7
%         (3) use "complex" data type, otherwise variance is 1/2
%         (4) turn off plotting, since it's uninformative when M >> M0
%
%    2. Post-A/D to pre-A/D noise variance scale factor
%
%       This relates the signal variance of the digitized WGN-like
%       component in the sampled data to the signal variance of the
%       analog pre-A/D signal.
%
%-------------------------------------------------------------------------------
%
%  References:
%
%  1. F. Bastide, et al., "Automatic Gain Control (AGC) as an Interference
%     Assessment Tool", Proc. ION GPS/GNSS 2003, pgs. 2042-2053.
%
%     pg. 2044:
%       "[t]o control the AGC gain ... map the ADC output signal power
%       [DSD:  meaning in post-ADC units] to the ADC input signal power.
%       Then using a look-up table, the receiver may estimate the input
%       noise standard deviation by estimating quantized signal power.
%     
%==============================================================================%

%% --- CLEAR VARIABLES & SET SIMULATION PARAMETERS ----------------------

home

clear;
fclose('all');
path(pathdef);

QUANTIZER_TYPE = 'centered';
% QUANTIZER_TYPE = 'non-centered';

% The code works for both real & complex data, although the scale
% factors are stored for complex data type only.
% signal_type = 'real';
signal_type = 'complex';

% Set input variance
%   This will get scaled-out in the simulation below, it's
%   included because it was a useful parameter for debugging.
input_variance = 13.5;

% The A/D converter peak-to-peak voltage
%   This will get scaled-out in the simulation below, it's
%   included because it was a useful parameter for debugging.
ADC_p2p = 1.800;

M0      =   5;      % the number of A/D levels to be excited by WGN
M       =  10;      % the number of full-scale A/D quantization levels

% The number of points to generate in the call to RANDN()
% npts = 1e7;     % slow, but more accurate
npts = 1e5;     % fast, but less accurate

%% --- A/D CONVERTER SETUP -------------------------------------------------

% The A/D converter specification function
[ADC,k,~,~,sf] = max8_adc_fn(QUANTIZER_TYPE, M0, M);

% Transcribe short-hand variables names to longer
% and more descriptive names.
ADC_OPTIMAL_RATIO               = k;
ADC_NOISE_VARIANCE_SCALE_FACTOR = sf;

%% --- PROCESSING ----------------------------------------------------------

switch signal_type
    
    case 'real'
        
        % Create signal vector
        s = sqrt(input_variance) * real(randn(1,npts));
        
        % Scale signal for quantization
        %
        %   The input variance scaling is by sqrt(input_variance)
        %
        switch QUANTIZER_TYPE
            case 'centered'
                t = ( s / sqrt(input_variance) ) ...
                  * ( (ADC_p2p/2) / ADC_OPTIMAL_RATIO ) ...
                  * ( (2^M0 - 3) / (2^M - 3) );
            case 'non-centered'
                t = ( s / sqrt(input_variance) ) ...
                  * ( (ADC_p2p/2) / ADC_OPTIMAL_RATIO ) ...
                  * ( (2^M0 - 2) / (2^M - 2) );
            otherwise
                error('Invalid quantization option!!');
        end
        
        % Quantization
        u = ADC(real(t),M,ADC_p2p/2);
        
    case 'complex'
        
        % Create signal vector
        s = sqrt(input_variance/2) * complex(randn(1,npts), ...
                                             randn(1,npts));
        
        % Scale signal for quantization
        %
        %   The input variance scaling is by sqrt(input_variance/2)
        %
        switch QUANTIZER_TYPE
            case 'centered'
                t = ( s / sqrt(input_variance/2) ) ...
                  * ( (ADC_p2p/2) / ADC_OPTIMAL_RATIO ) ...
                  * ( (2^M0 - 3) / (2^M - 3) );
            case 'non-centered'
                t = ( s / sqrt(input_variance/2) ) ...
                  * ( (ADC_p2p/2) / ADC_OPTIMAL_RATIO ) ...
                  * ( (2^M0 - 2) / (2^M - 2) );
            otherwise
                error('Invalid quantization option!!');
        end
        
        % Quantization
        u = complex(ADC(real(t),M,ADC_p2p/2), ...
                    ADC(imag(t),M,ADC_p2p/2));
        
    otherwise
        
        error('Invalid signal type option!!');
        
end

% Pre/Post-A/D variance
variance_pre_scaling = 1/length(s) * sum(s .* conj(s));
variance_pre_ADC     = 1/length(t) * sum(t .* conj(t));
variance_post_ADC    = 1/length(u) * sum(u .* conj(u));

post_to_pre_scale_factor_measured ...
    = ( variance_post_ADC/variance_pre_ADC * (ADC_p2p/2)^2 );

fprintf('\nQUANTIZER TYPE = "%s"\n',   QUANTIZER_TYPE);
fprintf(  'SIGNAL TYPE    =   "%s"\n\n', signal_type);
fprintf('    No. A/D levels for WGN = %2d\n', M0);
fprintf('    Full-scale A/D levels  = %2d\n\n', M);
fprintf('    Pre-scaling variance      = %6.3f\n',   variance_pre_scaling);
fprintf('    Pre-ADC variance          = %6.5f\n',   variance_pre_ADC);
fprintf('    Desired Post-ADC variance = %6.3f  <-----\n', variance_post_ADC);
fprintf('           10*log10(variance) = %6.3f dB\n\n', ...
       10*log10(variance_post_ADC));
fprintf('  Post-to-Pre Scale Factor --> predicted = %.2f\n', ...
        ADC_NOISE_VARIANCE_SCALE_FACTOR);
fprintf('                           --> measured  = %.2f\n\n', ...
        post_to_pre_scale_factor_measured);

%% --- PLOTTING ------------------------------------------------------------

%{
% Plotting can be time-consuming, so feel free to
% bypass this code block if you're looking for speedup.

figure(1);
clf(1);

box('on');
grid('on');
hold('on');

switch signal_type
    case 'real'
        plot(real(t), real(u), ...
             'bo', 'LineWidth', 1);
    case 'complex'
        plot([real(t), imag(t)], [real(u), imag(u)], ...
             'bo', 'LineWidth', 1);
    otherwise
        error('Invalid signal type option!!');
end

switch QUANTIZER_TYPE
    case 'centered'
        xlim(1.25*ADC_p2p/2*[-1 1]);
        ylim(1.2*(2^(M-1))*[-1 1]);
        line(-ADC_p2p/2*[1 1], ylim, ...
             'Color', 'r', 'LineWidth', 2);
        line(+ADC_p2p/2*[1 1], ylim, ...
             'Color', 'r', 'LineWidth', 2);
        line(xlim, (-2^(M-1)+1)*[1 1], ...
             'Color', 'r', 'LineWidth', 2);
        line(xlim, (+2^(M-1)-1)*[1 1], ...
             'Color', 'r', 'LineWidth', 2);
        line(xlim, (-2^(M0-1)+1)*[1 1], ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
        line(xlim, (+2^(M0-1)-1)*[1 1], ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
    case 'non-centered'
        xlim(1.25*ADC_p2p/2*[-1 1]);
        ylim(1.2*(2^M)*[-1 1]);
        line(-ADC_p2p/2*[1 1], ylim, ...
             'Color', 'r', 'LineWidth', 2);
        line(+ADC_p2p/2*[1 1], ylim, ...
             'Color', 'r', 'LineWidth', 2);
        line(xlim, (-2^M+1)*[1 1], ...
             'Color', 'r', 'LineWidth', 2);
        line(xlim, (+2^M-1)*[1 1], ...
             'Color', 'r', 'LineWidth', 2);
        line(xlim, (-2^M0+1)*[1 1], ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
        line(xlim, (+2^M0-1)*[1 1], ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
    otherwise
        error('Invalid quantization option!!');
end

xlabel('Input Sample Value [V]', 'FontSize', 12);
ylabel('Quantizer Output', 'FontSize', 12);
title('Quantizer Function', 'FontSize', 14);
drawnow();

figure(2);
clf(2);

box('on');
grid('on');
hold('on');

switch signal_type
    case 'real'
        h = histogram(real(u), ...
                      'BinLimits', [-2^M 2^M], ...
                      'BinMethod', 'integers', ...
                      'Normalization', 'probability');
    case 'complex'
        h = histogram([real(u),imag(u)], ...
                      'BinLimits', [-2^M 2^M], ...
                      'BinMethod', 'integers', ...
                      'Normalization', 'probability');
    otherwise
        error('Invalid signal type option!!');
end

switch QUANTIZER_TYPE
    case 'centered'
        xlim(1.2*(2^(M-1))*[-1 1]);
        ylim(h.Parent.YLim);
        line((-2^(M-1)+1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 2);
        line((+2^(M-1)-1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 2);
        line((-2^(M0-1)+1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
        line((+2^(M0-1)-1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
        if (M0 == 2)
            plot(0, 0.462, 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        elseif (M0 == 3)
            plot(0, 0.254, 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        elseif (M0 == 4)
            plot(0, 0.139, 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        elseif (M0 == 5)
            plot(0, 0.077, 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        end
    case 'non-centered'
        xlim(1.2*(2^M)*[-1 1]);
        ylim(h.Parent.YLim);
        line((-2^M+1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 2);
        line((+2^M-1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 2);
        line((-2^M0+1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
        line((+2^M0-1)*[1 1], h.Parent.YLim, ...
             'Color', 'r', 'LineWidth', 1, 'LineStyle', '--');
        if (M0 == 2)
            plot([-1 1], 0.342*[1 1], 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        elseif (M0 == 3)
            plot([-1 1], 0.222*[1 1], 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        elseif (M0 == 4)
            plot([-1 1], 0.132*[1 1], 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        elseif (M0 == 5)
            plot([-1 1], 0.074*[1 1], 'Color', 'k', ...
                 'LineStyle', 'none', 'Marker', 'o', ...
                 'LineWidth', 2, 'MarkerSize', 8);
        end
    otherwise
        error('Invalid quantization option!!');
end

xlabel('Quantized Value', 'FontSize', 12);
ylabel('Probability', 'FontSize', 12);
title('Quantized Signal Histogram', 'FontSize', 14);
drawnow();
%}

%% --- SPECTRUM ESTIMATES --------------------------------------------------

% Finagling factor going from median value to noise power
%
%   When estimating noise power based on the median value of
%   an averaged spectrum, your power calculation comes out too
%   large by a factor of ~(17/15) = 1.133... .
%
%   So, the finagling factor multiplies the number you get by (15/17)
%   in order to produce the number that you should have gotten.
%
NOISE_POWER_FINAGLING_FACTOR = (15/17);

fs    =  50e6;      % sampling frequency, Hz
fi    = 1590e6;     % intermediate frequency, Hz

nread = 128;

nsamp = nread;
nfft  = 2*nread;

% Hamming window function
%   (highest sidelobe level = -43 dB)
w = 0.54 - 0.46*cos(2*pi*(0:nsamp-1)/nsamp);
fft_win_fun_noise_bw = 1.36;

% Re-scale the window to preserve time-domain amplitude
w = (w / mean(w));

% FFT of the post-A/D signal 'u'
idx = 1:(floor(npts/nsamp)*nsamp);
S = fft([reshape(u(idx),nsamp,[]) ...
                 .*repmat(reshape(w,[],1),1,length(idx)/nsamp); ...
         complex(zeros((nfft-nsamp),length(idx)/nsamp))],[],1);

% Normalize by the number of input samples
S = S / sqrt(nsamp);

% Calculate & average the FFT magnitude
T = mean(abs(S),2);

% Median absolute value
median_value = median(T);

% Quick estimate of the noise power
noise_pwr = sqrt(2) * median_value^2 / fft_win_fun_noise_bw;

% Apply a frequency-domain noise power "finagling factor"
noise_pwr = noise_pwr * NOISE_POWER_FINAGLING_FACTOR;

fprintf('    Variance by spectrum est. = %6.4f\n\n', noise_pwr);

%% -----------------------------------------------------------------------------
return
%===============================================================================