%% Lezione del 2 dicembre 2009
%
% Si progetti un filtro risonatore discreto IIR del secondo ordine
%                   1
% H(z) = G -------------------------
%          1 + a1*z^(-1) + a2*z^(-2)
%
% (sez. 2.2.2 del libro di testo) che osserva le seguenti caratteristiche:
%
% r = 10^(-log10(dB/20)/(tau*Fs));
%
% a1 = -2*r*cos(2*pi*freq0/Fs);
% a2 = r^2;
%
% G = (1-r)*sqrt(1+a2-2*r*cos(4*pi*freq0/Fs));
%
% 1) si osservi 1 s della risposta all'impulso e il modulo della risposta in frequenza 
% campionando a Fs = 44100 Hz
%
% 2) si osservi la risposta a un rumore gaussiano di durata 1 s e il
% corrispondente modulo della risposta in frequenza 
%
% 3) si ascoltino le risposte calcolate ai punti 1) e 2)
%
% 4) verificata la distanza dall'origine dei poli complessi coniugati del
% risonatore, si collochino due zeri aventi lo stesso angolo dei due poli,
% rispettivamente a distanze uguali a 0.9990 e poi 1.0001 dall'origine: si
% osservi il modulo della risposta in frequenza nei due casi
%%

clear all

%%%%%%%%%% Simulation parameters

Fs = 44100;           % sampling frequency

n = [0:Fs-1];      % discrete-time domain (1s)

D = [0:Fs-1];          % frequency domain (Fs Hz)

%%%%%%%%%% Filter specs

freq0 = 2000;        % Hz
tau = .1;             % s
dB = 40;            % attenuation at tau

%%%%%%%%%% Input

x = [1,zeros(1,Fs-1)];   % impulse response
v = randn(1,Fs);         % 1s gaussian noise

%%%%%%%%%% Filter parameters

r = 10^(-log10(dB/20)/(tau*Fs));

a1 = -2*r*cos(2*pi*freq0/Fs);
a2 = r^2;

G = (1-r)*sqrt(1+a2-2*r*cos(4*pi*freq0/Fs));

%%%%%%%%%% Signal processing

y = zeros(1,Fs);       % malloc

y(1) = G*x(1);
y(2) = G*x(2) - a1*y(1);

for n = 3:Fs
    y(n) = G*x(n) - a1*y(n-1) - a2*y(n-2);
end

w = zeros(1,Fs);       % malloc

w(1) = G*v(1);
w(2) = G*v(2) - a1*w(1);

for n = 3:Fs
    w(n) = G*v(n) - a1*w(n-1) - a2*w(n-2);
end

%%%%%%%%%% Figure 1

figure(1)
plot(0:1/Fs:(Fs-1)/Fs,y);
xlabel('time (samples)'), ylabel('output'), title('Time response');

%%%%%%%%%% Figure 2

figure(2)
plot(0:1/Fs:(Fs-1)/Fs,w);
xlabel('time (samples)'), ylabel('output'), title('Time response');

%%%%%%%%%% Figure 3

Y = fft(y,Fs);    % FFT of y
W = fft(w,Fs);    % FFT of w

figure(3)
semilogx(0:Fs/2-1, 20*log10(abs(Y(1:Fs/2))), 0:Fs/2-1, 20*log10(abs(W(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency response');

%%%%%%%%%% Write wav files

wavwrite(y/max(y)-.01,Fs,'y.wav');
wavwrite(w/max(w)-.01,Fs,'w.wav');

%%%%%%%%%% Zeros

r_z = 0.9990;   % smaller than the two poles

b = conv([1 -r_z*exp(i*2*pi*freq0/Fs)],[1 -r_z*exp(-i*2*pi*freq0/Fs)]);  % numerator

y(1) = b(1)*x(1);
y(2) = b(1)*x(2) + b(2)*x(1) - a1*y(1);

for n = 3:Fs
    y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) - a1*y(n-1) - a2*y(n-2);
end

Y = fft(y,Fs);    % FFT of y

figure(4)
semilogx(0:Fs/2-1, 20*log10(abs(Y(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency response');

%%%%%%%%%% Zeros

r_z = 1.0001;   % larger than the two poles

b = conv([1 -r_z*exp(i*2*pi*freq0/Fs)],[1 -r_z*exp(-i*2*pi*freq0/Fs)]);  % numerator

y(1) = b(1)*x(1);
y(2) = b(1)*x(2) + b(2)*x(1) - a1*y(1);

for n = 3:Fs
    y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) - a1*y(n-1) - a2*y(n-2);
end

Y = fft(y,Fs);    % FFT of y

figure(5)
semilogx(0:Fs/2-1, 20*log10(abs(Y(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency response');
