function sampling_intro
%
% sampling_intro_in.eps: Example of sampled input to DSP (+quantization)
% sampling_intro_out.eps: Example of DSP output (+quantization)
%

%
% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%% 

% "Continuous-" and Discrete-Time Vectors
%t = linspace(0,20,10000);
%ts = linspace(0,20,100);
t = 0:2/1000:20-2/1000;
ts = 0:2/10:20-2/10;

% "Continuous-" and Discrete-Time Data
y = exp(-0.1*t) .* cos(2*pi*0.1*t);
ys = exp(-0.1*ts) .* cos(2*pi*0.1*ts);

%% Quantize the data
B = 3;
qlevels = (1/2^(B-1))*(-(2^(B-1)):(2^(B-1)-1));
yq = q( ys, B );

%%% Now filter the data
% z = tf( 'z', ts(2)-ts(1) );
% sys = 1.2/(z^2+0.2);


yfs = zeros(1,length(ts));
yfq = yfs;

% Must quantize the quantized signal; DSP's are quantized everywhere
yfq(1) = q(yfq(1),B);
yfq(2) = q(yfq(2),B);
for i = 3:length(ts)
    
    % Standard DT filter
    yfs(i) = 1.2*ys(i-2) - 0.2*yfs(i-2);
    
    % Quantize result of every operation
    yfq(i) = q( q(q(1.2,B)*yq(i-2),B) - q(q(0.2,B)*yfq(i-2),B), B ); 

end

%%% And sinc interpolate the ideal filtered signal (will ZOH the quantized)

yfssinc = sinc_interp_slow( yfs, ts, t );

%%%% All the crazy plots

%% Input plots (pre-DSP)

figure(1); clf; 

dtfigs = stem(ts, ys, 'filled'); set(dtfigs, 'ShowBaseLine', 'off');
hold on;
qlfigs = ...
    plot( repmat( xlim', 1, 2^(B)), [ qlevels; qlevels ], 'md--' ); 
    set( qlfigs, 'MarkerSize', 6, 'MarkerFaceColor', 'magenta');
qfigs = stem(ts, yq, '*'); 
    set(qfigs, 'MarkerSize', 8, 'MarkerFaceColor', 'blue');
ctfigs = plot(t, y);
hold off;

legfigs = legend( [ctfigs, dtfigs, qfigs, qlfigs(1)], ...
                  'Continuous-Time', ...
                  'Discrete-Time Samples', ...
                  [ num2str(B) '-bit Quantized Samples' ], ...
                  [ num2str(B) '-bit Quantization Levels' ], ...
                  'Location', 'SouthEast' );
set( legfigs, 'Interpreter', 'latex', 'FontSize', 14);

xlabel('Time (s)', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('Data', 'Interpreter', 'latex', 'FontSize', 14);
title('\quad Quantized (``Digital'''') Samples of Continuous-Time Signal', ...
      'Interpreter', 'latex', 'FontSize', 14);
ylim1 = ylim;


%% "Output" plots (post-DSP)

figure(2); clf; 

ctffigs = plot( t, yfssinc );
hold on;
qlfigs = ...
    plot( repmat( xlim', 1, 2^(B)), [ qlevels; qlevels ], 'md--' ); 
    set( qlfigs, 'MarkerSize', 6, 'MarkerFaceColor', 'magenta');
dtfigs = plot(ts, yfs, 'o');
    set(dtfigs, 'MarkerSize', 6, 'MarkerFaceColor', 'blue');
qrfigs = stairs(ts, yfq);
    set( qrfigs, 'LineWidth', 1.5 );
qfigs = plot(ts, yfq, '*'); 
    set(qfigs, 'MarkerSize', 8, 'MarkerFaceColor', 'blue');
hold off;

legfigs = legend( [dtfigs, ctffigs, qfigs, qrfigs, qlfigs(1)], ...
                  'Filtered from Discrete-Time', ...
                  'Perfectly Reconstructed', ...
                  [ 'Filtered from ' num2str(B) ...
                    '-bit Quantization' ], ...
                  [ 'Held from ' num2str(B) ...
                    '-bit Quantization' ], ...
                  [ num2str(B) '-bit Quantization Levels' ], ...
                  'Location', 'SouthEast' );
set( legfigs, 'Interpreter', 'latex', 'FontSize', 14);

xlabel('Time (s)', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('Data', 'Interpreter', 'latex', 'FontSize', 14);
title('\quad Filtered and Reconstructed Signal', ...
      'Interpreter', 'latex', 'FontSize', 14);
ylim2 = ylim;
  
ylimmax = [ min( [ ylim1 ylim2 ] ) max( [ ylim1 ylim2 ] ) ];


%%% Save the figures

figure(1); ylim( ylimmax ); 
set( gcf, 'PaperType', 'usletter', ...
          'PaperOrientation', 'portrait', ...
          'PaperPosition', [0.25 2.5 7.75 4.75]);
saveas( gcf, 'sampling_intro_in.eps', 'epsc2' );
%close(1);

figure(2); 
ylim( ylimmax ); 
set( gcf, 'PaperType', 'usletter', ...
          'PaperOrientation', 'portrait', ...
          'PaperPosition', [0.25 2.5 7.75 4.75]);
saveas( gcf, 'sampling_intro_out.eps', 'epsc2' ); 
%close(2);


%%% Helper functions  

    % Ideally "resamples" x vector from s to u by sinc interpolation
    %
    % NOTE: This *SLOW* version is provided for explanatory power.
    %       Actually use the fast sinc_interp below.
    %
    function y = sinc_interp_slow(x,s,u) 
        % Interpolates x sampled sampled at "s" instants
        % Output y is sampled at "u" instants ("u" for "upsampled")
       
        % Find the period of the undersampled signal
        T = s(2)-s(1);
      
        for i=1:length(u) 
            y( i ) = sum( x .* sinc( (1/T)*(u(i) - s) ) );
        end
    end

    % 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

    % Ideally "resamples" x vector from s to u by sinc interpolation
    %
    % NOTE: This is a fast version of sinc_interp_slow above.
    %       Its speed comes from all operations being "vectorized."
    %       That is, we avoid the use of "for" loops.
    %
    %       If you need an even faster version, just use the
    %       
    %          interpft
    %
    %       command. It does an FFT of the data, then it pads it with
    %       zeros, and then it does an IFFT to get the data back. It
    %       basically does the following (assuming length(x) is EVEN):
    %
    %         length( t )/length( ts ) * ...
    %         ifft( fftshift( ...
    %           [ zeros(1,length(u)/2) ...    % Negative padding
    %             fftshift( fft(x) ) ...
    %             zeros(1,length(u)/2) ] ...  % Positive padding
    %              ) )
    %
    %       The FFT algorithm is VERY FAST in general *AND* it's
    %       implemented as a MATLAB internal, so it's even faster.
    %       However, make sure that length(u)/length(s) is an
    %       INTEGER. Also make sure that u lands on the same points as
    %       s. That is, u should be a version of s with points ADDED.
    %       (e.g., s = 0:0.1:10-0.1; u = 0:0.01:10-0.01;).
    %
    function y = sinc_interp(x,s,u)
        % Interpolates x sampled sampled at "s" instants
        % Output y is sampled at "u" instants ("u" for "upsampled")
        
        % Find the period of the undersampled signal
        T = s(2)-s(1);

        % When generating this matrix, remember that "s" and "u" are
        % passed as ROW vectors and "y" is expected to also be a ROW
        % vector. If everything were column vectors, we'd do.
        %
        % sincM = repmat( u, 1, length(s) ) - repmat( s', length(u), 1 );
        %
        % So that the matrix would be longer than it is wide. 
        % Here, we generate the transpose of that matrix.
        sincM = repmat( u, length(s), 1 ) - repmat( s', 1, length(u) );
        
        % Equivalent to column vector math:
        %   y = sinc( sincM'/T )*x';
        y = x*sinc( sincM/T );
        
    end

end

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