%In punti:
clear all; close all; clc;
%% 1a. creare componenti di segnale utile, interferenza e rumore
% frequenza della sinusoide: a scelta, compresa tra 200 Hz e 1 kHz;
% frequenza del disturbo: a scelta, compresa tra 4 kHz e 7 kHz;
% ampiezza della sinusoide del disturbo: a scelta, compresa tra 10% e 100%
% del segnale utile;
% parametri delle componenti del segnale
f0=500; %Hz, frequenza del segnale utile
V0=1;  %Volt, ampiezza segnale utile
f0_dist=5000; %Hz, frequenza del segnale di disturbo sinusoidale
V0_dist=0.2*V0; %Volt, ampiezza del disturbo
SNR=60; % scelta del rapporto segnale rumore inteso rispetto al segnale utile
P_noise=((V0/sqrt(2)).^2)/SNR;

% parametri di acquisizione
Tw=0.05; %s, finestra di osservazione
fs=100000; %Hz
Ts=1/fs; %s
N_sample=round(Tw/Ts);
t_sample=[0:1:N_sample-1]*Ts;

% simulazione delle componenti del segnale
V_sample=V0*sin(2*pi*f0*t_sample);
V_sample_dist=V0_dist*sin(2*pi*f0_dist*t_sample);
noise=sqrt(P_noise)*randn(1,N_sample);

 %% 1b (alternativa a 1a) volendo sfruttare una funzione, 
% % vedere in fondo al codice la funzione "sinusoide"
% % parametri delle componenti del segnale e del campionamento
% f0=500; %Hz, frequenza del segnale utile
% V0=1;  %Volt, ampiezza segnale utile
% f0_dist=5000; %Hz, frequenza del segnale di disturbo sinusoidale
% V0_dist=0.2*V0; %Volt, ampiezza del disturbo
% Tw=0.1; %s, finestra di osservazione
% fs=200000; %Hz
% [N, t_sample,V_sample]=sinusoide(V0,f0,fs,Tw);
% [N, t_sample,V_sample_dist]=sinusoide(V0_dist,f0,fs,Tw);
% 
% SNR=60; % scelta del rapporto segnale rumore inteso rispetto al segnale utile
% P_noise=((V0/sqrt(2)).^2)/SNR;
% noise=sqrt(P_noise)*randn(1,N);

%% indipendentemente da che si scelga 1a o 1b,
% creare un grafico in cui siano sovrapposti, nell’ordine: 
%la somma di segnale utile, disturbo e rumore;
%la somma di segnale utile e disturbo;
%il solo segnale utile.

V_sample_sum_noise=V_sample+V_sample_dist+noise;
V_sample_sum=V_sample+V_sample_dist;

figure
plot(t_sample, V_sample_sum_noise, 'b');
hold on
plot(t_sample, V_sample_sum,'r');
hold on
plot(t_sample, V_sample,'k');
legend('Segnale con disturbo e interferenza','Segnale con interferenza','Segnale utile');
xlabel('Time (s)');
ylabel('Amplitude (V)');
hold off

%% 2. creare la function adc ideal() che quantizzi il segnale;
%vedi adc_ideal function in fondo al codice, righe identiche a quanto fatto
%la scorsa esercitazione
figure
N_bit=10;
V_fs=1;
[V_sample_q,Delta] = quantizzatore(N_bit,V_fs, V_sample_sum_noise);
hold on 
plot(t_sample, V_sample_sum_noise, 'r');
plot(t_sample,V_sample_q);
hold off
e_q=V_sample-V_sample_q;

%% 3. Determinare le specifiche del filtro 
% obiettivo: tenere solo il segnale utile e eliminare l'interferenza
% sinusoidale ad alta frequenza e il rumore a larga banda.
% Le specifiche del filtro una volta definite dovranno essere tradotte in  
% valori da assegnare ad altrettante variabili MatLab (indicate come Wp, Ws,
% Rs e Rp) che vengono passate alle funzioni per la sintesi del filtro.

% FREQUENZA DI BANDA PASSANTE e DI BANDA OSCURA
BP=1000;
BO=3000;
% ATTENZIONE: le frequenze devono essere normalizzate rispetto alla 
% frequenza di Nyquist, Fs/2. I valori di Wp = BP/(Fs/2) e Ws = BO/(Fs/2) 
% sono valori compresi tra 0 e 1, con Wp < Ws;
Wp=BP/(fs/2);
Ws=BO/(fs/2);

% ATTENUAZIONE (A)
A=100;
% per passare il valore a buttord trasformarlo in dB:
Rs=20*log10(A);

% ONDULAZIONE (ripple)IN BANDA PASSANTE (R)
R=0.05; % in percentuale (0.05 sta per il 5%)
% per passare il valore a buttord trasformo in Rp:
R1=20*log10(1+R);
R2=20*log10(1-R);
Rp=(R1-R2)/2;

%% 4. Calcolare i coefficienti del filtro di Butterworth
% Primo step: trovo ordine e banda passante del filtro che sia in grado di
% rispettare le specifiche richieste
[n,Wn] = buttord(Wp,Ws,Rp,Rs);

Wn_Hz=Wn*(fs/2); %rappresenta la banda passante

% Secondo step: trovo i vettori dei coefficienti della parte MA e AR. In
% ingresso a butter inserisco ordine e banda ottenuti con buttord e il tipo
% di filtro che si desidera tra "low" per bassa-basso, "high" per
% passa-alto e "bandpass" per passa banda. In questo caso essendo il
% disturbo a freq più alta scegliamo un passa basso
 [b,a] = butter(n,Wn,'low');

% Nei casi in cui dovessero emergere limitazioni dovute ad arrotondamenti,
% può essere utile sfruttare la forma
% [z,p,k] = butter(n,Wn,'low');
% [b,a] = zp2tf(z,p,k);
% % sos = zp2sos(z,p,k);
% hfvt = fvtool(b,a,sos,'FrequencyScale','log');
% legend(hfvt,'TF Design','ZPK Design')

%% 5. Verificare la risposta in frequenza del filtro 
% funzione freqz di matlab che, dati i coeffcienti del filtro contenuti nei 
% vettori a e b e la frequenza di campionamento Fs, calcola la risposta in 
% frequenza su K campioni (dove K deve essere una potenza di 2, p. es., 
% 1024). In uscita sono forniti i vettori f ed H, cioè l'asse delle frequenze
% e la risposta in frequenza calcolata.

K=1024;
[H, f] = freqz(b, a, K, fs);
figure
plot(f,20*log10(abs(H)));
xlabel('freq in Hz');
ylabel('gain (dB)');
%axis([0 1 -100 0]);

%% 6. Filtrare il segnale;
y = filter(b, a, V_sample_q);
% in input segnale quantizzato con componenti di rumore

%% 7. Confrontare i segnali nel dominio del tempo pre- e post-filtraggio.
% creando un grafico in cui siano sovrapposti:
% 1. il segnale ottenuto all'uscita del quantizzatore;
% 2. il solo segnale utile, creato all'inizio della simulazione;
% 3. il segnale filtrato,per verificare l'effetto del filtro e la qualità del 
% segnale ottenuto.
figure
plot(t_sample,V_sample_q,'c');
hold on
plot(t_sample,y,'r');
hold on
plot(t_sample,V_sample,'b');
legend('segnale con disturbo quantizzato','segnale filtrato','segnale utile');
hold off

%% 8. Valutazione al variare di SNR ed Nbit,Si puo quindi ripetere la 
% simulazione per diversi valori di Nbit, per ciascuno dei quali si possono
%valutare diversi valori di SNR, ad esempio: SNR circa 6*Nbit (es. impostare
% SNR a 48 per 8 bit scelti) o SNR<<6*Nbit (es. SNR=10 con 8 bit).
V_RMS=V0/sqrt(2);
ratio=(V_RMS^2)/((Delta^2)/12);
SNRQ_dB=10*log10(ratio)
ratio_check=6.02*N_bit;
matlab_snr_dB=snr(V_sample_q)
snr_dB=10*log10(SNR)

% N.B. attenzione a non confondere l'SNR legato esclusivamente alla 
%quantizzazione (SNRQ) e l'SNR complessivo che include anche il contributo
%non solo del rumore di quantizzazione ma anche dei disturbi.
%Se aumentate il numero di bit aumenterà il valore di SNRQ, ma mantenendo
%la stessa ampiezza di rumore aggiuntivo l'SNR resterà invariato


%% FUNZIONI

%% Function1 

function [N, t_sample,V_sample]=sinusoide(Amp,f0,fs,Tw)
Ts=1/fs; %s
N=round(Tw/Ts);
t_sample=[0:1:N-1]*Ts;
V_sample=Amp*sin(2*pi*f0*t_sample);
end

%% Function2

function [V_q,Delta] = quantizzatore(N_bit, V_FS, V)
N_liv=2^N_bit;                %numero livelli disponibili
Delta=2*V_FS/(N_liv-1);       %passo di quantizzazione
lev=round(V/Delta);           %vettore di livelli
for i=1:length(lev)           %ciclo for che descrive saturazione al di fuori del campo di ingresso
    if lev(i)> 2^(N_bit-1)-1
       lev(i)=2^(N_bit-1)-1;
    end
    if lev(i)< -2^(N_bit-1)
       lev(i)=-2^(N_bit-1);
    end
end
V_q=lev*Delta;
errore = V_q - V;
end