%% Matlab #4 / Problem 3 [M4005]
%
% Digital FIR filter design
% using equiripple FIR (Parks-McClellan)

%% Specifications
fT  = 22050;
F   = [    3300  4000      ];
A   = [0                  1];
DEV = [0.1             0.01];
figure(1); clf;
speksitFIR(F, A, DEV, fT);
grid on; 
xlabel('frequency (Hz)'); 
ylabel('|H(e^{j \omega})|');
title('Specifications for equiripple FIR filter: \delta_s = 0.1, \delta_p = 0.01');

%% Filter order 
%[  , ] = firpmord( , , ,  );  % add correct variables

%% Filter coefficients
% Transfer function H(z) = B(z)/A(z)
b   = 1;                   % * * * * * DELETE THIS LINE AWAY! * * * * 
%b   = firpm( fill here );  % add correct variables
a   = 1;                   % FIR: no feedback <=> A(z) == 1

filterOrder = max([length(b) length(a)])-1 % order of the filter

tf2latex(b, a)			% tf2latex.m from course web site

%% Compute and plot the amplitude response
[H, w] = freqz(b, a, 512, fT);
figure(1); hold on;
plot(w, abs(H));
grid on;
title('Specifications and magnitude response |H(e^{j\omega})|');

%% Plot the phase response
figure(2); clf;
plot(w, unwrap(angle(H)));
grid on;
xlabel('frequency (Hz)'); 
ylabel('\angle H(e^{j \omega})');
title('Phase response');

%% Group-delay of the filter
% -(d/dw) <H(e^jw)
figure(3); clf;
grpdelay(b, a, 512, fT);
xlabel('frequency (Hz)'); 
ylabel('\tau(e^{j \omega})');
title('Group delay');

%% Pole-zero plot
figure(4); clf;
zplane(b, a)

%% Filter kiisseli.wav
[x, dummy] = wavread('kiisseli.wav');
y  = filter(b, a, x);

%% Listen
soundsc(x, fT);
pause;                               % press any key to continue
soundsc(y, fT);

%% Spectrograms
figure(5); 
spectrogram(x, 512, 256, 512, fT, 'yaxis');
title('Input x[n]');

figure(6); 
spectrogram(y, 512, 256, 512, fT, 'yaxis');
title('Output y[n]');

%% Magnitude response in logarithmic scale
figure(7); 
plot(w, 20*log10(abs(H)/max(abs(H))));  % with scaling max to unity
grid on;
xlabel('frequency (Hz)'); ylabel('dB, |H(e^{j \omega})|^2');
title('Squared magnitude response in log. scale');

