% Example general purpose 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.00 $  $Date: 2005/02/13$
% 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 Filt to the MATLAB path
path(path,'C:\afd\DrawFilt')

kilo = 10^3;
mega = 10^6;
nano = 10^(-9);

% Given Q-factor and pole magnitude
fp = 2500;
wp = 2*pi*fp;
Qp = 1.5;

% Choose Cx and Rx
Cx = 10*nano;
Rx = 10*kilo;

disp(['wp = ' num2str(wp)])
disp(['Qp = ' num2str(Qp)])
disp(['Rx = ' num2str(Rx/kilo) ' kOhm'])
disp(['Cx = ' num2str(Cx/nano) ' nF'])
disp('-------------------------')

% Compute element values
C6 = Cx;
C8 = Cx;
Ro = 1/(wp*Cx);
R1 = Rx;
R3 = Rx;
R5 = Rx;
R7 = Rx;
R4 = Rx*(Rx/Ro)^2;
R2 = Rx*(Qp*(1+R4/Rx)/sqrt(R4/Rx)-1);
R9 = 2*Rx;
R10= Rx;
R11= Rx;


% Round element values to 3 digits of precision
R1 = str2num(num2str(R1,3)); 
R2 = str2num(num2str(R2,3));
R3 = str2num(num2str(R3,3)); 
R4 = str2num(num2str(R4,3));
R5 = str2num(num2str(R5,3)); 
R7 = str2num(num2str(R7,3));
C6 = str2num(num2str(C6,3)); 
C8 = str2num(num2str(C8,3));

disp(['R1 = ' num2str(R1/kilo) ' kOhm'])
disp(['R2 = ' num2str(R2/kilo) ' kOhm'])
disp(['R3 = ' num2str(R3/kilo) ' kOhm'])
disp(['R4 = ' num2str(R4/kilo) ' kOhm'])
disp(['R5 = ' num2str(R5/kilo) ' kOhm'])
disp(['R7 = ' num2str(R7/kilo) ' kOhm'])
disp(['C6 = ' num2str(C6/nano) ' nF'])
disp(['C8 = ' num2str(C8/nano) ' nF'])

% Generate transfer function in terms of element values
numLP = [0 0 R2*(R3+R4)];
numHP = [C6*C8*R2*(R3+R4)*R5*R7 0 0];
numBP = [0 C8*R2*(R3+R4)*R7 0];
numN = [C6*C8*R2*(R3+R4)*R5*R7*R11/R9 0 R2*(R3+R4)*R11/R10];

den = [(C6*C8*R1*R3*R5*R7+C6*C8*R2*R3*R5*R7) ...
       (C8*R1*R3*R7+C8*R1*R4*R7) ...
       (R1*R4+R2*R4)];
f = 0:fp/50:3.2*fp;
HLP = freqs(numLP,den,2*pi*f);
HHP = freqs(numHP,den,2*pi*f);
HBP = freqs(numBP,den,2*pi*f);
HN = freqs(numN,den,2*pi*f);

drawkhn0(0,0,4,5,9,'b')

figure
plot(f,abs(HLP),'b',f,abs(HHP),'r',f,abs(HBP),'g',f,abs(HN),'m:')
text(200, 3.5, ['C_{6}= C_{8}=' num2str(Cx/nano) ' nF'])
text(200, 3, ['R_{1}=R_{3}=R_{5}=' num2str(Rx/kilo) ' k\Omega'])
text(200, 2.5, ['R_{2}=' num2str(R2/kilo) ' k\Omega'])
text(200, 2.0, ['R_{4}=' num2str(R4/kilo) ' k\Omega'])
text(200, 1.5, ['R_{7}=' num2str(Rx/kilo) ' k\Omega'])

xlabel('Frequency (Hz)')
ylabel('Magnitude')
legend('LP','HP','BP','Notch')
