% this example shows sampling theorem and interpolation
% time -domain function is sinc^2(2*pi*t), 
% Fourier domain function is triangle
%
% sinc^2(2*pi*t) has Nyquist frequency of 4, so delta_t is pi/2
%
%build a high resolution function as the standard to compare against
thi = [-10*pi:pi/100:10*pi-pi/100];
shi = sinc(thi/pi).^2; % Matlab defines sinc(x) as sin(pi*x)/(pi*x)!
figure(1);hold;plot(thi,shi,'g');
whi = [-2:0.001:1.999];
fhi = [2+whi(1:2000),2-whi(2001:4000)];
figure(2);hold;plot(whi,fhi,'g');
%
% First choose T = 10*pi
T = 10*pi;
delta_t = pi/2;
ts1 = [-T/2:delta_t:T/2-delta_t];
% sample
ss1 = sinc(ts1/pi).^2; % Matlab defines sinc(x) as sin(pi*x)/(pi*x)!
%display
figure(1);plot(ts1,ss1,'b+');
% Fourier domain
w0 = 2*pi/T;
W = 2*pi/delta_t;
fs1 = fft(ss1);
ws1 = [-W/2:w0:W/2-w0];
figure(2);plot(ws1,fftshift(abs(fs1)),'b+');
% 
% Now Keep delta_t fixed, choose T = 20*pi
T = 20*pi;
delta_t = pi/2;
ts2 = [-T/2:delta_t:T/2-delta_t];
% sample
ss2 = sinc(ts2/pi).^2; % Matlab defines sinc(x) as sin(pi*x)/(pi*x)!
%display
figure(1);plot(ts2,ss2,'r+');
%Fourier Domain
w0 = 2*pi/T;
W = 2*pi/delta_t;
fs2 = fft(ss2);
ws2 = [-W/2:w0:W/2-w0];
figure(2);plot(ws2,fftshift(abs(fs2)),'r+');
%
% Now Keep T as 10*pi, half delta_t 
T = 10*pi;
delta_t = pi/4;
ts3= [-T/2:delta_t:T/2-delta_t];
% sample
ss3 = sinc(ts3/pi).^2; % Matlab defines sinc(x) as sin(pi*x)/(pi*x)!
%display
figure(1);plot(ts3,ss3,'m+');
%Fourier Domain
w0 = 2*pi/T;
W = 2*pi/delta_t;
fs3 = fft(ss3);
ws3 = [-W/2:w0:W/2-w0];
figure(2);plot(ws3,0.5*fftshift(abs(fs3)),'m+'); % needs scaling WHY?