% Main program to calculate group delay of a butterworth bessel
% Author: mik81@de.wikipedia.org
% bessel especially butterworth with 3rd Order

clear all;

omega = (-2.:0.01:2.);
omega = 10..^omega;
time = (0.:0.01:20.);
step = zeros(length(time),1);
for i=200:1:length(step)
  step(i) = 1;
end

%
% Functions
%

besselThirdOrder = 1 ./ \
   ( 1 .+ ( omega .*  1i ) \
      + ( omega .*  1i).^2 .* (2 ./ 5) \
      + ( omega .*  1i).^3 .* (1 ./ 15) \
   );
besselFourthOrder = 1 ./ \
  ( 1 .+ ( omega .* 1i ) \
    + ( omega .* 1i).^2 .* (3 ./ 7) \
    + ( omega .* 1i).^3 .* (2 ./ 21) \
    + ( omega .* 1i).^4 .* (1 ./ 105) \
  );
besselFifthOrder = 1 ./ \
  ( 1 .+ ( omega .* 1i ) \
    + ( omega .* 1i ).^2 .* (4 ./ 9) \
    + ( omega .* 1i ).^3 .* (1 ./ 9) \
    + ( omega .* 1i ).^4 .* (1 ./ 63) \
    + ( omega .* 1i ).^5 .* (1 ./ 945) \
  );

butterworthThirdOrder = 1 ./ \
   ( (1 .+ (omega .* (1/1.8) .* 1i)) \
     .* ( 1 .+ (omega .* (1/1.8) .* 1i) .+ (omega .* (1/1.8) .* 1i).^2 ));
% butterworthThirdOrderIIR_y_a = 1 ./ (1+ y_a_1 )
% butterworthThirdOrderIIR_y_b = ( 1 ./ ( 1 .+ y_b_1 .+ y_b_2^2 ))
threeFilterFirstOrder = 1 ./ \
  ( (1 .+ ( 1/3.6 .* omega .* 1i)) \
     .* (1 .+ ( 1/3.6 .* omega .* 1i)) \
     .* (1 .+ ( 1/3.6 .* omega .* 1i)) \
  );

%
% Calculate and plot gain
%

xlabel("Normalized omega / omega_0");
ylabel("Gain");
loglog( \
   omega, abs(butterworthThirdOrder), ";Butterworth 3. Ordnung;" \
   , omega, abs(besselThirdOrder),  ";Besselfilter 3. Ordnung;" \
%   , omega, abs(besselFourthOrder),  ";Besselfilter 4. Ordnung;" \
  , omega, abs(threeFilterFirstOrder), ";3 Tiefpaesse 1. Ordnung;" \
   , [0.01, 100.], [1/sqrt(2), 1/sqrt(2)], ";-3dB;");
axis ([0.01, 100., 0.00001, 5.]);
grid ("on");

print ("butterworthFrequency.png", "-dpng", "-landscape");

%
% Calculate and plot phase
%

% Butterworth
butterworthThirdOrderAngle = angle(butterworthThirdOrder);
% correct roll over of phase
shift = 0;
lastAngle = butterworthThirdOrderAngle(1);

for i=2:1:length(butterworthThirdOrderAngle)
  if ( lastAngle < 0 && butterworthThirdOrderAngle(i) > 0)
    shift = shift - 2*pi;
  end
  lastAngle = butterworthThirdOrderAngle(i);
  butterworthThirdOrderAngle(i) = butterworthThirdOrderAngle(i) + shift;
end

% Bessel

besselThirdOrderAngle = angle(besselThirdOrder);

% correct roll over of phase
shift = 0;
lastAngle = besselThirdOrderAngle(1);
for i=2:1:length(besselThirdOrderAngle)
  if ( lastAngle < 0 && besselThirdOrderAngle(i) > 0)
    shift = shift - 2*pi;
  end
  lastAngle = besselThirdOrderAngle(i);
  besselThirdOrderAngle(i) = besselThirdOrderAngle(i) + shift;
end
besselFourthOrderAngle = angle(besselFourthOrder);
% correct roll over of phase
shift = 0;
lastAngle = besselFourthOrderAngle(1);

for i=2:1:length(besselFourthOrderAngle)
  if ( lastAngle < 0 && besselFourthOrderAngle(i) > 0)
    shift = shift - 2*pi;
  end
  lastAngle = besselFourthOrderAngle(i);
  besselFourthOrderAngle(i) = besselFourthOrderAngle(i) + shift;
end

% three filter first order

threeFilterFirstOrderAngle = angle(threeFilterFirstOrder);

% correct roll over of phase
shift = 0;
lastAngle = threeFilterFirstOrderAngle(1);
for i=2:1:length(threeFilterFirstOrderAngle)
  if ( lastAngle < 0 && threeFilterFirstOrderAngle(i) > 0)
    shift = shift - 2*pi;
  end
  lastAngle = threeFilterFirstOrderAngle(i);
  threeFilterFirstOrderAngle(i) = threeFilterFirstOrderAngle(i) + shift;
end

xlabel("Normalized omega / omega_0");
ylabel("phase in rad");
% !!! unavoidable !!!
__gnuplot_set__ ytics pi/2;
semilogx( \
   omega, butterworthThirdOrderAngle, ";Butterworthfilter 3. Ordnung;"
   , omega, besselThirdOrderAngle, ";Besselfilter 3. Ordnung;"
%   , omega, besselFourthOrderAngle, ";Besselfilter 4. Ordnung;" );
   , omega, threeFilterFirstOrderAngle, ";3 Tiefpaesse 1. Ordnung;");
% legend (, );
axis ([0.01, 100., -pi*9/4, pi/4]);
grid ("on");

print ("butterworthAngle.png", "-dpng", "-landscape");

%
% Calculate and plot group delay with phase
%

%  -deltaAngle/deltaOmega

% butterworth
lastAngle = butterworthThirdOrderAngle(1);
lastOmega = omega(1);
butterworthThirdOrderGroupDelay = zeros(length(butterworthThirdOrder),1);
for i=2:1:length(butterworthThirdOrderAngle)
  butterworthThirdOrderGroupDelay(i) = \
     -1. * (butterworthThirdOrderAngle(i)-lastAngle)/(omega(i)-lastOmega);
  lastOmega = omega(i);
  lastAngle = butterworthThirdOrderAngle(i);
end
butterworthThirdOrderGroupDelay(1) = butterworthThirdOrderGroupDelay(2);

% bessel
lastAngle = besselThirdOrderAngle(1);
lastOmega = omega(1);
besselThirdOrderGroupDelay = zeros(length(besselThirdOrder),1);
for i=2:1:length(besselThirdOrderAngle)
  besselThirdOrderGroupDelay(i) = \
     -1. * (besselThirdOrderAngle(i)-lastAngle)/(omega(i)-lastOmega);
  lastOmega = omega(i);
  lastAngle = besselThirdOrderAngle(i);
end
besselThirdOrderGroupDelay(1) = besselThirdOrderGroupDelay(2);

lastAngle = besselFourthOrderAngle(1);
lastOmega = omega(1);
besselFourthOrderGroupDelay = zeros(length(besselFourthOrder),1);
for i=2:1:length(besselFourthOrderAngle)
  besselFourthOrderGroupDelay(i) = \
     -1. * (besselFourthOrderAngle(i)-lastAngle)/(omega(i)-lastOmega);
  lastOmega = omega(i);
  lastAngle = besselFourthOrderAngle(i);
end
besselFourthOrderGroupDelay(1) = besselFourthOrderGroupDelay(2);

% three first order
lastAngle = threeFilterFirstOrderAngle(1);
lastOmega = omega(1);
threeFilterFirstOrderGroupDelay = zeros(length(threeFilterFirstOrder),1);
for i=2:1:length(threeFilterFirstOrderAngle)
  threeFilterFirstOrderGroupDelay(i) = \
     -1. * (threeFilterFirstOrderAngle(i)-lastAngle)/(omega(i)-lastOmega);
  lastOmega = omega(i);
  lastAngle = threeFilterFirstOrderAngle(i);
end
threeFilterFirstOrderGroupDelay(1) = threeFilterFirstOrderGroupDelay(2);

xlabel("Normalized omega / omega_0");
ylabel("Groupdelay");
% !!! unavoidable !!!
__gnuplot_set__ ytics auto;

semilogx( \
   omega, butterworthThirdOrderGroupDelay, ";Butterworth 3. Ordnung;"
   , omega, besselThirdOrderGroupDelay, ";Bessel 3. Ordnung;"
%   , omega, besselFourthOrderGroupDelay, ";Besselfilter 4. Ordnung;");
  , omega, threeFilterFirstOrderGroupDelay, ";3 Tiefpaesse 1. Ordnung;");
axis ("auto");
grid ("on");

print ("butterworthGroupDelay.png", "-dpng", "-landscape");