% Example lowpass 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
Fp = 3400;
Fs = 5000;
Ap = 1;
As = 20;

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

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

% Decompose transfer function into biquads
p = sort(p);
Q = -abs(p)./(2*real(p));
w = abs(p);
D1 = [1 w(1)/Q(1) w(1)^2];
D2 = [1 w(3)/Q(3) w(3)^2];
N1 = [0 0 w(1)^2];
N2 = [0 0 w(3)^2];

% Plot attenuation
Fplot = 6000;
f = 0:Fplot/200:Fplot;
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(500, 25, ['\omega_1 = ' num2str(round(w(1))) ',  Q_1 = ' num2str(Q(1))],'color','r')
text(500, 18, ['\omega_2 = ' num2str(round(w(3))) ',  Q_2 = ' num2str(Q(3))],'color','b')

% Choose realization from album
figure 
drawlplq0(0,2,4,5,8,'r','\bf Biquad #1')
drawlpmq0(11,0,4,5,8,'b','\bf Biquad #2')
drawfout(32,6,4,5,8,'m')
drawfin(2,8,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 0 < K < 1 
K  = 1-(1e-14);

% Choose C4
C4 = 3.3*nano;

% Optimize P to achieve C2 = C2opt (nF)
C2opt = 12;
Pinit = 3; % initial guess
y = inline('Qp*Qp*C4*(1+P)^2/P-C2opt','P','C4','Qp','C2opt');

% MATLAB 7 syntax
% P = fzero(y,Pinit,3,C4/nano,Qp,C2opt);

% MATLAB 6.5 syntax
[P, objectiveFunctionValue] = fzero(y,Pinit,[],C4/nano,Qp,C2opt);

if isnan(P); P=Pinit; end;

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

% Compute element values
C2 = 2*Qp^2*C4*(1+(1+P^2)/(2*P));
R1 = 1/(wp*C4*Qp*(1+P));
R3 = P*R1;
R11 = R1/K;
R12 = R1/(1-K);

% Large values of R12 are interpreted as open-circuited branches
if abs(R12)>10^7; R12 = 10^7; end
% Large R12 does not appear in implementation!

% Round element values to 3 digits of precision
R11 = str2num(num2str(R11,3)); 
R12 = str2num(num2str(R12,3));
R3  = str2num(num2str(R3,3));
C2  = str2num(num2str(C2,3)); 
C4  = str2num(num2str(C4,3));

disp(['R11 = ' num2str(R11/kilo) ' kOhm'])
disp(['R12 = ' num2str(R12/mega) ' MOhm'])
disp(['R3 = ' num2str(R3/kilo) ' kOhm'])
disp(['C2 = ' num2str(C2/nano) ' nF'])
disp(['C4 = ' num2str(C4/nano) ' nF'])

% Validate realization
% Generate transfer function in terms of element values
num1 = R12;
den1 = [(C2*C4*R11*R12*R3) (C4*R11*R12+C4*R11*R3+C4*R12*R3) (R11+R12)];
H1r = freqs(num1,den1,2*pi*f);

figure
subplot(2,1,1)
drawlplq0(0,0,4,5,10,'r','Biquad #1, LP-LQ')
drawfout(17,5,4,5,10,'m')
drawfin(2,6,4,5,10,'m')

subplot(2,1,2)
plot(f(1:10:end),abs(H1(1:10:end)),'mx',f,abs(H1r),'r')
legend('ideal','realized')
text(2300, 0.7, ...
    ['R_{11} = ' num2str(R11/kilo) ' k\Omega, ' ...
     'R_{12} = ' num2str(R12/mega) ' M\Omega, ' ...
     'R_{3} = ' num2str(R3/kilo) ' k\Omega'])
text(3500, 0.4, ...
    ['C_{2} = ' num2str(C2/nano) ' nF, ' ...
     'C_{4} = ' num2str(C4/nano) ' nF'])
xlabel('Frequency (Hz)')
ylabel('Magnitude')

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

% Choose K < 1 
K  = 1-(1e-14);

% Choose 0.1 < P < 10
P  = 0.5;

% Choose C2, C4, R5
C2 = 12*nano;
C4 = 3.3*nano;
R5 = 10*kilo;

disp(['wp = ' num2str(wp)])
disp(['Qp = ' num2str(Qp)])
disp(['K = ' num2str(K)])
disp(['P = ' num2str(P)])
disp('-------------------------')

% Compute element values
R1 = 1/(wp*sqrt(C2*C4*P));
R3 = P*R1;
R6 = R5*((1+P)*C4/C2-sqrt(P*C4/C2)/Qp);
Ko = 1+R6/R5;
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
R11 = str2num(num2str(R11,3)); 
R12 = str2num(num2str(R12,3));
R3  = str2num(num2str(R3,3));
R5  = str2num(num2str(R5,3));  
R6  = str2num(num2str(R6,3));
C2  = str2num(num2str(C2,3));  
C4  = str2num(num2str(C4,3));

disp(['R11 = ' num2str(R11/kilo) ' kOhm'])
disp(['R12 = ' num2str(R12/kilo) ' kOhm'])
disp(['R3 = ' num2str(R3/kilo) ' kOhm'])
disp(['R5 = ' num2str(R5/kilo) ' kOhm'])
disp(['R6 = ' num2str(R6/kilo) ' kOhm'])
disp(['C2 = ' num2str(C2/nano) ' nF'])
disp(['C4 = ' num2str(C4/nano) ' nF'])

% Validate realization
% Generate transfer function in terms of element values
num2 = R12*(R5+R6);
den2 = [(C2*C4*R11*R12*R3*R5) ...
       (C4*R11*R12*R5+C4*R11*R3*R5+C4*R12*R3*R5-C2*R11*R12*R6) ...
       (R11*R5+R12*R5)];
H2r = freqs(num2,den2,2*pi*f);

figure
subplot(2,1,1)
drawlpmq0(0,0,4,5,9,'b','           Biquad #2, LP-MQ')
drawfout(21,6,4,5,11,'m')
drawfin(6,7,4,5,11,'m')

subplot(2,1,2)
plot(f,abs(H2r),f(1:10:end),abs(H2(1:10:end)),'x')
text(200, 3.5, ...
    ['C_{2} = ' num2str(C2/nano) ' nF, ' ...
     'C_{4} = ' num2str(C4/nano) ' nF'])
text(200, 2.8, ...
    ['R_{11} = ' num2str(R11/kilo) ' k\Omega, ' ...
     'R_{12} = ' num2str(R12/kilo) ' k\Omega'])
text(200, 2.1, ...
    ['R_{5} = ' num2str(R5/kilo) ' k\Omega, ' ...
     'R_{6} = ' num2str(R6/kilo) ' k\Omega'])
text(200, 1.4, ...
    ['R_{3} = ' num2str(R3/kilo) ' k\Omega'])
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:5:end),A1(1:5:end)+A2(1:5:end),'mx',...
     f,A1+A2-(A1r+A2r),'r-')
axis(1.1*[0 Fp -Ap Ap])
legend('realized','ideal','error')
xlabel('Frequency (Hz)')
ylabel('Attenuation (dB)')
