% Example filter specification and design
% Course: Analog and Digital Filter Design
% Lecture: Specification of an analog filter
% 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/01/30$
% 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, format short e

% spec

Fp = 3400;
Fs = 5000;
Ap = 1;
As = 20;
Fplot = 6000;

plot([0  Fp Fp],[Ap Ap As],'r--',...
     [0 Fs Fs Fplot],[-Ap -Ap As As],'r--')
axis([0 Fplot -2*Ap 1.2*As])
ytick1 = [0 Ap As]; set(gca,'YTick',ytick1)
xtick1 = [0 Fp Fs]; set(gca,'XTick',xtick1)
xlabel('f (Hz)')
ylabel('A (dB)')
text(Fp/10,As/2,['S = [' num2str([Fp Fs Ap As]) ']'])
title('Filter specification')

% step 1
% select approximation Chev
% find order, poles and zeros, tf
[n, Wn] = cheb1ord(Fp, Fs, Ap, As, 's');
[z,p,k] = cheb1ap(n, Ap);

Fn = Fp; 
num = 10^(-Ap/20); den = 1;
for i = 1:length(p)
    if imag(p(i)) >0
        Q(i) = (-abs(p(i))./(2*real(p(i))));
        w(i) = round(abs(p(i))*(2*pi*Fn));
        D(i,:) = [1 w(i)/Q(i) w(i)^2];
        N(i,1) = w(i)^2;
        num = conv(num,N(i,:));
        den = conv(den,D(i,:));
    end
end

f = 0:Fplot/200:Fplot;
H = freqs(num,den,2*pi*f);
M = -20*log10(abs(H));
figure
plot(f,M,[0 Fp Fp],[Ap Ap As],'r--',[0 Fs Fs Fplot],[-Ap -Ap As As],'r--')
xlabel('f (Hz)')
ylabel('A (dB)')
text(500, As-5, ['\omega_1 = ' num2str(w(1)) ',  Q_1 = ' num2str(Q(1))])
text(500, As, ['\omega_2 = ' num2str(w(2)) ',  Q_2 = ' num2str(Q(2))])
title('Design step 1: approximation')

% step 2: realization
% Chose cascade of 2nd order programmable analog filters
figure
drawschematic(0,0,4,5,9,'b')
text(5, 15,['\omega_1 = ' num2str(w(1)) ',  Q_1 = ' num2str(Q(1)) ',   '...
            '\omega_2 = ' num2str(w(2)) ',  Q_2 = ' num2str(Q(2))])
title('Design step 2: realization')

% step 3: study of imperfections
% Sensitivity analysis
figure
tolerance = 4;
num1 = 10^(-Ap/20); den1 = 1;
num2 = 10^(-Ap/20); den2 = 1;
Q1 = (1 + tolerance/100)*Q;
w1 = (1 + tolerance/100)*w;
Q2 = (1 - tolerance/100)*Q;
w2 = (1 - tolerance/100)*w;
for i = 1:length(Q)
    num1 = conv(num1,w1(i)^2);
    den1 = conv(den1,[1 w1(i)/Q1(i) w1(i)^2]);
    num2 = conv(num2,w2(i)^2);
    den2 = conv(den2,[1 w2(i)/Q2(i) w2(i)^2]);
end
H1 = freqs(num1,den1,2*pi*f);
M1 = -20*log10(abs(H1));
H2 = freqs(num2,den2,2*pi*f);
M2 = -20*log10(abs(H2));
plot(f,M,f,M1,'r--',f,M2,'r:',[0 Fp Fp],[Ap Ap As],'r--',[0 Fs Fs Fplot],[-Ap -Ap As As],'r--')
xlabel('f (Hz)')
ylabel('A (dB)')
text(500, As-5, ['\omega_1 = ' num2str(w(1)) ',  Q_1 = ' num2str(Q(1))])
text(500, As, ['\omega_2 = ' num2str(w(2)) ',  Q_2 = ' num2str(Q(2))])
title('Design step 3: study of imperfections')

% step 1 re-design
% select new approximation Chev
% increase order, find poles and zeros, tf
n = n+2;
ap = Ap/25;
Fn = 1.12*Fp; 

[z,p,k] = cheb1ap(n+2, ap);
num = 10^(-ap/20); den = 1;
for i = 1:length(p)
    if imag(p(i)) >0
        Q(i) = (-abs(p(i))./(2*real(p(i))));
        w(i) = round(abs(p(i))*(2*pi*Fn));
        D(i,:) = [1 w(i)/Q(i) w(i)^2];
        N(i,1) = w(i)^2;
        num = conv(num,N(i,:));
        den = conv(den,D(i,:));
    end
end

f = 0:Fplot/200:Fplot;
H = freqs(num,den,2*pi*f);
M = -20*log10(abs(H));
figure
plot(f,M,[0 Fp Fp],[Ap Ap As],'r--',[0 Fs Fs Fplot],[-Ap -Ap As As],'r--')
xlabel('f (Hz)')
ylabel('A (dB)')
text(500, As-5, ['\omega_1 = ' num2str(w(1)) ',  Q_1 = ' num2str(Q(1))])
text(500, As, ['\omega_2 = ' num2str(w(2)) ',  Q_2 = ' num2str(Q(2))])
text(500, As+5, ['\omega_3 = ' num2str(w(3)) ',  Q_3 = ' num2str(Q(3))])
title('Re-design step 1: approximation')

% step 3: study of imperfections
% Sensitivity analysis
figure
num1 = 10^(-ap/20); den1 = 1;
num2 = 10^(-ap/20); den2 = 1;
Q1 = (1 + tolerance/100)*Q;
w1 = (1 + tolerance/100)*w;
Q2 = (1 - tolerance/100)*Q;
w2 = (1 - tolerance/100)*w;
for i = 1:length(Q)
    num1 = conv(num1,w1(i)^2);
    den1 = conv(den1,[1 w1(i)/Q1(i) w1(i)^2]);
    num2 = conv(num2,w2(i)^2);
    den2 = conv(den2,[1 w2(i)/Q2(i) w2(i)^2]);
end
H1 = freqs(num1,den1,2*pi*f);
M1 = -20*log10(abs(H1));
H2 = freqs(num2,den2,2*pi*f);
M2 = -20*log10(abs(H2));
plot(f,M,f,M1,'g--',f,M2,'g:',[0 Fp Fp],[Ap Ap As],'r--',[0 Fs Fs Fplot],[-Ap -Ap As As],'r--')
xlabel('f (Hz)')
ylabel('A (dB)')
text(500, As-5, ['\omega_1 = ' num2str(w(1)) ',  Q_1 = ' num2str(Q(1))])
text(500, As, ['\omega_2 = ' num2str(w(2)) ',  Q_2 = ' num2str(Q(2))])
text(500, As+5, ['\omega_3 = ' num2str(w(3)) ',  Q_3 = ' num2str(Q(3))])
title('Re-design step 3: study of imperfections')

% step 4: implementation
% Chose cascade of 2nd order programmable analog filters
figure
drawschematic3(0,0,4,5,8,'b')
text(0, 15,['\omega_1 = ' num2str(w(1)) ',  Q_1 = ' num2str(Q(1)) ',   '...
            '\omega_2 = ' num2str(w(2)) ',  Q_2 = ' num2str(Q(2)) ',   '...
            '\omega_3 = ' num2str(w(3)) ',  Q_3 = ' num2str(Q(3))])
title('Design step 4: implementation')
