%% Matlab #4 / Problem 2 [M4002]
%
% Digital IIR filter design
%
% (i)   Determine the specifications of the filter.
% (ii)  Choose the type of the filter.
% (iii) Estimate the order of the filter.
% (iv)  Compute the coefficients of transfer function.
% (v)   Check that the specifications are fulfilled!

%% Read and analyze the signal
[x, fT]  = wavread('M4002.wav');
t = [0 : length(x)-1]/fT;
figure(1); clf;
plot(x);		   % NOTE! Index values at x-axis, not seconds
plot(t,x);                 % NOTE! Seconds at x-axis, not index values

figure(2); clf;              
spectrogram(x, 512, 256, 512, fT, 'yaxis');
title('Input x[n]');

%% (i) Determine specifications
fT                           % sampling frequency of signal & system
Wp=0.4; Ws=0.5; Rp=1.5; Rs=18; %  * * * * * DELETE THIS LINE AWAY! * * * * 
%fp = ;                       % fill in correct value!
%fs = ;                       % fill in correct value! 
%Wp = 2*fp/fT;                % scaled frequency 0 < Wp < 1
%Ws = 2*fs/fT;                % scaled frequency 0 < Ws < 1
%Rp = ;                       % fill in correct value!
%Rs = ;			     % fill in correct value!
figure(3); clf;
speksit(Wp, Ws, Rp, Rs, '', fT); % draw the specs!
grid on; xlabel('Hz'); ylabel('dB'); title('|H(e^{j \omega})|^2');

%% (ii), (iii) Estimate filter order, help ellipord
%[N, Wn] = ellipord(  );      % add correct variables

%% (iv) Compute filter coefficients, help ellip
B=1; A=1;                    %  * * * * * DELETE THIS LINE AWAY! * * * *
%[B, A] = ellip(  );	      % add correct variables

filterOrder = max([length(B) length(A)])-1 % order of the filter

tf2latex(B, A)               % tf2latex.m from course web site

% pole-zero plot:
figure(4); clf;
zplane(B, A);                % pole-zero plot

%% (v) Check if specifications are fulfilled
[H, w] = freqz(B, A, 512, fT);
magdB  = 20*log10(abs(H)/max(abs(H))); % squared magnitude in dB
figure(3); hold on;          % draw into the same figure No. 1
plot(w, magdB);              % amplitude response in decibels

%% (vi) Use the filter and test it
y  = filter(B, A, x);        % apply X -> H(z) -> Y

%% Analyze the input and output signal
figure(1); hold on;
py = plot(y);
set(py, 'Color', [1 0 0], 'LineStyle', '-.', 'LineWidth', 2);

% Draw spectrogram
figure(5), clf;
spectrogram(y, 512, 256, 512, fT, 'yaxis');
title('Output y[n]');

xF = 20*log10(abs(fft(x))); 
yF = 20*log10(abs(fft(y)));
figure(6), clf;
plot(xF); hold on;           % NOTE! Index values at x-axis, not Hertz
plot(yF,'r'); title('Spectra');  % See Matlab 2 / Problem 3 for Hertz

%% ... listen to the signals
soundsc(x, fT); pause;       % original; press any key to continue
%
soundsc(y, fT);              % filtered
