% Esercitazione 2: estrazione di una sinusoide dal rumore AWGN
% tramite filtro IIR passabanda sintetizzato col metodo
% del piazzamento di poli e zeri

clear all;

% parametri generali
Fs= 16000;  % [Hz] freq. di campionamento



% SEGNALE SINUSOIDALE

% parametri della sinusoide
f0= 3500;   % [Hz] freq. sinusoide
tfin= 2;    % [s] durata sinusoide
t= (0:1/Fs:tfin);
A=0.2;      % ampiezza
x= A*cos(2*pi*f0*t);

Ns= floor(length(x)/2);	% numero campione da cui iniziare la visualizzazione
Nc= 200;	% numero campioni da visualizzare per grafici in t
% grafico segnale pulito (primi Nc campioni)
figure(1); clf;
stem(t(Ns:Ns+Nc),x(Ns:Ns+Nc));    % grafico segnale pulito
title('Sinusoide pulita');
xlabel('t [s]');
ylabel('x(t)');
axis([t(Ns) t(Ns+Nc) -1 1]);

% calcolo e visualizzazione spettro sinusoide pulita
fvec= linspace(-Fs/2,Fs/2,length(x));
X= fftshift(fft(x));  % spettro bilatero con freq. negative a sinistra (fftshift)
figure(2); clf;
plot(fvec,(1/length(X))*abs(X),'r');  % stampa spettro d'ampiezza
title('Spettro sinusoide pulita');
xlabel('f [Hz]');
ylabel('|X(f)|');

disp('schiaccia un tasto per ascoltare la sinusoide pulita')
pause;
sound(x,Fs);    % ascolto




% AGGIUNTA DEL RUMORE

Pn= 1.0;	% valore quadratico medio (= potenza media = varianza)
% campioni di rumore Gaussiano indipendenti a media nulla
n= sqrt(Pn)*randn(1,length(x));  % vettore riga campioni AWGN

y= x+n; % segnale rumoroso

% grafico segnale rumoroso e suo spettro
figure(3); clf;
stem(t(Ns:Ns+Nc),y(Ns:Ns+Nc));    % grafico segnale rumoroso
title('Sinusoide rumorosa (il segnale sin. è mascherato dal rumore)');
xlabel('t [s]');
ylabel('y(t)=x(t)+n(t)');
axis([t(Ns) t(Ns+Nc) -1 1]);

Y= fftshift(fft(y));  % spettro
figure(4); clf;
plot(fvec,(1/length(X))*abs(Y),'r');  % stampa spettro d'ampiezza
title('Spettro sinusoide rumorosa (il sin. è visibile nello spettro)');
xlabel('f [Hz]');
ylabel('|Y(f)|');

disp('schiaccia un tasto per ascoltare la sinusoide rumorosa')
pause;
sound(y,Fs);    % ascolto



% FILTRAGGIO PER ELIMINARE IL RUMORE

% sintesi filtro passabanda con piazzamento poli/zeri
DelF= 10;  % [Hz] banda a 3dB del filtro (parametro di progetto
zer= [1; -1];    % zeri in f/Fs= 0 e +-Fs/2
r= 1-pi*DelF/Fs;
if (r<=0.9), disp('ATTENZIONE: banda a 3dB troppo larga'), end;
pol= [r*exp(j*2*pi*f0/Fs);  r*exp(-j*2*pi*f0/Fs)];  % poli coniugati a freq. +-f0

% controllo della posizione di poli e zeri nel piano complesso z
figure(5); clf;
zplane(zer,pol);

% sintesi dei polinomi a numeratore e denominatore del filtro IIR
% occorre speficicare il "guadagno" k come costante moltiplicativa di H(z) razionale
k= 2e-3;    % guadagno asintotico (z-> inf): è un parametro libero
[b,a]= zp2tf(zer,pol,k);

% visualizzazione risposta in ampiezza, risposta in fase e ritardo di gruppo
figure(6); clf;
freqz(b,a,128,Fs);
% la risposta in ampiezza si può traslare (in dB) variando k
% quella in fase rimane così com'è: non si può controllare
figure(7); clf;
grpdelay(b,a,128,Fs);   % ritardo di gruppo= -d(fase)/dw

disp('schiaccia un tasto per applicare il filtro al segnale rumoroso')
pause;

% Generazionesegnale filtrato
yout= filter(b,a,y);

% visualizzazione segnale filtrato e suo spettro
figure(8); clf;
stem(t(Ns:Ns+Nc),yout(Ns:Ns+Nc));    % grafico segnale filtrato
title('Sinusoide rumorosa filtrata');
xlabel('t [s]');
ylabel('y_{out}(t)');
axis([t(Ns) t(Ns+Nc) -1 1]);

Yout= fftshift(fft(yout));  % spettro bilatero con freq. negative a sinistra (fftshift)
figure(9); clf;
plot(fvec,(1/length(Yout))*abs(Yout),'r');  % stampa spettro d'ampiezza segnale filtrato
title('Spettro sinusoide rumorosa');
xlabel('f [Hz]');
ylabel('|Y_{out}(f)|= |Y(f)| |H(f)|');

disp('schiaccia un tasto per ascoltare la sinusoide rumorosa filtrata')
pause;
sound(yout,Fs);    % ascolto