% Questa esercitazione ? speculare rispetto a quella sulla finestratura.
% L'obiettivo qui ? ridurre di un fattore r=3 la frequenza di campionamento di un segnale
% attraverso la tecnica della "decimazione"
% la quale richiede lo scarto di campioni (si mantiene solo 1 campione ogni r)
% e deve essere preceduta da un filtraggio passabasso per evitare aliasing.
% Useremo un filtro equiripple progettato con l'Alg. di Parks-McClellan (funzione "remez")
clear all;

% 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]

% scarto dei campioni, se ne tiene solo 1 su 3, a partire dal primo
r= 2;	% fattore di downsample
ydn= downsample(y,r);
Fsdn= Fs/r;	% freq. di campionamento per il segnale ydn
tdnfin= length(ydn)/Fsdn;	% durata invariata (=tfin)
tdnvec= linspace(0,tdnfin,length(ydn));	% [s]

% parametri per la visualizzazione di segnali e spettri
t1= 2.6; t2= 2.605;% [s]
vx= 2.0;	% ampiezza segnali
VX= 800;	% 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(tdnvec,ydn,'r');
xlabel('t [s]'); ylabel('y_{dn}'); title('Segnale downsampled');
axis([t1 t2 -vx vx]);
disp('premi un tasto per ascoltare il segnale originario'); pause;
soundsc(y,Fs);

% Effetti sullo spettro dello scarto di campioni
% Aliasing: ? come aver sottocampionato il segnale analogico originario
fvec= linspace(-Fs/2,Fs/2,length(y));	% vettore frequenze
fdnvec= linspace(-Fsdn/2,Fsdn/2,length(ydn));	% freq. di Nyquist divisa per r
figure(3); clf;
plot(fvec,fftshift(abs(fft(y))),'b');
xlabel('f [Hz]'); ylabel('Y'); title('Spettro di ampiezza Segnale originario');
axis([-Fs/2 Fs/2 0 VX]);
figure(4); clf;
plot(fdnvec,fftshift(abs(fft(ydn))),'r');
xlabel('f [Hz]'); ylabel('Y_{dn}'); title('Spettro di ampiezza Segnale downsampled');
axis([-Fs/2 Fs/2 0 VX]);
% NOTARE il ripiegamento dello spettro (aliasing)
%        e che non ci sono campioni oltre Fsdn/2
% NOTA: lo spettro di un segnale campionato ? proporzionale a Fs
% che qui ? dimezzata e cos? deve essere
% Il fattore di guadagno del filtro antialiasing dovr? dunque essere 1
disp('premi un tasto per ascoltare il segnalee downsampled'); pause;
soundsc(ydn,Fsdn);

% Dobbiamo far precedere il downsample da un filtro antialiasing

% Progetto del filtro per la decimazione del segnale:
% passabasso di guadagno 1 e banda passante (monolatera) Fs/2/r
% Metodo Equiripple
fcut= (Fs/2)/r;	% frequenza di taglio del filtro
% scegliamo una banda di transizione del 20% della banda passante
fp= fcut*0.9;  fs= fcut*1.1;
ap= 1;  as= 0;	% guadagni in banda passante e stopband
deltap= 0.001; deltas= 0.01;	% massime deviazioni dal guadagno ideale
[M,f,m,w]= remezord([fp fs], [ap as], [deltap deltas], Fs);	% frequenze non normalizzate poich? fornisco Fs
[h, delta]= remez(M,f,m,w);	% calcolo filtro e massimo ripple
sprintf('Freq. di campionamento %f Hz; fp= %f Hz; fs= %f Hz',Fs,fp,fs)
sprintf('Ordine filtro= %d; max ripple da specifiche= %f; max ripple raggiunto= %f',M,max(deltap,deltas),delta)
% visualizzazione del filtro
fvtool(h,1);
% click destro sull'asse verticale x selezionare "Magnitude"
% - controllo guadagno =1
% - controllo rispetto deltap e deltas assegnati; max ripple delta raggiunto
% - controllo della propriet? equiripple e delle alternanze
% - controllo della derivata nulla (punto estremo) in f=0 o Fs/2 (se M ? pari)
% - conteggio numero di alternanze = M/2+2 (M= ordine del filtro)
% visualizzare l'ascissa in Hz inserendo la freq. di campionamento Fs (dal Menu Insert)
% e facendo click destro sull'asse orizzontale x selezionare "Linear Frequency"
% - controllare la banda di transizione: punti estremi in fp e fs
%      guadagno 1-deltap  in fp;  guadagno deltas   in fs
% - controllo fase lineare e ritardo di gruppo
disp('premi un tasto per proseguire'); pause;


% Quanta complessit? avremmo avuto usando una finestra di Kaiser?
wp= fp/(Fs/2);  ws= fs/(Fs/2);
[Mk, Wn, beta, ftype]= kaiserord([wp ws], [ap as], [deltap deltas]);
hk= fir1(Mk,Wn,ftype,kaiser(Mk+1,beta));	% includo il guadagno r
sprintf('Con finestra Kaiser avresti avuto ordine %d ; col metodo ottimo %d', Mk,M)
fvtool(h,1,hk,1);	% confronto Ottimo/Kaiser
% NOTARE che Kaiser ha il massimo ripple vicino alla transizione e che il ripple ? uguale nelle 2 bande
%        dunque ? assolutamente sovradimensionato, nella maggior parte delle frequenze
disp('premi un tasto per proseguire'); pause;


% segnale filtrato (filtro ottimo) e poi decimato
yfil= filter(h,1,y);
figure(7); clf;
stem(tvec,yfil,'k');
xlabel('t [s]'); ylabel('y_{fil}'); title('Segnale filtrato (Filtro Ottimo)');
axis([t1+M/2/Fs t2+M/2/Fs -vx vx]);
% NOTARE il ritardo di M/2 campioni nel tempo
% cambio gli estremi temporali poich? M ? grande e non si vedrebbe il confronto con Fig. 1
figure(8); clf;
plot(fvec,fftshift(abs(fft(yfil))),'k');
xlabel('f [Hz]'); ylabel('Y_{fil}'); title('Spettro di ampiezza Segnale filtrato');
axis([-Fs/2 Fs/2 0 VX]);
disp('premi un tasto per ascoltare il segnale filtrato (Filtro Ottimo) con Fs originaria'); pause;
soundsc(yfil,Fs);

% segnale decimato
ydec= downsample(yfil,r);
figure(9); clf;
stem(tdnvec,ydec,'g');
xlabel('t [s]'); ylabel('y_{dec}'); title('Segnale decimato');
axis([t1+M/2/Fs t2+M/2/Fs -vx vx]);
% NOTARE il ritardo di ? stato introdotto prima, nel filtraggio
figure(10); clf;
plot(fdnvec,fftshift(abs(fft(ydec))),'g');
xlabel('f [Hz]'); ylabel('Y_{dec}'); title('Spettro di ampiezza Segnale decimato');
axis([-Fs/2 Fs/2 0 VX]);
% NOTARE che sono spariti campioni FFT dove prima (yfil) c'erano campioni nulli
% ma che non c'? pi? l'aliasing di Fig. 4
disp('premi un tasto per ascoltare il segnale decimato, con Fs ridotta'); pause;
soundsc(ydec,Fsdn);


% Matlab esegui la decimazione di un fattore intero r automaticamente con la funzione
ydecM= decimate(y,r);	% qui usa un filtro IIR di Chebychev tipo I, di ordine 8
% NOTA: usa un "filtraggio doppio a fase zero"  (vedi funzione "filtfilt" di Matlab)
% ydecMfir= decimate(y,r,'fir');	% cos? usa invece un filtro FIR di ordine 30
figure(11); clf;
stem(tdnvec,ydecM,'m');
xlabel('t [s]'); ylabel('y_{decM}'); title('Segnale decimato da Matlab');
axis([t1 t2 -vx vx]);
% NOTA non c'? ritardo temporale poich? Matlab esegue un "Filtraggio doppio a fase nulla"
figure(12); clf;
plot(fdnvec,fftshift(abs(fft(ydecM))),'m');
xlabel('f [Hz]'); ylabel('Y_{decM}'); title('Spettro di ampiezza Segnale decimato da Matlab');
axis([-Fs/2 Fs/2 0 VX]);
% NOTARE che, impiegando un filtro peggiore, la funzione decimate
% usa una frequenza di taglio minore e produce uno spettro con valori nulli vicino a Fs/2/r
disp('premi un tasto per ascoltare il segnale decimato da Matlab'); pause;
soundsc(ydecM,Fsdn);