%   Esercizio assegnato relativamente alla
%   reiezione di un disturbo sinusoidale nella amplificazione
%   audio dovuto all' effetto Larsen.
%   Soluzione sviluppata da Francesco Zanini.
%   22 / 10 / 2004.
clear all;
load('handel');      % Fs= 8192 Hz : freq. di campionamento
f0= 2000;   % [Hz] freq. sinusoide
tfin= 5;    % [s] durata sinusoide
t= (0:1/Fs:tfin);
A= 1.0;      % ampiezza
x= A*cos(2*pi*f0*t);
fvec= linspace(-Fs/2,Fs/2,length(x));

% parametri per la visualizzazione di segnali e spettri
t1= 2.5; t2= 2.8;% [s]
f1= -Fs/2; f2= Fs/2;	% [Hz]
vx= 2.0;	% ampiezza segnali

X= fftshift(fft(x));  % spettro disturbo sinusoidale
figure(1); clf;
subplot(1,2,1);
plot(t,x,'r');
axis([t1 t2 -vx vx]);
xlabel('t [s]'); ylabel('x(t)'); title('disturbo sinusoidale');
subplot(1,2,2);
plot(fvec,abs(X),'r');  % stampa spettro d'ampiezza
temp= axis; axis([f1 f2 0 temp(4)]);
xlabel('f [Hz]'); ylabel('!X(f)|'); title('spettro disturbo sinusoidale');
disp('premi un tasto per ascoltare'); pause;
soundsc(x,Fs);

y=y(1:length(x));
Y= fftshift(fft(y));  % spettro handel pulito
figure(2); clf;
subplot(1,2,1);
plot(t,y,'b');
axis([t1 t2 -vx vx]);
xlabel('t [s]'); ylabel('y(t)'); title('segnale pulito');
subplot(1,2,2);
plot(fvec,abs(Y),'b');  % stampa spettro d'ampiezza
temp= axis; axis([f1 f2 0 temp(4)]);
xlabel('f [Hz]'); ylabel('!Y(f)|'); title('spettro segnale pulito');
disp('premi un tasto per ascoltare'); pause;
soundsc(y,Fs);

yn=y+x';
Yn= fftshift(fft(yn));  % spettro handel con larsen
figure(3); clf;
subplot(1,2,1);
plot(t,yn,'k');
axis([t1 t2 -vx vx]);
xlabel('t [s]'); ylabel('y_n(t)'); title('segnale con Larsen');
subplot(1,2,2);
plot(fvec,abs(Yn),'k');  % stampa spettro d'ampiezza
temp= axis; axis([f1 f2 0 temp(4)]);
xlabel('f [Hz]'); ylabel('|Y_n(f)|'); title('spettro segnale con Larsen');
disp('premi un tasto per ascoltare'); pause;
soundsc(yn,Fs);

% sintesi filtro reiettabanda con piazzamento poli/zeri

% Gli zeri sopprimono la frequenza incriminata
zer= [exp(j*2*pi*f0/Fs);exp(-j*2*pi*f0/Fs)];    % zeri in +-f0/Fs
% i poli servono a "tirare su" la risposta in ampiezza
% nell'intorno della freq. f0 (provare con solo zeri per riscontro)
DelF= 50;  % [Hz] banda a 3dB del filtro
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
figure(4); clf;
zplane(zer,pol);  % poli-zeri del filtro sul piano z
title('poli-zeri del filtro sul piano z');

k= 0.981;    % in modo da avere guadagno H(z=0)=1: fissato per tentativi
[b,a]= zp2tf(zer,pol,k);
figure(5); clf;
freqz(b,a,128,Fs);   % risposta in frequenza del filtro
title('risposta in frequenza del filtro');

figure(6); clf;
grpdelay(b,a,128,Fs);   % ritardo di gruppo
title('ritardo di gruppo del filtro');
disp('premi un tasto per continuare'); pause;


% segnale filtrato ottenuto

yout= filter(b,a,yn);
Yout= fftshift(fft(yout));  % spettro
figure(7); clf;
subplot(1,2,1);
plot(t,yout,'m');
axis([t1 t2 -vx vx]);
xlabel('t [s]'); ylabel('y_{out}(t)'); title('segnale filtrato');
subplot(1,2,2);
plot(fvec,abs(Yout),'m');  % stampa spettro d'ampiezza segnale filtrato
temp= axis; axis([f1 f2 0 temp(4)]);
xlabel('f [Hz]'); ylabel('|Y_{out}(f)|'); title('spettro del segnale filtrato');
disp('premi un tasto per ascoltare'); pause;
soundsc(yout,Fs);    % ascolto

% segnale Errore
% NOTARE CHE L'ERRORE E' CONCENTRATO SULLE FREQUENZE PROSSIME AL NOTCH DEL FILTRO
Dy= yout-y;
DY= fftshift(fft(Dy));  % spettro
figure(8); clf;
subplot(1,2,1);
plot(t,Dy,'g');
axis([t1 t2 -vx vx]);
xlabel('t [s]'); ylabel('\Deltay(t)'); title('segnale errore y_{out}(t)-y(t)');
subplot(1,2,2);
plot(fvec,abs(DY),'g');  % stampa spettro d'ampiezza segnale errore
temp= axis; axis([f1 f2 0 temp(4)]);
xlabel('f [Hz]'); ylabel('|\Delta Y(f)|'); title('spettro del segnale errore');
disp('premi un tasto per ascoltare'); pause;
soundsc(Dy,Fs);    % ascolto