%----------------------------------------------%
%-- 1) Evaluating a 2nd-order feedback filter--%
%----------------------------------------------%
 
%---------- b: Plotting the poles and zeros on the z-plane -------------%
 
b = [1 -1.414 1];                  % coefficients of the numerator           
a = [1 -1.273 0.81];               % coefficients of the denominator
B = roots(b);                      % zeros of the transfer function
A = roots(a);                      % poles of the transfer function
figure();
zplane(B,A);                       % plotting the pole-zero plot
 
%------- c: Converting the zeros and poles to polar coordinates --------%
 
zeros_mag = abs(B);                % magnitude of the zeros
poles_mag = abs(A);                % magnitude of the poles
zeros_ang = angle(B)*(180/pi);     % angle of the zeros (degrees)
poles_ang = angle(A)*(180/pi);     % angle of the poles (degrees)
 
% To represent the poles and zeros in polar coordinates we write the
% magnitude followed by the angle
 
%------- d: Evaluating this transfer function --------%
 
fs = 48000;  % note this is NOT 44100! 
freq = linspace(0, fs/2, 1025);
freq = freq(1:(end-1)); 
w = (freq*2*pi)/fs; 
z = exp(1i*w);
Hz = (1-1.414*z.^(-1)+z.^(-2))./(1-1.273*z.^(-1)+0.81*z.^(-2));  % the transfer function
 
Hz_mag = abs(Hz);                  % the magnitude of Hz
Hz_phase = angle(Hz);              % the angle of Hz
 
figure();
subplot(2,1,1);
plot(freq,Hz_mag);
title('Magnitude Response of H(z)');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
 
subplot(2,1,2);
plot(freq,Hz_phase);
title('Phase Response of H(z)');
xlabel('Frequency (Hz)');
ylabel('Phase Angle (Rad)');
 
%------- e: Plotting the frequency response --------%
 
[H, W] = freqz(b,a,1024);       % H is the frequency response vector
W = W*fs/(2*pi);                % W is the n-point frequency vector
                                % converted to hertz
figure();                       % Plotting the magnitude response using
subplot(2,1,1);                 % W (converted to hertz) for the x-axis 
plot(W,abs(H));                 % and the magnitude of the frequency response
title('Magnitude Response');    % (the absolute value of H) for the y-axis
xlabel('Frequency (Hz)');
ylabel('Magnitude');
 
subplot(2,1,2);                 % Plotting the phase response using W (hz) 
plot(W,angle(H));               % for the x-axis and the phase response of
title('Phase Response');        % the frequency response for the y-axis
xlabel('Frequency (Hz)');
ylabel('Phase Angle (Rad)');
 
%----------------------------------------------%
%-- 2) Designing a 2nd-order resonant filter ---%
%----------------------------------------------%
 
% Specs:
%   fc = 2000Hz
%   BW = 1000Hz
 
%---------- a: A 2nd-order resonant bandpass filter -------------%
 
b = 0.040;                      % coefficients of the numerator
a = [1 -1.778 0.863];           % coefficients of the denominator
B = roots(b);                   % the zeros
A = roots(a);                   % the poles
figure()
zplane(B,A);                    % plotting the pole-zero plot
 
[H, W] = freqz(b,a,1024);
W = W*fs/(2*pi);                % the vector of frequencies must be in degrees
 
figure();                       % plotting the magnitude response
subplot(2,1,1);
plot(W,abs(H));
axis([0 5000 0 1.5]);
title('Magnitude Response');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
 
subplot(2,1,2);                 % plotting the phase response
plot(W,angle(H));
axis([0 5000 -2.5 2.5]);
title('Phase Response');
xlabel('Frequency (Hz)');
ylabel('Phase Angle (Rad)');
 
 
%---------- b: Calculating the filter's poles and zeros -------------%
 
zeros_mag = abs(B);                 % magnitude of the zeros
poles_mag = abs(A);                 % magnitude of the poles
zeros_angle = angle(B)*(180/pi);    % angle of the zeros (degrees)
poles_angle = angle(A)*(180/pi);    % angle of the poles (degrees)
 
% To represent the poles and zeros in polar coordinates we write the
% magnitude followed by the angle
 
 
%---------- c: Using the filter to filter a wav file -------------%
 
[y, fs, nbits] = wavread('Alexei_20SULTANOV_20_20Chopin_20Sonata_20_3_20Fina',[10400000 12150000]);
yfiltered = filter(b,a,y);
 
sound(y,fs);
sound(yfiltered,fs);
 
wavwrite(y,fs,'NatePaternoster-Lab5-2original.wav');
wavwrite(yfiltered,fs,'NatePaternoster-Lab5-2filtered.wav');