% M09-bpfilt.diary
% Simple bandpass filter example
% 2007-10-18

% w_c = 0.4 pi, bandwidth B = 0.1 pi
wc = 0.4*pi;
B = 0.1*pi;
% hence...
beta = cos(wc)
%beta =
%    0.3090
% quadratic in alpha for B
roots([cos(B) -2 cos(B)])
%ans =
%    1.3764
%    0.7265
% |alpha| must be < 1 for stable system
alpha = ans(2)
%alpha =
%    0.7265
% Hence numerator and denominator coefficients
b = (1-alpha)/2*[1 0 -1];
a = [1, -beta*(1+alpha), alpha];
% Plot the response
freqz(b,a)
% Adjust the axes for more detail in passband
axis([0 1 -20 0])
% Peak at 0.4pi, -3 dB bandwidth about 0.1pi

% Apply the filter to a signal
[x,sr]=wavread('mpgr1_sx419.wav');
soundsc(x,sr);
% Apply the filter
y = filter(b,a,x);
soundsc(y,sr);
% Bandpassed sound!  

% Visualize it: take spectrum of first 1 sec of sound
subplot(311)
X = fft(x,16384);
plot([0:8192]/8192, 20*log10(abs(X(1:8193))));  % plot in dB, 1st half
% Make all the dB axes the same; 80 dB (-60..20) is plenty of range
axis([0 1 -60 20])
grid
% compare to spectrum of output
subplot(312)
Y = fft(y,16384);
plot([0:8192]/8192, 20*log10(abs(Y(1:8193))));
axis([0 1 -60 20])
grid
% .. and compare to spectrum of filter
[H,W] = freqz(b,a);
subplot(313)
plot(W/pi, 20*log10(abs(H)));
axis([0 1 -60 20])
grid
% Output is clearly input spectrum 'shifted' by magnitude response
% (multiplication becomes an additive shift with a db (log) axis)

% .. or try spectrograms to see whole sound
subplot(211)
specgram(x,512,sr);
subplot(212)
specgram(y,512,sr);