% Vogliamo qui raddoppiare la frequenza di campionamento di un segnale
% campionato con Fs bassa (<= 11000Hz) attraverso la tecnica dell'interpolazione
% la quale richiede l'inserzione di valori nulli tra i campioni,
% seguìta da un il filtraggio passabasso.

% caricare qui il file audio che si desidera e inizializzare la variabile Fs
load handel;	% Fs= 8192 Hz; campioni nel vettore y
tfin= length(y)/Fs;	% durata brano
tvec= linspace(0,tfin,length(y));	% vettore valori temporali [s]

% inserzione degli zeri: sono campioni fittizi che raddoppiano artificiosamente
% la frequenza di campionamento
r= 2;	% fattore di upsample
yup= upsample(y,r);
Fsup= r*Fs;	% freq. di campionamento per il segnale yup
tupfin= length(yup)/Fsup;	% durata invariata (=tfin)
tupvec= linspace(0,tupfin,length(yup));	% [s]

% parametri per la visualizzazione di segnali e spettri
t1= 2.6; t2= 2.605;% [s]
vx= 2.0;	% ampiezza segnali
VX= 1500;	% ampiezza spettri

figure(1); clf;
stem(tvec,y,'b');
xlabel('t [s]'); ylabel('y'); title('Segnale originario');
axis([t1 t2 -vx vx]);
figure(2); clf;
stem(tupvec,yup,'r');
xlabel('t [s]'); ylabel('y_{up}'); title('Segnale upsampled');
axis([t1 t2 -vx vx]);
disp('premi un tasto per ascoltare il segnale originario'); pause;
soundsc(y,Fs);

% Effetti sullo spettro dell'inserzione degli zeri
% compressione dello spettro periodico
fvec= linspace(-Fs/2,Fs/2,length(y));	% vettore frequenze
fupvec= linspace(-Fsup/2,Fsup/2,length(yup));	% freq. di Nyquist raddoppiata
figure(3); clf;
plot(fvec,fftshift(abs(fft(y))),'b');
xlabel('f [Hz]'); ylabel('Y'); title('Spettro di ampiezza Segnale originario');
axis([-Fsup/2 Fsup/2 0 VX]);
figure(4); clf;
plot(fupvec,fftshift(abs(fft(yup))),'r');
xlabel('f [Hz]'); ylabel('Y_{up}'); title('Spettro di ampiezza Segnale upsampled');
axis([-Fsup/2 Fsup/2 0 VX]);
% NOTA: lo spettro di un segnale campionato è proporzionale a Fs
% che noi stiamo raddoppiando...
% QUINDI dovremo inserire un fattore 2 di guadagno
% poiché non è corretto avere gli stessi valori nello spettro d'ampiezza
disp('premi un tasto per ascoltare il segnalee upsampled'); pause;
soundsc(yup,Fsup);


% Progetto del filtro per l'interpolazione del segnale:
% passabasso di guadagno 2 e banda passante (monolatera) Fs/2
% Metodo della finestratura
fcut= Fs/2;	% frequenza di taglio del filtro
M= 32;	% ordine del filtro (tipo I): scelto arbitrariamente...

% FUNZIONE fir1(M,wcut,winvec);	% di default usa winvec= hamming(M+1);
% h= r*fir1(M,fcut/(Fsup/2));	% NOTA: freq. di taglio normalizzata alla freq. di Nyquist
% che è equivalente a:
% h= r*fir1(M,fcut/(Fsup/2), hamming(M+1));

% ALTERNATIVA: CALCOLIAMO h ESPLICITAMENTE. CONTROLLARE CHE IL RISULTATO E' LO STESSO!
nvec= (-M/2:M/2);	% vettore simmetrico di M+1 punti (M pari)
% hd= r*sin(2*pi*(fcut/Fsup)*nvec)./(pi*nvec);	% risposta impulsiva passabasso desiderato,
% esprimibile anche come...
hd= r*(2*(fcut/Fsup)*sinc(2*(fcut/Fsup)*nvec));	   % così si evita l'indeterminazione 0/0
% così com'è, hd è troncata con finestra rettangolare
wh= hamming(M+1)';	% finestra di Hamming (simmetrica); trasposto poiché è un vettore colonna
h= hd.*wh;


% visualizzazione di entrambe i filtri (finestra rect + finestra Hamming)
fvtool(hd,1,h,1);	% confronto Fin. Rect., Fin. Hamming
% click destro sull'asse verticale x selezionare "Magnitude"
% - controllo guadagno =r
% - controllare che la banda di transizione (ma max a min ripple) sia
%        = 8/M (*pi) x Hamming; = 4/(M+1) (*pi) x rect.
% In particolare, x Hamming si ha
wp= 0.375; ws= 0.609;	% si useranno dopo
% - controllo deviazione 20*log10((Amaxpass-r)/r)=20*log10(Aminstop/r) ~= -21 (rett.) o -53 (Hamming)
% In particolare, x Hamming si ha
deltap= 0.0037; deltas= 0.00372;	% si useranno dopo
% - controllo fase lineare e ritardo di gruppo
% Si può visualizzare l'ascissa in Hz
% inserendo la freq. di campionamento Fsup (dal Menu Insert)
% e facendo click destro sull'asse orizzontale x selezionare "Linear Frequency"
disp('premi un tasto per proseguire'); pause;

% Quanta complessità avremmo risparmiato usando una finestra di Kaiser?
[Mk, Wn, beta, ftype]= kaiserord([wp ws], [1 0], [deltap/r deltas/r]);
hk= r*fir1(Mk,Wn,ftype,kaiser(Mk+1,beta));	% includo il guadagno r
sprintf('Con finestra Kaiser avresti avuto ordine %d anziché %d', Mk,M)
fvtool(hk,1,h,1);	% confronto Fin. Kaiser, Fin. Hamming
% NOTARE che sono praticamente identiche, a fronte di 4 coeffricienti in meno
disp('premi un tasto per proseguire'); pause;

% segnale interpolato (Hamming)
yint= filter(h,1,yup);
figure(6); clf;
stem(tupvec,yint,'k');
xlabel('t [s]'); ylabel('y_{int}'); title('Segnale interpolato (Fin. Hamming)');
axis([t1 t2 -vx vx]);
% NOTARE il ritardo di M/2 campioni nel tempo
figure(7); clf;
plot(fupvec,fftshift(abs(fft(yint))),'k');
xlabel('f [Hz]'); ylabel('Y_{int}'); title('Spettro di ampiezza Segnale interpolato');
axis([-Fsup/2 Fsup/2 0 VX]);
disp('premi un tasto per ascoltare il segnale interpolato (Hamming)'); pause;
soundsc(yint,Fsup);

% segnale interpolato (finestra Rect)
yintR= filter(hd,1,yup);
figure(8); clf;
stem(tupvec,yintR,'g');
xlabel('t [s]'); ylabel('y_{int}'); title('Segnale interpolato (Fin. Rect.)');
axis([t1 t2 -vx vx]);
% NOTARE il ritardo di M/2 campioni nel tempo
figure(9); clf;
plot(fupvec,fftshift(abs(fft(yintR))),'g');
xlabel('f [Hz]'); ylabel('Y_{int}'); title('Spettro di ampiezza Segnale interpolato');
axis([-Fsup/2 Fsup/2 0 VX]);
% NOTARE il rigonfiamento dello spettro nella banda di transizione, rispetto alla finestra rect.
disp('premi un tasto per ascoltare il segnale interpolato (Rect.)'); pause;
soundsc(yintR,Fsup);


% Matlab esegui l'interpolazione di un fattore intero r automaticamente con la funzione
[yintM,b]= interp(y,r);	% in b restituisce i coeff. del filtro
figure(10); clf;
stem(tupvec,yintM,'m');
xlabel('t [s]'); ylabel('y_{intM}'); title('Segnale interpolato da Matlab');
axis([t1 t2 -vx vx]);
% NOTA non c'è ritardo temporale poiché (presumibilmente) Matlab esegue un "Filtraggio doppio a fase nulla"
figure(11); clf;
plot(fupvec,fftshift(abs(fft(yintM))),'m');
xlabel('f [Hz]'); ylabel('Y_{int}'); title('Spettro di ampiezza Segnale interpolato');
axis([-Fsup/2 Fsup/2 0 VX]);
% NOTARE le "code" dello spettro, che sono più estese....infatti (vedi dopo)
disp('premi un tasto per ascoltare il segnale interpolato da Matlab'); pause;
soundsc(yintM,Fsup);

% caratteristiche del filtro usato da "interp"
fvtool(b,1,h,1);	% confronto filtro di "Interp" - filtro Fin. Hamming
% Il filtro usato da "interp" ha banda di transizione + larga, che
% lascia sopravvivere delle "code di aliasing" più grandi, come si vede nello spettro