% Example bandpass filter design
% Course: Analog and Digital Filter Design
% Lecture: Op amp RC filter realizations
% Authors: Miroslav D. Lutovac, Dejan V. Tosic
% Email: lutovac@kondor.etf.bg.ac.yu
% Web: http://kondor.etf.bg.ac.yu/~lutovac/
% Copyright (c) 2005 by Lutovac & Tosic
% $Revision: 1.01 $  $Date: 2005/02/24$
% See also: Miroslav D. Lutovac, Dejan V. Tosic, Brian L. Evans
%  Filter Design for Signal Processing Using MATLAB and Mathematica
%  Prentice Hall (c)2001. ISBN 0-201-36130-2

clear all, close all, clc

% Add DrawFilt to the MATLAB path
path(path,'C:\afd\DrawFilt')

% Give specification
Fs1 = 1200;
Fp1 = 2500;
Fp2 = 3000;
Fs2 = 5500;
Ap = 1;
As = 20;
Fp  = sqrt(Fp1*Fp2);

% Choose approximation: Chebyshev Type I filter
[n, Wn] = cheb1ord(2*pi*[Fp1 Fp2], 2*pi*[Fs1 Fs2], Ap, As,'s');

% Compute transfer function
[z,p,K] = cheby1(n, Ap, Wn, 's');

% Decompose transfer function into biquads
p = sort(p);
Q = -abs(p)./(2*real(p));
w = abs(p);
Q = str2num(num2str(Q,3));
D1 = [1 w(1)/Q(1) w(1)^2];
D2 = [1 w(3)/Q(3) w(3)^2];
N1 = [0 w(1)/Q(1) 0];  
N2 = [0 w(3)/Q(3) 0];
k1 = abs(polyval(N1,i*2*pi*Fp)/polyval(D1,i*2*pi*Fp));
k2 = abs(polyval(N2,i*2*pi*Fp)/polyval(D2,i*2*pi*Fp));
N1 = N1/k1; N2 = N2/k2;

% Plot attenuation
Fplot1 = 1*1000;  Fplot2 = 6*1000;
f = Fplot1:(Fplot2-Fplot1)/200:Fplot2;
H1 = freqs(N1,D1,2*pi*f);
A1 = -20*log10(abs(H1));
H2 = freqs(N2,D2,2*pi*f);
A2 = -20*log10(abs(H2));

plot(f,A1+A2,'g-',f,A1,'r--',f,A2,'b:')
xlabel('Frequency (Hz)')
ylabel('Attenuation (dB)')
legend('A_1+A_2','A_1','A_2')
text(2000, 42, ['\omega_1 = ' num2str(round(w(1))) ',  Q_1 = ' num2str(Q(1))],'color','r')
text(2000, 32, ['\omega_2 = ' num2str(round(w(3))) ',  Q_2 = ' num2str(Q(3))],'color','b')

% Choose realization from album
figure 
drawbpmq0(0,4,4,5,8,'r','\bf Biquad #1')
drawbpmq0(16,3,4,5,8,'b','\bf Biquad #2')
drawfout(34,12,4,5,8,'m')
drawfin(2,14,4,5,8,'m')

disp('======== Biquad 1 ========')
kilo = 10^3;
mega = 10^6;
nano = 10^(-9);

% Given a second-order lowpass transfer function identify
wp = w(1);
Qp = Q(1);

% Choose K
K  = 1/k1;

% Choose P
P  = 1; 

% Choose C2, C3, R6
C2 = 3.3*nano;
C3 = 3.3*nano;
R6 = 10*kilo;

disp(['wp = ' num2str(wp)])
disp(['Qp = ' num2str(Qp)])
disp(['K = ' num2str(K)])
disp(['C2 = ' num2str(C2/nano) ' nF'])
disp(['C3 = ' num2str(C3/nano) ' nF'])
disp(['P = ' num2str(P)])
disp('-------------------------')

% Compute element values
R1 = 1/(wp*sqrt(C2*C3*P));
R4 = P*R1;
R5 = R6*((1+C2/C3)/P-sqrt(C2/(P*C3))/Qp);
Ko = Qp*(1+R5/R6)*sqrt(P*C3/C2);
R11 = R1*Ko/K;
R12 = R1*Ko/(Ko-K); 

if abs(R12)>10^(7); R12 = 10^7; end

% Round element values to 3 digits of precision
C2 = str2num(num2str(C2,3));   C3  = str2num(num2str(C3,3));
R5  = str2num(num2str(R5,3));  R6  = str2num(num2str(R6,3));
R4  = str2num(num2str(R4,3));
R11 = str2num(num2str(R11,3)); R12 = str2num(num2str(R12,3));

disp(['C2 = ' num2str(C2/nano) ' nF'])
disp(['C3 = ' num2str(C3/nano) ' nF'])
disp(['R11 = ' num2str(R11/kilo) ' kOhm'])
disp(['R12 = ' num2str(R12/mega) ' MOhm'])
disp(['R4 = ' num2str(R4/kilo) ' kOhm'])
disp(['R5 = ' num2str(R5/kilo) ' kOhm'])
disp(['R6 = ' num2str(R6/kilo) ' kOhm'])

% Validate realization
% Generate transfer function in terms of element values
num = [0 -C3*R12*R4*(R5+R6) 0];
den = [(C2*C3*R11*R12*R4*R6) ...
       -(C3*R11*R4*R5+C3*R12*R4*R5-C2*R11*R12*R6-C3*R11*R12*R6) ...
       (R11*R6+R12*R6)];
H1r = freqs(num,den,2*pi*f);

figure
subplot(2,1,1)
drawbpmq0(0,0,4,5,9,'b','Biquad #1, BP-MQ                ')
drawfout(18,9,4,5,9,'m')
drawfin(2,10,4,5,9,'m')

subplot(2,1,2)
plot(f,abs(H1r),f(1:10:end),abs(H1(1:10:end)),'x')
text(2700, 1.5, ...
    ['R_{5} = ' num2str(R5/kilo) ' k\Omega, ' ...
     'R_{6} = ' num2str(R6/kilo) ' k\Omega'])
text(3000, 1, ...
    ['R_{11} = ' num2str(R11/kilo) ' k\Omega, ' ...
     'R_{12} = ' num2str(R12/kilo) ' k\Omega, ' ...
     'R_{4} = ' num2str(R4/kilo) ' k\Omega'])
text(3300, 0.5, ...
    ['C_{2} = ' num2str(C2/nano) ' nF, ' ...
     'C_{3} = ' num2str(C3/nano) ' nF'])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
legend('realized','ideal')


disp('======== Biquad 2 ========')
% Given a second-order lowpass transfer function identify
wp = w(3);
Qp = Q(3);

% Choose K
K  = 1/k2;

% Choose P
P  = 1; 

% Choose C2, C3, R6
C2 = 3.3*nano;
C3 = 3.3*nano;
R6 = 10*kilo;

y = inline('-(sqrt(C2)*C3*(1+K)*abs(P))+C2*sqrt(C3*abs(P))*Qp+C3*sqrt(C3*abs(P))*(1+abs(P))*Qp',...
           'P','C2','C3','Qp','K');

disp(['wp = ' num2str(wp)])
disp(['Qp = ' num2str(Qp)])
disp(['K = ' num2str(K)])
disp(['C2 = ' num2str(C2/nano) ' nF'])
disp(['C3 = ' num2str(C3/nano) ' nF'])
disp(['P = ' num2str(P)])
disp('-------------------------')

% Compute element values
R1 = 1/(wp*sqrt(C2*C3*P));
R4 = P*R1;
R5 = R6*((1+C2/C3)/P-sqrt(C2/(P*C3))/Qp);
Ko = Qp*(1+R5/R6)*sqrt(P*C3/C2);
R11 = R1*Ko/K;
R12 = R1*Ko/(Ko-K); if abs(R12)>10^(7); R12 = 10^7; end

% Round element values to 3 digits of precision
C2 = str2num(num2str(C2,3));   C3  = str2num(num2str(C3,3));
R5  = str2num(num2str(R5,3));  R6  = str2num(num2str(R6,3));
R4  = str2num(num2str(R4,3));
R11 = str2num(num2str(R11,3)); R12 = str2num(num2str(R12,3));

disp(['C2 = ' num2str(C2/nano) ' nF'])
disp(['C3 = ' num2str(C3/nano) ' nF'])
disp(['R11 = ' num2str(R11/kilo) ' kOhm'])
disp(['R12 = ' num2str(R12/mega) ' MOhm'])
disp(['R4 = ' num2str(R4/kilo) ' kOhm'])
disp(['R5 = ' num2str(R5/kilo) ' kOhm'])
disp(['R6 = ' num2str(R6/kilo) ' kOhm'])

% Validate realization
% Generate transfer function in terms of element values
num = [0 -C3*R12*R4*(R5+R6) 0];
den = [(C2*C3*R11*R12*R4*R6) ...
       -(C3*R11*R4*R5+C3*R12*R4*R5-C2*R11*R12*R6-C3*R11*R12*R6) ...
       (R11*R6+R12*R6)];
H2r = freqs(num,den,2*pi*f);

figure
subplot(2,1,1)
drawbpmq0(0,0,4,5,9,'b',' Biquad #2, BP-MQ                ')
drawfout(18,9,4,5,9,'m')
drawfin(2,10,4,5,9,'m')

subplot(2,1,2)
plot(f,abs(H2r),f(1:10:end),abs(H2(1:10:end)),'x')
text(1100, 1.4, ...
    ['R_{5} = ' num2str(R5/kilo) ' k\Omega, ' ...
     'R_{6} = ' num2str(R6/kilo) ' k\Omega'])
text(3900, 1, ...
    ['R_{11} = ' num2str(R11/kilo) ' k\Omega, ' ...
     'R_{12} = ' num2str(R12/kilo) ' k\Omega'])
text(1100, 0.9, ...
    ['R_{4} = ' num2str(R4/kilo) ' k\Omega'])
text(3900, 0.5, ...
    ['C_{2} = ' num2str(C2/nano) ' nF, ' ...
     'C_{3} = ' num2str(C3/nano) ' nF'])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
legend('realized','ideal')

% Validate realization overall response
figure
A1r = -20*log10(abs(H1r));
A2r = -20*log10(abs(H2r));
plot(f,A1r+A2r,'b-',...
     f(1:1:end),A1(1:1:end)+A2(1:1:end),'mx',...
     f,A1+A2-(A1r+A2r),'r-')
axis([Fp1 Fp2 -1.2*Ap 1.2*Ap])
legend('realized','ideal','error')
xlabel('Frequency (Hz)')
ylabel('Attenuation (dB)')

figure
plot(f,A1r+A2r,'b-',...
     f(1:5:end),A1(1:5:end)+A2(1:5:end),'mx',...
     f,A1+A2-(A1r+A2r),'r-')
axis([Fs1 Fs2 -1.2*Ap 1.2*As])
legend('realized','ideal','error')
xlabel('Frequency (Hz)')
ylabel('Attenuation (dB)')
