function aliasing(example,m,k)
%This worksheet illustrates several phenomena connected to sampling and
%aliasing for Fourier series. 
%example = 1,2,3 chooses one of the examples below
%m is an integer between 1 and 10, it determines the sample spacing:
%Delta x=16*2^m/N
%k is an optional parameter which determines the smoothness of the
%functions in examples 2 and 3, larger k means smoother examples.
%
%
M=2^m;
%To effect sampling we work directly in the Fourier domain.
%
if nargin <3
    k=2
end
%The maximum frequency sampled is about N/4*pi
N=8192*2;
%We use N/M frequecies in the reconstruction:
%M=4;
%This means that the effective spacing is 16*M/N
smp_spcing=16*M/N
smp_rate=1/smp_spcing
%
%The "fully" sampled sample spacing:
x=linspace(-4,4,N);
%This is the sampling that we end up using in the reconstructions:
xl=linspace(-4,4,N/M);
%
%example = 1;
if example == 1
    %Characteristic function, which is very unbandlimited
h=(abs(x)<1);
hsl=(abs(xl)<1);
elseif example == 2
    %Smoother function with bounded support. The larger k is the faster the
    %Fourier coefficients decay.
    %k=4;
    h=(abs(x)<1).*(1-x.*x).^k;
    hsl=(abs(xl)<1).*(1-xl.*xl).^k;
elseif example == 3
    %Here we modulate the previous and use sampling to demodulate.
     %k=2;
     %Modulation frequency
     mod=100;
    h=(abs(x)<1).*(1-x.*x).^k;
    %We modulate it here
    h=real(h.*exp(mod*pi*i*x));
    hsl=(abs(xl)<1).*(1-xl.*xl).^k;
     %We modulate it here
     hsl=real(hsl.*exp(mod*pi*i*xl));
end
figure, 
plot(x,h)
title('Original Function')
Fh=fft(h);
%We seperate out the low frequencies:
%In the ordering used by the FFT these are the samples (in the frequency domain)
%with indices (1:N/M) along with (N-N/M+1:N)
Fhl=[Fh(1:N/(2*M)) Fh((N-N/(2*M)+1):N)];
%This gives the low pass filtered version of h:
hl=real(ifft(Fhl)/M);
figure,
subplot(1,4,1)
plot(xl,hl)
title('Low frequency approximation')
subplot(1,4,2)
plot(xl,hsl-hl)
title('Low frequency approx error')
%We now collect the high frequency samples into aliased frequencies
FhAl=zeros(1,N/(2*M));
FhAh=zeros(1,N/(2*M));
for j=1:M-1
    %Here we shift the high frequencies to lower frequencies and add them
    %up as occurs in aliasing:
    FhAl=Fh((1+j*N/(2*M)):((1+j)*N/(2*M)))+FhAl;
    FhAh=Fh((N/2+(j-1)*(N/(2*M))+1):(N/2+j*(N/(2*M))))+FhAh;
end
%We compute the pure aliasing contribution:
ha=(ifft([FhAl FhAh])/M);
subplot(1,4,3)
plot(xl,real(ha),xl,imag(ha))
title('Pure aliasing artifact')
%Finally we show the reconstruction that results from the Shannon-Whitakker
%with the sample spacing 16*M/N:
subplot(1,4,4)
plot(xl,real(ha)+hl)
title('Full aliased reconstruction')

