%==============================================================================%
% $Id: adc_saturation_study.m,v 1.6 2023/07/26 21:49:41 ddelore Exp $
%                                                                              %
%  Plot Sinusiodal RFI & WGN Histograms                                        %
%  David S. De Lorenzo                                                         %
%  February 05, 2021                                                           %
%                                                                              %
%==============================================================================%
%#ok<*NOPTS>

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

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

M  = 10;            % total # of A/D quantization bits
M0 =  4;            % # of A/D bits allocated for WGN
ADC_p2p = 1.8;      % A/D converter peak-to-peak range, [V]

N = 1e5;            % number of data points in test vectors

% Look up the optimal ratio 'k' of the converter's maximum
% quantization threshold 'L' to the input signal's noise
% standard deviation 'sigma'.
[~,k,~,~,~] = max8_adc_fn(QUANTIZER_TYPE, M0, M);

%% This is a sinusoidal RFI signal -----------------------------------------

% Signal power, in dBm
% P_desired_dBm = -3.02;
P_desired_dBm = -93.85;

% Convert to RMS voltage, in volts
%
%   V = sqrt(power x resistance)
%
%                     dBm          mW to
%                    to mW         Watts   50 ohms
%              ________|__________    |      |
%             /                   \   v      v
v_rms = sqrt((10^(P_desired_dBm/10)/1000) * 50);

% Convert to peak-to-peak voltage, in mV
%
%            RMS    max          convert
%             to     to           to RMS
%            max   p-to-p       millivolts
%             |      |              |
%             v      v              v
v_p2p_mV = sqrt(2) * 2 * (v_rms * 1000)

% Create a sinusoidal signal of correct peak-to-peak amplitude
%
%              peak-to-peak        just
%                back to            one
%               amplitude          cycle
%                   |             ___|___
%                   v            /       \
s_rfi_mV = v_p2p_mV/2 * sin(2*pi*(0:N-1)/N);

% Calculate the power of the RFI signal, in Watts
%
%   P = V^2 / R
%
%            calculate      calculate     this is             mV^2 to     50
%             average          sum       squared mV           volts^2    ohms
%         _______|_______       |    _________|___________       |        |
%        /               \      v   /                     \      v        v
P_rfi = (1/length(s_rfi_mV)) * sum((s_rfi_mV.*conj(s_rfi_mV))/(1000^2)) / 50;

% Convert power to dBm -- this should match your value 'P'
%
%            mW to   Watts   input
%             dBm    to mW   Watts
%            __|___    |      |
%           /      \   v      v
P_rfi_dBm = 10*log10(1000 * P_rfi)

%% This is a WGN signal ----------------------------------------------------
%
% The WGN signal has the optimal ratio of the analog-to-digital
% converter's maximum quantization threshold to the RFI signal's
% noise standard deviation.

% What is the 1-sigma voltage for the WGN
%
%   This is RFI 'RMS voltage' divided by the ADC ratio 'k'
%
%             ADC      p-to-p     amplitude
%            ratio     back to        to
%             'k'     amplitude      RMS
%              |      ____|___        |
%              v     /        \       v
v_wgn_std = (1/k) * (v_p2p_mV/2) / sqrt(2);

% Create a WGN signal of correct amplitude
%
%   This has the WGN +/- tails at the RFI peak-to-peak level
%
%         scale by        unit
%         standard      variance
%         deviation    WGN signal
%              |           |
%              v           v
s_wgn_mV = v_wgn_std * randn(1,N);

% Calculate the power of the WGN signal, in Watts
%
%   P = V^2 / R
%
%            calculate      calculate     this is             mV^2 to     50
%             average          sum       squared mV           volts^2    ohms
%         _______|_______       |    _________|___________       |        |
%        /               \      v   /                     \      v        v
P_wgn = (1/length(s_wgn_mV)) * sum((s_wgn_mV.*conj(s_wgn_mV))/(1000^2)) / 50;

% Convert power to dBm -- this can compare to your RFI power 'P'
%
%            mW to   Watts   input
%             dBm    to mW   Watts
%            __|___    |      |
%           /      \   v      v
P_wgn_dBm = 10*log10(1000 * P_wgn)

%-------------------------------------------------------------------------------
% Now, let's calculate the desired WGN power that actually "fills"
% the lower M0 bits of an M-bit analog-to-digital converter.

% The number of millivolts per quantization bit
%
%                     volts      # of
%                     to mV      bins
%                        |       _|_
%                        v      /   \
mv_per_bit = (ADC_p2p * 1000) / (2^M);

% The "target" quantizer threshold for M0 bits
%
%                          peak-to-peak
%                  # of      back to
%                  bins     one-sided
%                  _|__      /
%                 /    \    v
L = (mv_per_bit * (2^M0)) / 2;

% What is the 1-sigma voltage for the WGN
%
%   This is the target quantizer threshold divided by the ADC ratio 'k'
%
%  max.   ADC
%  ADC   ratio
% value   'k'
%   \     /
%    v   v
v = (L / k)

% Create a WGN signal of correct amplitude
s = v * randn(1,N);

% Calculate the power of the WGN signal, in Watts
P = (1/length(s)) * sum((s.*conj(s))/(1000^2)) / 50;

% Convert power to dBm
P_dBm = 10*log10(1000 * P)

%% Plot --------------------------------------------------------------------

figure(1);
clf(1);

subplot(2,1,1);
box('on');
grid('on');
hold('on');
histogram(s_rfi_mV,40,'DisplayName','RFI Signal at A/D Converter');
histogram(s_wgn_mV,40,'DisplayName','WGN of Comparable Range');
xlabel('Signal Level [mV]','FontSize',12)
title('Histograms for RFI & WGN','FontSize',14);
legend('show','Location','EastOutside');

subplot(2,1,2);
box('on');
grid('on');
hold('on');
histogram(s,40,'DisplayName','Desired WGN Quantization');
xlabel('Quantization Value','FontSize',12)
legend('show','Location','EastOutside');

%% The AGC Gain vs. Noise Figure relation from Dave Attard -----------------

%{
% AGC Gain [dB]
x = [-62.45 -54.83 -47.57 -40.66 -34.08 -27.83 -21.88 -16.23 -10.85 -5.75, ...
     -0.90 3.69 8.06 10.79 15.97 20.03 24.11 27.50 30.70 33.73 36.43 38.99, ...
     41.36 43.59 45.88 47.97 50.07 51.99 54.00 55.82 57.62 59.50 61.22 62.90, ...
     64.47 65.98 67.50 68.80 70.09 71.37 72.65 73.77 74.86 75.76 76.70 77.54, ...
     78.47 79.50 80.48 81.39 82.25 83.05 83.80 84.48 85.11 85.68 86.19 86.65, ...
     87.05 87.39 87.67 87.90 88.07 88.18 88.24 88.24 88.25 88.26 88.29];

% AGC N.F. [dB]
y = [63.00 62.27 61.01 59.15 56.87 54.31 51.55 48.71 45.83 42.98 40.20 37.52, ...
     34.95 33.14 30.94 27.83 25.84 23.61 21.62 20.07 18.81 17.76 16.65 15.74, ...
     14.94 13.99 13.42 12.97 11.62 11.10 10.30 9.84 9.20 8.76 8.59 8.08 7.92, ...
     7.58 7.43 7.05 6.99 6.96 6.92 6.91 7.09 7.12 7.13 7.15 7.16 7.17 7.19, ...
     7.20 7.22 7.23 7.24 7.26 7.27 7.28 7.30 7.31 7.32 7.32 7.32 7.32 7.32 7.32, ...
     7.32 7.32 7.32];

% Fit a 3rd-order polynomial to the data
p = polyfit(x,y,3);

xf = -62:88;
yf = p(1)*xf.^3 + p(2)*xf.^2 + p(3)*xf + p(4);
%}

%{
% LT5546
% AGC Gain [dB]
x = [-6.00 -4.00 -2.00 0.00 2.00 4.00 6.00 8.00 10.00 12.00, ...
     14.00 16.00 18.00 20.00 22.00 24.00 26.00 28.00 30.00, ...
     32.00 34.00 36.00 38.00 40.00 42.00 44.00 46.00 48.00, ...
     50.00 52.00 54.00 56.00 56.00 56.00 56.00 56.00 56.00];

% AGC N.F. [dB]
y = [48.00 45.64 43.35 41.13 38.98 36.89 34.86 32.90 31.01, ...
     29.19 27.43 25.74 24.11 22.55 21.06 19.63 18.27 16.98, ...
     15.75 14.59 13.49 12.47 11.50 10.61 9.78 9.01 8.32 7.69, ...
     7.12 6.63 6.19 5.83 5.53 5.30 5.13 5.03 5.00];

 % Fit a 3rd-order polynomial to the data
p = polyfit(x,y,3);

xf = -6:56;
yf = p(1)*xf.^3 + p(2)*xf.^2 + p(3)*xf + p(4);
%}

% Plot
% figure(2);
% clf(2);
% box('on');
% grid('on');
% hold('on');
% plot(x, y,  'b-',  'LineWidth', 3);
% plot(xf,yf, 'r--', 'LineWidth', 2);

%% Finally:  percent saturation vs. "excess gain" --------------------------
%
%   Nominal gain will saturate to the desired level 'k' for the
%   specified analog-to-digital converter; "excess gain" will
%   drive the ADC further into saturation.  We wish to characterize
%   that excess gain vs. the percent of quantized samples that
%   exceed the converter's maximum quantization threshold.

% The maximum quantizer threshold for M bits, mV
L = (ADC_p2p/2) * 1000;

% Create a sinusoidal signal, mV
%
%  amplitude      just one cycle
%      |             ___|___
%      v            /       \
s_cw = L * sin(2*pi*(0:N-1)/N);

% Look up the optimal ratio 'k' of the converter's maximum
% quantization threshold 'L' to the input signal's noise
% standard deviation 'sigma'.
%
%                        we use 'M' as the value for the number of
%                        noise bits since that will fully saturate
%                             the analog-to-digital converter
%                                         |
%                                         v
[~,k,~,~,~] = max8_adc_fn(QUANTIZER_TYPE, M, M);

% Create a WGN signal, mV
%
%      max.   ADC
%      ADC   ratio
%     value   'k'
%       \     /
%        v   v
s_wgn = (L / k) * randn(1,N);

gg  = linspace(-3, 17, 21);
pwr_cw  = NaN(size(gg));
pct_cw  = NaN(size(gg));
pwr_wgn = NaN(size(gg));
pct_wgn = NaN(size(gg));

for ii = 1:length(gg)
    
    %--- Sinusoidal ------------------------------------------------------------
    %
    % Scale the signal by the gain value
    s_temp = s_cw * sqrt(10^((gg(ii))/10));
    
    % Calculate the power of the CW signal, in Watts
    P = (1/length(s_temp)) * sum((s_temp.*conj(s_temp))/(1000^2)) / 50;
    
    % Convert power to dBm
    pwr_cw(ii) = 10*log10(1000 * P);

    % Record the percent outside the quantization limits
    pct_cw(ii) = sum(abs(s_temp) > L) / length(s_temp);
    
    %--- WGN -------------------------------------------------------------------
    %
    % Scale the signal by the gain value
    s_temp = s_wgn * sqrt(10^(gg(ii)/10));
    
    % Calculate the power of the WGN signal, in Watts
    P = (1/length(s_temp)) * sum((s_temp.*conj(s_temp))/(1000^2)) / 50;
    
    % Convert power to dBm
    pwr_wgn(ii) = 10*log10(1000 * P);
    
    % Record the percent outside the quantization limits
    pct_wgn(ii) = sum(abs(s_temp) > L) / length(s_temp);
    
end

figure(3);
clf(3);

subplot(2,1,1);
box('on');
grid('on');
hold('on');
plot(gg,pct_cw,'DisplayName','CW');
plot(gg,pct_wgn,'DisplayName','WGN');
xlabel('Gain [dB]','FontSize',12)
ylabel('Percent Outsize Quantization Limits','FontSize',12);
legend('show','Location','EastOutside');

subplot(2,1,2);
box('on');
grid('on');
hold('on');
plot(gg,pwr_cw,'DisplayName','CW');
plot(gg,pwr_wgn,'DisplayName','WGN');
xlabel('Gain [dB]','FontSize',12)
ylabel('Power [dBm]','FontSize',12);
legend('show','Location','EastOutside');

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