% Ricostruzione
% Author C. Prati
% March 24, 2025

clear
close all

%--------------------------------------------------
% 1
% Let us generate and plot  a "continuous" signal x=cos(2*pi*1*t)-sin(2*pi*0.3*t);
% in the time interval -T0<t<T0. (the infinitesimal dt=1/1000 seconds).
% 
%--------------------------------------------------

dt=1/1000;
T0=4;
t=(-T0:dt:T0);
x=cos(2*pi*1*t)-sin(2*pi*0.3*t);
figure(1)
plot(t,x)
L=length(t);
grid on
xlabel('time')
title('time continuous signal x(t)')
pause

%--------------------------------------------------
% 2
% Let us sample the signal x(t) with a sampling interval DT=1/4 
% obtaining the discrete signal xm.
% The samples of xm are superimposed to the plot of x(t)
%--------------------------------------------------

DT=1/10;
T=(-T0:DT:T0);
xm=x(1:DT/dt:L);
Lxm=length(xm);
figure(1)
hold on
stem(T,xm,'r','fill')
title('Original sampling')
hold off
pause


%--------------------------------------------------
% 3
% Let us interpolate M:1 the discrete signal xm with a sinc_N
% interpolator. The length of the sinc is set by the odd number LSINC
% LSINC is and odd number so that sinc goes from -LS to +LS where
% LS=(LSINC-1)/2
% The cumulative sum of the interpolator sinc((-LS:LS)/M) are progressively
% plotted until all the original samples of xm are explored. 
%---------------------------------------------

M=4;
LSINC=201;
LS=(LSINC-1)/2;
DT2=DT/M;
T2=(-T0-(LS*DT2):DT2:T0+(LS*DT2));
xn=zeros(1,length(T2));
xx=zeros(1,length(T2));
figure(2)
subplot(211)
axis([-8 8 -2 2])
grid on
title('Reconstruction sinc')
figure(2)
subplot(212)
axis([-8 8 -2 2])
grid on
title('Sum of the Reconstruction sinc')
pause
figure(2)
subplot(211)
    plot(t,x,'b')
    hold on
    plot(T,xm,'or')
    pause(.1)
for m=1:Lxm

    xx=zeros(1,length(T2));
    xx(M*(m-1)+1:M*(m-1)+LSINC)=xm(m)*sinc((-LS:LS)/M);
    subplot(211)
    plot(t,x,'b')
    hold on
    plot(T2,xx,'black')
    plot(T,xm,'xr')
    
    
    hold off
   pause(.5)
    subplot(212)
    xn(M*(m-1)+1:M*(m-1)+LSINC)=xn(M*(m-1)+1:M*(m-1)+LSINC)+xm(m)*sinc((-LS:LS)/M);
   
    plot(T2,xn,'black')
    grid on
    xlabel('time')
    title('Reconstructed')
    axis([ -8 8 -2.1 2.1])
    hold on
    plot(t,x,'b')
    plot(T,xm,'xr')
    hold off
pause(.5)
    
end



