function sampling_demo(nu)
% SAMPLING_DEMO illustrates the sampling theorem for sinc^2(t),
% which is a band-limited function. The parameter nu is the 
% Nyquist rate.
%
% Francis J. Narcowich--6/6/00 (Revised 3/25/04)
%

if nargin==0,
nu=3; 
end
t_samp=(-60:60)/nu;
s=sinc(t_samp).^2; 
%t=-5:0.01:5; 
t=linspace(-5,5,4096);
y=sampling_thm(t,nu,s);

%[yhat,omega]=nstdfft(y,-5,5);
[yhat,omega]=ft_approx(y,-5,5);

indx=find(abs(omega)<2);

%[sig_hat,omega_sig]=nstdfft(sinc(t).^2,-5,5); 
[sig_hat,omega_sig]=ft_approx(sinc(t).^2,-5,5);

yhat2=zeros(size(omega_sig));

% Find the indices where omega_sig is between -nu/2 and nu/2.
indx0=find(abs(omega_sig)<=nu/2); 

% This is the integer period for the signal.
period_sig=indx0(end)-indx0(1); 

 % Form aliased FT for signal.
yhat2(indx0)=sig_hat(indx0)+sig_hat(indx0-period_sig)+ ...
    sig_hat(indx0+period_sig);

if nu<1,
    yhat2=yhat;
end


subplot(2,2,1),plot(t,sinc(t).^2), axis tight, grid, title('sinc^2(t)')
     subplot(2,2,2),plot(t,y-sinc(t).^2,'g'), grid, title('Error')
     subplot(2,2,3),plot(t,y,'r'), axis tight,grid, title('Reconstruction')
     subplot(2,2,4),plot(omega(indx),real(yhat2(indx)),'ko',...
                         omega(indx),real(sig_hat(indx))), ... 
            grid, legend('FT/Rec.','FT/orig.'), title('Fourier Transforms')

function y=sampling_thm(time,frequency,samples)
nu=frequency;
t=nu*time;
s=samples; L=length(s);
N=floor(L/2);
y=s*sinc(ones(L,1)*t - ((-N):(L-N-1))'*ones(size(t)));

function y=sinc(t)
y=ones(size(t));
nn=find(t~=0);
y(nn)=sin(pi*t(nn))./(pi*t(nn));


function [yhat,omega]=ft_approx(y,a,b)
N=floor(length(y)/2);
dt=(b-a)/(2*N);
del_omega=1/(b-a);
omega=((-N):(N-1))*del_omega;
phase_factor=exp(-i*a*omega*2*pi);
yhat=phase_factor.*fftshift(fft(y))*(dt/sqrt(2*pi));