function sampling_ambiguity
%
% sampling_ambiguity.eps: Example of energy spreading via sampling
%

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


%%% Generate ambiguity plot

t = 0:1/1000:0.25;
ts = 0:1/10:0.25;

% Generate a signal
y = sin(2*pi*6*t);

% Sample it
ys = sin(2*pi*6*ts);

% Construct different aliases for it
y1 = sin(2*pi*-4*t);
y2 = sin(2*pi*-14*t);
y3 = sin(2*pi*16*t);

figure(1); clf;

pfigs = plot(t,y1,':',t,y,'-.',t,y2,'-',t,y3,'--'); 
    set( pfigs, 'LineWidth', 2 );
hold on;
sfigs = stem(ts,ys,'filled'); 
    set( sfigs, 'MarkerSize', 8 );
hold off;

lfig = legend([pfigs; sfigs], ... 
              '$-4$ Hz', '$6$ Hz', '$-14$ Hz', '$16$ Hz', ...
              'Samples @ 10 S/s', ...
              'Location', 'SouthEast');
    set(lfig, 'Interpreter', 'latex', 'FontSize', 14);
         
xlabel( 'Time (s)', 'Interpreter', 'latex', 'FontSize', 14 );
ylabel( 'Data', 'Interpreter', 'latex', 'FontSize', 14 );
title( '\quad Sampling Ambiguity', 'Interpreter', 'latex', 'FontSize', 14 );


%%% Save the first figure

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


%% Helper functions

    % Ideally "resamples" x vector from s to u by sinc interpolation
    % (this function will not be used in this script; it is provided here
    %  as an example)
    %
    % NOTE: This *SLOW* version is provided for explanatory power.
    %       Actually use the fast sinc_interp below.    
    %
    function y = sinc_interp_slow(x,s,u,N) 
        % Interpolates x sampled sampled at "s" instants
        % Output y is sampled at "u" instants ("u" for "upsampled")
        % Optionally, uses the Nth sampling window where N=0 is DC
        %     (so non-baseband signals have N = 1,2,3,...)
        
        if nargin < 4
            N = 0;
        end
        
        % Find the period of the undersampled signal
        T = s(2)-s(1);
      
        for i=1:length(u) 
            y( i ) = ...
                sum( x .* ...
                     ( (N+1)*sinc( ((N+1)/T)*(u(i) - s) ) - ...
                       N*sinc( (N/T)*(u(i) - s) ) ) );
        end
    end

    % Ideally "resamples" x vector from s to u by sinc interpolation
    % (this function will not be used in this script; it is provided here
    %  as an example)
    %
    % 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 N=0, it's MUCH faster to use the
    %
    %         interpft( x, length(u) )
    %
    %       command. 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;).
    %
    %       If N>0, you can use fft on x and pad the result with zeros
    %       until you form the desired fft of your desired signal. That
    %       is, you're shifting the FFT into the appropriate aliasing
    %       window. Once complete, do an ifft, and the result will have
    %       as many points as the number of zeros you've added. This
    %       process is identical to what INTERPFT does, but INTERPFT
    %       adds all of its zeros to the "ends" of the fft. When doing
    %       your zero padding, the "fftshift" function may be handy.
    %
    function y = sinc_interp(x,s,u,N) 
        % Interpolates x sampled sampled at "s" instants
        % Output y is sampled at "u" instants ("u" for "upsampled")
        % Optionally, uses the Nth sampling window where N=0 is DC
        %     (so non-baseband signals have N = 1,2,3,...)
        
        if nargin < 4
            N = 0;
        end
        
        % 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'(N+1)/T )*x';
        y = x*( (N+1)*sinc( sincM*(N+1)/T ) - N*sinc( sincM*N/T ) );

    end

end

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