close all;

fp = 1e5;
wp = 2*pi*fp;

f = logspace(4, 7, 400);
w = 2*pi*f;
flin = linspace(0, 2*fp, 400);
wlin = 2*pi*flin;

[ N, wn ] = buttord(wp, 10*wp, 3, 80, 's');
[b, a] = butter(N, wn, 's');
but = tf(b, a);

[ N, wn ] = cheb2ord(wp, 10*wp, 3, 80, 's');
[b, a] = cheby2(N, 80, wn, 's');
cheby4 = tf(b, a);

[ N, wn ] = ellipord(wp, 10*wp, 3, 80, 's');
[b, a] = ellip(N, 3, 80, wn, 's');
elliptic = tf(b, a);

[b, a] = besself(4, 0.76*1.52*wp);
bessel4 = tf(b, a);

[b, a] = besself(6, 0.75*1.75*wp);
bessel6 = tf(b, a);

[b, a] = besself(8, 0.74*1.95*wp);
bessel8 = tf(b, a);

% allpass (from E. Swanson)

fp1 = 14523;    Q1 = 0.56707;
fp2 = 30495;    Q2 = 1.9341;

fp1 = 5*fp1;    fp2 = 5*fp2;
s = tf('s');

w1 = 2*pi*fp1;  w2 = 2*pi*fp2;

allpass1 = (1 - s/(w1*Q1) + (s/w1)^2) / (1 + s/(w1*Q1) + (s/w1)^2);
allpass2 = (1 - s/(w2*Q2) + (s/w2)^2) / (1 + s/(w2*Q2) + (s/w2)^2);
allpass = allpass1 * allpass2;
cheby4allpass = cheby4*allpass;

% butter and cheby

bodemaghz(but, w);
hold on;
bodemaghz(cheby4, w);
bodemaghz(elliptic, w);
legend('Butterworth', 'Chebychev 2', 'Elliptic', 3);
fixfig;
hold off;

figure;
t = linspace(0, 2e-5, 100);
step(but, t);
hold on;
step(cheby4, t);
step(elliptic, t);
legend('Butterworth', 'Chebychev 2', 'Elliptic', 4);
fixfig;
hold off;

% magnitude response

figure;
bodemaghz(cheby4, w);
hold on;
bodemaghz(bessel4, w);

axis([ min(f) max(f) -100 5 ]);
% axis([ 0.9e5 1.1e5 -5 -1 ]);
legend('4th Order Chebychev 2', '4th Order Bessel', 3);
fixfig;
hold off;

% phase response

figure;
bodephhz(cheby4, wlin);
hold on;
bodephhz(bessel4, wlin);

legend('4th Order Chebychev 2', '4th Order Bessel', 3);
fixfig;
hold off;

% step response

figure;
t = linspace(0, 2e-5, 100);
step(cheby4, t);
hold on;
step(bessel4, t);
legend('4th Order Chebychev 2', '4th Order Bessel');
fixfig;
hold off;

% group delay

figure;
d1 = grpdelay(cheby4,  w);
d2 = grpdelay(bessel4, w);
d3 = grpdelay(bessel6, w);
d4 = grpdelay(bessel8, w);
semilogx(f, 1e6*d1, f, 1e6*d2);
legend('4th Order Chebychev 2', '4th Order Bessel');
xlabel('Frequency [Hz]');
ylabel('Group Delay [\mu s]');
fixfig;

% normalized delay

figure;
d1 = d1 / d1(1);
d2 = d2 / d2(1);
d3 = d3 / d3(1);
d4 = d4 / d4(1);
semilogx(f, d1, f, d2);
legend('4th Order Chebychev 2', '4th Order Bessel');
xlabel('Frequency [Hz]');
ylabel('Group Delay [normalized]');
fixfig;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ALLPASS magnitude response

figure;
bodemaghz(cheby4, w);
hold on;
bodemaghz(allpass, w);
bodemaghz(cheby4allpass, w);
bodemaghz(bessel8, w);

axis([ min(f) max(f) -100 5 ]);
legend('4th Order Chebychev 2', 'Allpass', 'Cheby+Allpass', '8th Order Bessel', 3);
fixfig;
hold off;

% ALLPASS phase response

figure;
bodephhz(cheby4, wlin);
hold on;
bodephhz(allpass, wlin);
bodephhz(cheby4allpass, wlin);
bodephhz(bessel8, wlin);

legend('4th Order Chebychev 2', 'Allpass', 'Cheby+Allpass', '8th Order Bessel', 3);
fixfig;
hold off;

% ALLPASS step response

figure;
t = linspace(0, 2e-5, 100);
step(cheby4, t*1.5);
hold on;
step(allpass, t*1.5);
step(cheby4allpass, t*1.5);
step(bessel8, t*1.5);
legend('4th Order Chebychev 2', 'Allpass', 'Cheby+Allpass', '8th Order Bessel', 4);
fixfig;
hold off;

% ALLPASS normalized delay

figure;
d5 = grpdelay(allpass, w);
d6 = grpdelay(cheby4allpass, w);
d5 = d5 / d5(1);
d6 = d6 / d6(1);
semilogx(f, d1, f, d6, f, d4);
legend('4th Order Chebychev 2', 'Cheby+Allpass', '8th Order Bessel');
xlabel('Frequency [Hz]');
ylabel('Group Delay [normalized]');
fixfig;
