% Example highpass 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 = 5000;
Fs = 3400;
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, 'high', 's');

% Decompose transfer function into biquads

% MATLAB 7 sort supports 'descend'
% p = sort(p,'descend');

% Here is the code for older versions:
pSortAscend = sort(p);
p = pSortAscend(end:-1:1);

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 = [1 0 0];
N2 = [1 0 0];

% Plot attenuation
Fplot = 16*1000;
f = 2000: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(4000, 35, ['\omega_1 = ' num2str(round(w(1))) ',  Q_1 = ' num2str(Q(1))],'color','r')
text(4000, 22, ['\omega_2 = ' num2str(round(w(3))) ',  Q_2 = ' num2str(Q(3))],'color','b')

% Choose realization from album
figure 
drawhplq0(0,4,4,5,8,'r','\bf Biquad #1')
drawhpmq0(15,0,4,5,8,'b','\bf Biquad #2')
drawfout(32,7,4,5,8,'m')
drawfin(2,9,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-10^(-15);

% Choose C1
C1 = 3.3*nano;

% Choose P to achieve C12 = 0
P  = 4*Qp^2; 

disp(['wp = ' num2str(wp)])
disp(['Qp = ' num2str(Qp)])
disp(['K = ' num2str(K)])
disp(['C1 = ' num2str(C1/nano) ' nF'])
disp(['P = ' num2str(P) ',   0.1 < P < 10'])
disp('-------------------------')

% Compute element values
C11 = C1*K;
C12 = C1-C11;
C3 = C1*(P-2*Qp^2 - sqrt(P^2-4*P*Qp^2))/(2*Qp^2);
R2 = 1/(wp*(C1+C3)*Qp);
R4 = 1/(C1*C3*R2*wp^2);

% Small values of C12 are interpreted as open-circuited branches
if abs(C12)<10^(-12); C12 = 0; end
% Small C12 does not appear in implementation!

% Round element values to 3 digits of precision
C11 = str2num(num2str(C11,3)); 
C12 = str2num(num2str(C12,3));
C3  = str2num(num2str(C3,3));
R2  = str2num(num2str(R2,3)); 
R4  = str2num(num2str(R4,3));

disp(['C11 = ' num2str(C11/nano) ' nF'])
disp(['C12 = ' num2str(C12/nano) ' nF'])
disp(['C3 = ' num2str(C3/nano) ' nF'])
disp(['R2 = ' num2str(R2/kilo) ' kOhm'])
disp(['R4 = ' num2str(R4/kilo) ' kOhm'])

% Validate realization
% Generate transfer function in terms of element values
num = [(C11*C3*R2*R4) 0 0];
den = [(C11*C3*R2*R4+C12*C3*R2*R4) (C11*R2+C12*R2+C3*R2) 1];
H1r = freqs(num,den,2*pi*f);

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

subplot(2,1,2)
plot(f,abs(H1r),'r',f(1:10:end),abs(H1(1:10:end)),'rx')
text(8000, 0.5, ...
    ['C_{11} = ' num2str(C11/nano) ' nF, ' ...
     'C_{12} = ' num2str(C12/nano) ' nF, ' ...
     'C_{3} = ' num2str(C3/nano) ' nF'])
text(8000, 0.2, ...
    ['R_{2} = ' num2str(R2/kilo) ' k\Omega, ' ...
     'R_{4} = ' num2str(R4/kilo) ' k\Omega'])
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 > 1 
K  = 1.01;

% Choose P
Pinit  = 10; % initial guess

% Choose C1, C3, R5
C1 = 3.3*nano;
C3 = 3.3*nano;
R5 = 6800;

% Optimize P to achieve C12 = 0
y = inline('-C1*P*sqrt(C1*C3)+sqrt(abs(P))*(Qp*C1^2+C1*C3*Qp+C1*C3*P*Qp-C1*C3*K*P*Qp)',...
           'P','C1','C3','Qp','K');

% MATLAB 7 syntax
% P = fzero(y,Pinit,10,C1,C3,Qp,K);

% MATLAB 6.5 syntax
[P, objectiveFunctionValue] = fzero(y,Pinit,[],C1,C3,Qp,K);

if isnan(P); P=Pinit; end

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

% Compute element values
R2 = 1/(wp*sqrt(C1*C3*P));
R4 = P*R2;
R6 = R5*((1+C1/C3)/P-sqrt(C1/(P*C3))/Qp);
Ko = 1+R6/R5;
C11 = C1*K/Ko;
C12 = C1-C11; if abs(C12)<10^(-12); C12 = 0; end

% Round element values to 3 digits of precision
C11 = str2num(num2str(C11,3)); 
C12 = str2num(num2str(C12,3));
C3  = str2num(num2str(C3,3));
R2  = str2num(num2str(R2,3));  
R4  = str2num(num2str(R4,3));
R5  = str2num(num2str(R5,3));  
R6  = str2num(num2str(R6,3));

disp(['C11 = ' num2str(C11/nano) ' nF'])
disp(['C12 = ' num2str(C12/nano) ' nF'])
disp(['C3 = ' num2str(C3/nano) ' nF'])
disp(['R5 = ' num2str(R5/kilo) ' kOhm'])
disp(['R6 = ' num2str(R6/kilo) ' kOhm'])
disp(['R2 = ' num2str(R2/kilo) ' kOhm'])
disp(['R4 = ' num2str(R4/kilo) ' kOhm'])

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

figure
subplot(2,1,1)
drawhpmq0(0,0,4,5,9,'b','Biquad #2, HP-MQ')
drawfout(17,7,4,5,11,'m')
drawfin(2,8,4,5,11,'m')

subplot(2,1,2)
plot(f,abs(H2r),'b',f(1:10:end),abs(H2(1:10:end)),'bx')
text(4200, 0.95, ...
    ['C_{11} = ' num2str(C11/nano) ' nF, ' ...
     'C_{12} = ' num2str(C12/nano) ' nF, ' ...
     'C_{3} = ' num2str(C3/nano) ' nF'])     
text(4200, 0.4, ...
    ['R_{2} = ' num2str(R2/kilo) ' k\Omega, ' ...
     'R_{4} = ' num2str(R4/kilo) ' k\Omega, ' ...
     'R_{5} = ' num2str(R5/kilo) ' k\Omega, ' ...
     'R_{6} = ' num2str(R6) ' \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([Fs Fplot -1.2*Ap 1.2*Ap])
legend('realized','ideal','error')
xlabel('Frequency (Hz)')
ylabel('Attenuation (dB)')
