% Recovering a signal from sampled data using the sinc reconstruction filter

clc             % clear the command window
clear           % clear all variables
close all       % close any existing windows
format compact	% prevent extra blank lines in the output

% Original signal: 15 periods of a 1 Hz sinewave, plotted with fine resolution
T = 1;
f = 1/T;
w = 2*pi/T;
Ts = 15*T;
dt = 0.01;
t = 0 : dt : Ts;
F = sin(w*t);
plot (t, F);
xlabel('t');
ylabel('sin(2\pit)');

% Signal digitally sampled at just over the Nyquist frequency (twice the signal frequency)
fs = 2.1*f;
dts = 1/fs;
Ns = Ts/dts;
i = 0 : Ns;
t_sampled = i * dts;
F_sampled = sin(w*t_sampled);
figure;
plot (t_sampled, F_sampled);
xlabel('t');
ylabel('F');

% sinc function
figure;
t = -6 : 0.1 : 6;
plot (t, sinc(t));
title('sin(\pit)/(\pit)');
xlabel('t');
ylabel('sinc(\pit)');
grid on;

% Signal recreated from the sampled data using the sinc reconstruction filter
% file: sinc_reconstruct.m
% function [f] = sinc_reconstruct(t, t_sampled, F_sampled, Ns, dts)
% % Function to calculate a sinc resconstruction filter
%   f = 0;
%   for i = 1 : Ns+1
%     f = f + F_sampled(i)*sinc((t-t_sampled(i))/dts);
%   end
% end
N = Ts/dt;
i = 0 : N;
t_reconstruct = i*dt;
F_reconstruct = sinc_reconstruct(t_reconstruct, t_sampled, F_sampled, Ns, dts);
figure;
plot (t_reconstruct, F_reconstruct);
xlabel('t');
ylabel('f');
hold on;
plot (t_reconstruct, sin(w*t_reconstruct), 'k:');
legend('reconstructed', 'original');

disp 'Hit Enter to close all windows and quit'
pause
close all