%  Sinc function interpolation example upsamples a rectangular pulse
%  through interpolation.  In the final plot note the effects of the 
%  bandlimited nature of the sinc function interpolation (recall it
%  is the FT pair with the ideal low pass filter).  The overshoot
%  and undershoot around the edges is sometimes referred to as Gibbs
%  phenomenon.
%
%    written by Kevin D. Donohue (donohue@engr.uky.edu)  June 2002
%  This script demonstrates upsampling using the sinc function/ideal filter


%  Create a simple set of samples for interpolation
g = [0 0 0 0 0 0 0  0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0];   %  Square wave
%g = [0 0 0 1 2 3 4 5 0 -5 -4 -3 -2 -1 0 0 0];  %  Triangle wave - odd symetric
%g = [0 0 0 1 2 3 4 5 6 5 4 3 2 1 0 0 0];  %  Triangle wave - even symetric


%  Define original sampling rate and sampling rate of upsampled signal 
fs = 1;   %  Original sampling frequency
taxis = [0:length(g)-1]/fs;  %  Original time axis
fsup = fs*100;   %  Upsampled frequency
taxup = min(taxis)+[0:floor((length(g))*fsup/fs)-1]/fsup;  % Upsampled time axis
upsig = zeros(size(taxup));   %  Initialize array to accumulated interplolated points

%  Big loop steps through each new point to be interpolated
for k = 1:length(taxup)
    % Inner loop weights all original signal points with sinc values to get
    % and adds up the products to get the new point 
    for n = 1:length(g)
        upsig(k) = upsig(k)+g(n)*sinc(fs*(taxup(k)-taxis(n)));
    end
end
% Plot original and upsmpled signal
plot(taxis,g,'r',taxup,upsig,'b')
xlabel('seconds')
ylabel('Amplitude')
title('Original signal in red, upsampled in blue')