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