% Bessel filters
freq = 10.^(-3:0.001:3);

% Begin with a clean graphical slate...
for fig=1:6 %bodeplot creates 6 figures.
   figure(fig)
   hold off
endfor


% 2-pole Bessel
Q = [0.577];
fn = [1.273];
Bessel2 = multistage(freq, Q, fn);
bodeplot(freq, Bessel2, 'r')

% 4-pole Bessel
Q = [0.522, 0.806];
fn = [1.432, 1.606];
Bessel4 = multistage(freq, Q, fn);
bodeplot(freq, Bessel4, 'g')

% 6-pole Bessel
Q = [0.510, 0.611, 1.024];
fn = [1.607, 1.692, 1.908];
Bessel6 = multistage(freq, Q, fn);
bodeplot(freq, Bessel6, 'b')

% 8-pole Bessel
Q = [0.506, 0.560, 0.711, 1.225];
fn = [1.781, 1.835, 1.956, 2.192];
Bessel8 = multistage(freq, Q, fn);
bodeplot(freq, Bessel8, 'm')

legendtext=['2-pole';'4-pole';'6-pole';...
   '8-pole'];

figure(1)
title('Gain of Bessel filters','fontsize',18)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Gain (dB)','fontsize',16)
legend(legendtext)
axis([0.01,1000,-160,2])
print('BesselGain.pdf')
axis([0.1,10,-20,5])
print('BesselGainZoom.pdf')

figure(2)
title('Phase of Bessel filters','fontsize',18)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Phase (degrees)','fontsize',16)
legend(legendtext)
print('BesselPhase.pdf')

figure(3)
title('Time delay of Bessel filters','fontsize',18)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Time delay (fc = 1)','fontsize',16)
legend(legendtext)
print('BesselTimeDelay.pdf')

figure(4)
title('Time delay vs. attenuation of Bessel filters','fontsize',18)
xlabel('Attenuation (dB)','fontsize',16)
ylabel('Time delay (s)','fontsize',16)
legend(legendtext)
axis([-3,10])
print('BesselTimeDelay-Attenuation.pdf')

figure(5)
title('Step response of Bessel filters','fontsize',18)
xlabel('Time (s)','fontsize',16)
ylabel('Response','fontsize',16)
legend(legendtext)
axis([-1,2])
print('BesselStepResponse.pdf')

figure(6)
title('Impulse response of Bessel filters','fontsize',18)
xlabel('Time (s)','fontsize',16)
ylabel('Response','fontsize',16)
legend(legendtext)
axis([-1,2])
print('BesselImpulseResponse.pdf')
