% Campiona il segnale x(t)=sin(t*2*pi*3) con frequenza di campionamento
% a scelta ottenendo il segnale campionato xn.
% Mostra l'effetto dell'alias in frequenza quando non si rispetta il teorema 
% del campionamento che dice che fs>6.
% Si provi a vedere cosa succede se si utilizza fs=24, fs=12, fs=4, fs=3.


clear all, close all, clc

%% Parametri
BeginTime = -2;
EndTime  = 2;
TrueF = 440; % Frequenza Sinusoide
TrueF2 = 220;
Amplitude = 0.5;

flagPlay = 0;
flagPlayFiltered = 1;

fs = 2000; % Frequenza Campionamento

F_tot = lcm(TrueF, fs); % N campioni a secondo per approssimare il tempocontinuo
F_tot_Increment = 1;

F_tot = F_tot*F_tot_Increment;

ToPlotPeriods = 4;
ToPlotBegin = 0;

TranslationBand = 0;

%% Campionamento
dt = (F_tot)^-1;
t=(BeginTime:dt:EndTime); % Array dei tempi
k=length(t); % vedi significato inglese
x(1:k)=Amplitude*sin(2*pi*TrueF*t);% Costruisco un vettore con al suo interno i campioni della sinusoide
x = x + Amplitude*sin(2*pi*TrueF2*t);
figure(1), % Apro la figura 1
idxToPlotContinuo = (t>=ToPlotBegin & t<=ToPlotPeriods*(TrueF^-1));
plot(t(idxToPlotContinuo),x(idxToPlotContinuo))
grid on
xlabel('time')
title(sprintf('sinusoid frequency = %d Hz', TrueF))
axis([min(t(idxToPlotContinuo)), max(t(idxToPlotContinuo)), -1.2, 1.2])
hold on
pause
% fs=input('frequenza di campionamento fs? '); % Stampo a video e aspetto la risposta, che salvo in fs
%% Inserimento Filtraggi, vediamo a che punto arriviamo

Num = [1];
Den = [1];

T=1/fs;
nT=(BeginTime:T:EndTime); % Vettore Tempo Campionato
idxToPlotDiscreto = ((nT>=ToPlotBegin) & (nT<=ToPlotPeriods*(TrueF^-1)));
kk=length(nT);
xn=Amplitude*sin(2*pi*TrueF*nT);
xn=xn + Amplitude*sin(2*pi*TrueF2*nT);
figure(1)
stem(nT(idxToPlotDiscreto),xn(idxToPlotDiscreto),'r','filled')
pause

%% Ricostruzione Segnale

SincLengthInPeriods = 16;

SincTime = -(SincLengthInPeriods/2):dt*fs:(SincLengthInPeriods/2);
idxMaxSinc = (numel(SincTime)+1)/2;

SincFilter = sinc(SincTime) .* cos(2*pi*SincTime*TranslationBand);
SincTime = SincTime/fs;

figure(2), plot(SincTime, SincFilter)

ToReconstruct = zeros(numel(xn)*F_tot/fs, 1);
ToReconstruct(1:F_tot/fs:end) = xn;

Ricostruito = conv(ToReconstruct, SincFilter, 'same');

figure(1)
% plot(t(idxToPlotContinuo),Ricostruito(idxToPlotContinuo), 'b--')
plot([t(:);(t(end) + dt:dt:t(end) + dt*(F_tot/fs - 1)).'],Ricostruito(:), 'b--', 'LineWidth', 2)

H=audioplayer(Ricostruito, F_tot);

if(flagPlay)
    H.play
end
pause

%% Inserimento Filtraggi, vediamo a che punto arriviamo

% Num = [-0.5 1 -0.5]/2;
if(fs > 2* max([TrueF, TrueF2]))
    BandStop = TrueF/fs;% + 0.2; % Normalizzata
    Num = fir1(100, 2*[max(BandStop-0.05,0), min(BandStop+0.05,1)], 'stop');
    Den = [1];

    XFiltrato = filter(Num,Den,xn);
    ToReconstruct = zeros(numel(XFiltrato)*F_tot/fs, 1);
    ToReconstruct(1:F_tot/fs:end) = XFiltrato;

    RicostruitoFiltrato = conv(ToReconstruct, SincFilter, 'same');

    figure(3), % Apro la figura 3
    hold on, grid on,
    plot(t(idxToPlotContinuo),x(idxToPlotContinuo)), hold on, grid on
    stem(nT(idxToPlotDiscreto),xn(idxToPlotDiscreto),'r','filled')
    title('Segnale campionato')

    hold off
    axis([min(t(idxToPlotContinuo)), max(t(idxToPlotContinuo)), -1.2, 1.2])
    pause
    
    plot(t(idxToPlotContinuo),x(idxToPlotContinuo)), hold on, grid on    
    stem(nT(idxToPlotDiscreto),XFiltrato(idxToPlotDiscreto),'r','filled'), hold on, grid on
    idxToPlotContinuo = (t>=ToPlotBegin & t<=ToPlotPeriods*(TrueF^-1));
    axis([min(t(idxToPlotContinuo)), max(t(idxToPlotContinuo)), -1.2, 1.2])
    title('Segnale Filtrato')
    pause
    plot([t(:);(t(end) + dt:dt:t(end) + dt*(F_tot/fs - 1)).'],RicostruitoFiltrato(:), 'b--', 'LineWidth', 2)
    title('Segnale filtrato e ricostruito')    

    if(flagPlayFiltered)
        H1=audioplayer(RicostruitoFiltrato, F_tot);
        H1.play;
    end
    
end