function quantization
%
% quantization_levels.eps: Example of noise added from quantization
% quantization_noise.eps: "Equivalent" noise effect
%

%
% Copyright (c) 2008 by Theodore P. Pavlic
%
% This work is licensed under the Creative Commons
% Attribution-Noncommercial 3.0 United States License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc/3.0/us/
% or send a letter to Creative Commons, 171 Second Street, Suite 300,
% San Francisco, California, 94105, USA.
%

% Upper-case    A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
% Lower-case    a b c d e f g h i j k l m n o p q r s t u v w x y z
% Digits        0 1 2 3 4 5 6 7 8 9
% Exclamation   !           Double quote "          Hash (number) #
% Dollar        $           Percent      %          Ampersand     &
% Acute accent  '           Left paren   (          Right paren   )
% Asterisk      *           Plus         +          Comma         ,
% Minus         -           Point        .          Solidus       /
% Colon         :           Semicolon    ;          Less than     <
% Equals        =           Greater than >          Question mark ?
% At            @           Left bracket [          Backslash     \
% Right bracket ]           Circumflex   ^          Underscore    _
% Grave accent  `           Left brace   {          Vertical bar  |
% Right brace   }           Tilde        ~


%%%%%%%%%%%%%%%%%%%%%%%%%%% Start of Script %%%%%%%%%%%%%%%%%%%%%%%%%%%

B = 3;
LSB=2/2^B;
qlevels = ( -1:LSB:(1-LSB) )';

t = linspace(0,1,1000);

%large_signal = sin(2*pi*t-pi/8);
large_signal = .5*sin(2*pi*t-pi/8) + .5*sin(2*pi*2*t-pi/8);

%small_signal = 0.25*LSB*sin(2*pi*25*t);
small_signal = 0.25*LSB*sin(2*pi*25*exp(-2*t));

% Combine large and small signal; small signal gets lost in quantization
y = large_signal + small_signal;

% Add some noise to compare to quantization effect 
% (ideally we add 1 LSB^2 of noise power, but that looks ugly, so we tune
%  for pretty plots)
%ynoise = y + LSB^(1.5)/sqrt( (1/3)-(1/2)^2 )*(rand(1,length(t))-0.5);
ynoise = y + LSB^(1.55)*randn(1,length(t));


%%% Plot quantization stuff
figure(1); close(1); figure(1);

% Plot quantization noise region
hpatch = patch( [0;0;1;1], [0;LSB;LSB;0] - LSB, 0.65*[1 1 1] );
%set( hpatch, 'FaceAlpha', 0.5 );
hold on;

% Plot quantization levels
hql = plot( [0 1], [ qlevels qlevels ], 'md--' );
set( hql, 'MarkerSize', 6, 'MarkerFaceColor', 'magenta' );

% Plot continuous and discrete-value signal
[ax, hy, hq] = plotyy( t, y, t, q(y-LSB/2) );
set( hq, 'LineWidth', 2 );

% Make sure we're done holding
hold( ax(1), 'off' );
hold( ax(2), 'off' );

% Adjust axes to make quantization levels clear
xlim( [0 1] );
ylim( ax(1), [-1.1 1.1] ); 
ylim( ax(2), [-1.1 1.1] );

% Label left axis
ylabel( ax(1), 'Continuous Values', ...
               'Interpreter', 'latex', 'FontSize', 14 );
set( ax(1), 'YTick', [ qlevels; 1 ] );           
           
% Label right axis
ylabel( ax(2), 'Quantization Levels', ...
               'Interpreter', 'latex', 'FontSize', 14 );
set( ax(2), 'YTick', qlevels );
set( ax(2), 'YTickLabel', dec2bin( 0:(length(qlevels)-1), B) )
           
% Label time axis
xlabel( 'Time(s)', 'Interpreter', 'latex', 'FontSize', 14 );

% Title the plot
title( [ '\quad ' ...
         num2str(B) '-bit Quantization Noise Covers Up Small Signals' ], ...
       'Interpreter', 'latex', 'FontSize', 14 );

% Add a legend
legfig = legend( [hy hq hql(1) hpatch], ...
                 'Continuous-Value Signal', ...
                 'Quantized Signal', ...
                 'Quantization Levels', ...
                 ['A ``Noise Region'''''] );
set( legfig, 'Interpreter', 'latex', 'FontSize', 14 );

%% Save the figure

set( gcf, 'PaperType', 'usletter', ...
          'PaperOrientation', 'portrait', ...
          'PaperPosition', [0.25 2.5 7.75 4.75]);
saveas( gcf, 'quantization_levels.eps', 'epsc2' ); 



%%% Plot noise "equivalent"
figure(2); close(2); figure(2);
% Plot quantization noise region
hpatch = patch( [0;0;1;1], [0;LSB;LSB;0] - LSB, 0.65*[1 1 1] );
%set( hpatch, 'FaceAlpha', 0.5 );
hold on;

% Plot quantization levels
hql = plot( [0 1], [ qlevels qlevels ], 'md--' );
set( hql, 'MarkerSize', 6, 'MarkerFaceColor', 'magenta' );

% Plot continuous and discrete-value signal
[ax, hy, hn] = plotyy( t, y, t, ynoise );
delete(hy);

% Make sure we're done holding
hold( ax(1), 'off' );
hold( ax(2), 'off' );

% Adjust axes to make quantization levels clear
xlim( [0 1] );
ylim( ax(1), [-1.1 1.1] ); 
ylim( ax(2), [-1.1 1.1] );

% Label left axis
ylabel( ax(1), 'Continuous Values', ...
               'Interpreter', 'latex', 'FontSize', 14 );
set( ax(1), 'YTick', [ qlevels; 1 ] );           
           
% Label right axis
ylabel( ax(2), 'Quantization Levels', ...
               'Interpreter', 'latex', 'FontSize', 14 );
set( ax(2), 'YTick', qlevels );
set( ax(2), 'YTickLabel', dec2bin( 0:(length(qlevels)-1), B) )
           
% Label time axis
xlabel( 'Time(s)', 'Interpreter', 'latex', 'FontSize', 14 );

% Title the plot
title( [ '\quad ' ...
         'A Perceived-Noise Model of ' num2str(B) '-bit Quantization' ], ...
       'Interpreter', 'latex', 'FontSize', 14 );

% Add a legend
legfig = legend( [hn hql(1) hpatch], ...
                 'Signal Plus Noise', ...
                 'Quantization Levels', ...
                 'A ``Noise Region''''' );
set( legfig, 'Interpreter', 'latex', 'FontSize', 14, 'Color', [1 1 1] );

%% Save the figure

set( gcf, 'PaperType', 'usletter', ...
          'PaperOrientation', 'portrait', ...
          'PaperPosition', [0.25 2.5 7.75 4.75]);
saveas( gcf, 'quantization_noise.eps', 'epsc2' ); 




%% Helper functions

    % Calculates B-bit double-ended quantization (-V to V) of x
    function xq = q(x, B, V)

        if nargin < 3
            V = 1;
            
            if nargin < 2
                B = 3;
            end
        end
        
        xq = max( min( quant( x, 1/2^(B-1) ), V - V/2^(B-1) ), -V );

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%