%% M4005.m, FIR filter design

%% Read in demo signal
[x, fT] = wavread('M4005.wav');
soundsc(x, fT);

figure(1); clf;
spectrogram(x, 512, 256, 512, fT, 'yaxis');
colorbar;

%% Filter specifications
fT                             % from the signal!
F   = [     4500 5000   5800 6300      ];
A   = [  1            0            1   ];
DEV = [ 0.1          0.01         0.1  ];
figure(2); clf;
speksitFIR(F, A, DEV, fT);
xlabel('Frequency (Hz)'); ylabel('|H(e^{j \omega})|');
title('FIR bandstop filter');
grid on;

%% Filter order 
[ fill here ] = firpmord( fill here );

%% Filter coefficients
% Transfer function H(z) = B(z)/A(z)
B   = firpm( fill here );
A   = 1;                   % FIR: no feedback <=> A(z) == 1

%% Compute and plot the amplitude response
[H, wF] = freqz(B, A, 512, fT);
figure(2); hold on;
plot(wF, abs(H));
grid on;

%% Plot the phase response
figure(3); clf;
plot(wF, unwrap(angle(H)));
grid on;
%phasez(B, A, 512, fT);
title('Phase response \angle H(e^{j \omega})');

%% Group-delay of the filter
% -(d/dw) <H(e^jw)
figure(4); clf;
grpdelay(B, A, 512, fT);
grid on;
title('Group delay \tau(e^{j \omega})');

%% Pole-zero plot
figure(5); clf;
zplane(B, A)

%% Filter demo signal
y  = filter(B, A, x);

%% Output spectrogram
figure(6); clf;
spectrogram(y, 512, 256, 512, fT, 'yaxis');
colorbar;

%% Listen to filtered signal
soundsc(y, fT);

%% Magnitude response in decibels
figure(7); clf;
magdB = 20*log10(abs(H)/max(abs(H)));
plot(wF, magdB);
grid on;
xlabel('Frequency (Hz)');
title('Magnitude response in decibels');
ylabel('dB');
 