% Computer HW #3
% T Kite May 1998


% Set up the problem

fs=44100;					% sampling frequency
ts=1/fs;					% sampling period
al=2/ts;					% bilinear alpha...
al2=al*al;					% ...and alpha squared
c=1e-6;						% choose a capacitor
q=2;						% filter Q


% Compute the numerators and denominators of the digital filters

cf=[50 200 800 3200 12800]*2*pi;		% center frequencies
w0=2*fs*tan(cf/(2*fs))';			% prewarped
l=1./(c*w0.^2);					% inductors
r=sqrt(l/c)/q;					% resistors
beta=1+al*c*r+al2*l*c;
b=al*c*r./beta*[1 0 -1];
a=[ones(5,1) 2*(1-al2*l*c)./beta (1-al*c*r+al2*l*c)./beta];


% Plot the responses of the five filters

figure(1); clf
for t=1:5
  [h,w]=freqz(b(t,:),a(t,:),8192,fs);
  semilogx(w,abs(h))
  axis([20 20000 0 1])
  hold on
end
grid
xlabel 'Frequency (Hz)'
ylabel 'Magnitude'


% Plot the system response with the given slider positions

levdB=[-3 12 6 -12 6];				% gain in each band
K=10.^(levdB/20)-1;				% compute gains

ti=1:16384;					% generate chirp
x=sin(ti.*ti/16384*pi/2);

y=x;						% straight through output
for t=1:5					% add each filter output
  y=y+filter(b(t,:),a(t,:),x)*K(t);
end
h=fft(y)./fft(x);				% compute transfer function

figure(2)
f=linspace(0,fs/2,8193);
semilogx(f,20*log10(abs(h(1:8193))))
axis([20 20000 -12 12])
set(gca,'ytick',[-12 -9 -6 -3 0 3 6 9 12])
grid
xlabel 'Frequency'
ylabel 'Gain (dB)'
