function shanwhit
%This worksheet demonstrates the classical Shannon-Whittaker formula.
%Choose an example 1, 2 or 3. In 1 and 3 there are parameters you can
%adjust to control the effective bandlimit of the function.
%
%We use 2^(ks+1)+1 samples
ks=6;
%We interpolate onto a grid with 2^ki+1 samples
ki=10;
%
%
%Choose a function to sample and interpolate
example = 3;
if example == 1
    %l controls the smoothness. This can be any non-negative number
    %Try numbers between 0 and 1 to get less smooth functions, that is less
    %bandlimited functions.  Try 1/4, 1/2, 1, 2 (for example).
    l=1/2;
    f = @(x) (1-x.*x).^l.*(abs(x)<= 1);
elseif example == 2
    f = @(x) (abs(x)<= 0.5);
elseif example == 3
    p=2;
    f = @(x) (abs(x)<=1).*cos((2*p+1)*pi*x/2);
end

%
%Grid for sampling
%
xs=linspace(-2,2,2^(ks+1)+1);
%Sample spacing
d=2^(1-ks);
%Grid for interpolation
xi=linspace(-2,2,2^ki+1);
%
%Build interprolation matrix. This is a (2^ki+1)x(2^(ks+1)+1) matrix
for j=1:2^(ks+1)+1
    intm(1:2^ki+1,j)=sinc((xi/d)+2^ks-j+1);
end
intm=intm;
%Sample f on the coarse grid
fs=f(xs)';
%Interpolate from the samples
fi=intm*fs;
%Sample f on the fine grid
ff=f(xi);
%Compare the results
figure
subplot(1,2,1)
plot(xi,fi,'r',xi,ff,'g')
title('Shannon Whittaker interpolant (red)')
%Let's also look at the magnitude of the Fourier transform in a semilogy
%plot
ffi=fft(fi);
subplot(1,2,2)
semilogy(abs(ffi))
title('log_{10} modulus of FFT')
xlim([1,length(fi)]);

