The Complete Matlab/Octave Source Code for Design of Unassisted Alignments

function params = calc_vented_box_params(qt, ql)
  %CALC_VENTED_BOX_PARAMS Compute h and alpha given QT and QL
  theta_a = pi / 8;
  theta_b = 3 * pi / 8;
  % System order
  n = 4;
  % Compute Butterworth QT = qtb
  qtb = 1 / (2 * (cos(theta_a) + cos(theta_b)) - 1 / ql);
  if qt >= qtb
    % Butterworth and Chebyshev alignments
    params = find_system_params_cheby_butter(qt, ql, theta_a, theta_b, n);
  else
    % Quasi-Butterworth alignment
    params = find_system_params_quasi_butter(qt, ql);
  end
end

function params = find_system_params_cheby_butter(qt, ql, theta_a, theta_b, n)
  a = 2 * (cos(theta_a) + cos(theta_b));
  b = 2 * (cos(theta_a) * cos(theta_b)^2 + cos(theta_b) * cos(theta_a)^2);
  c = 2 * (cos(theta_a) * sin(theta_b)^2 + cos(theta_b) * sin(theta_a)^2);
  d = cos(theta_a)^2 + cos(theta_b)^2 + 4 * cos(theta_a) * cos(theta_b);
  e = sin(theta_a)^2 + sin(theta_b)^2;
  f = cos(theta_a)^2 * cos(theta_b)^2;
  g = cos(theta_a)^2 * sin(theta_b)^2 + sin(theta_a)^2 * cos(theta_b)^2;
  h = sin(theta_a)^2 * sin(theta_b)^2;
  guess_for_k = 0.95;
  k = fzero(@(k)func_to_solve_cheby_butter(k, qt, ql, a, b, c, d, e, f, g, h), guess_for_k, get_fzero_opts());
  coeffs = calc_tf_coeffs_cheby_butter(k, a, b, c, d, e, f, g, h);
  a1 = coeffs.a1;
  a2 = coeffs.a2;
  a3 = coeffs.a3;
  sqrt_h = calc_sqrt_h(qt, ql, a1, a3);
  hval = sqrt_h^2;
  alpha = calc_alpha(qt, ql, a2, hval);
  params.h = hval;
  params.alpha = alpha;
  if k < 1
    epsi = 1 / sinh(n * atanh(k));
    params.dBripple = 10 * log10(1 + epsi^2);
    mu = 1 / cosh(1 / n * asinh(1 / epsi));
    w3LP = mu * cosh(1 / n * acosh(sqrt(2 + 1 / epsi^2)));
    w1 = coeffs.w1;
    params.f3_norm = w1 * sqrt_h / w3LP;
  else
    % Should only get here in Butterworth case
    params.dBripple = 0;
    params.f3_norm = 1.0;
  end
end

function options = get_fzero_opts()
  % Forces the fzero() termination tolerance for
  % Octave to be the same as for Matlab
  options = optimset('fzero');
  options.TolX = eps;
end

function to_zero = func_to_solve_cheby_butter(k, qt, ql, a, b, c, d, e, f, g, h)
  coeffs = calc_tf_coeffs_cheby_butter(k, a, b, c, d, e, f, g, h);
  a1 = coeffs.a1;
  a3 = coeffs.a3;
  sqrth = calc_sqrt_h(qt, ql, a1, a3);
  to_zero = sqrth^2 - a1 * ql * sqrth + ql / qt;
end

function coeffs = calc_tf_coeffs_cheby_butter(k, a, b, c, d, e, f, g, h)
  w1 = (k^4 * f + k^2 * g + h)^0.25;
  coeffs.a1 = (k^3 * b + k * c) / w1^3;
  coeffs.a2 = (k^2 * d + e) / w1^2;
  coeffs.a3 = k * a / w1;
  coeffs.w1 = w1;
end

function params = find_system_params_quasi_butter(qt, ql)
  guess_for_a2 = 3.4142136; %a2 value for fourth order Butterworth
  a2 = fzero(@(a2)func_to_solve_quasi_butter(a2, qt, ql), guess_for_a2, get_fzero_opts());
  a1 = calc_a1_quasi_butter(a2);
  a3 = calc_a3_quasi_butter(a2);
  sqrt_h = calc_sqrt_h(qt, ql, a1, a3);
  hval = sqrt_h^2;
  alpha = calc_alpha(qt, ql, a2, hval);
  params.h = hval;
  params.alpha = alpha;
  params.dBripple = 0;
  params.f3_norm = find_f3_norm(a1, a2, a3, hval);
end

function to_zero = func_to_solve_quasi_butter(a2, qt, ql)
  a1 = calc_a1_quasi_butter(a2);
  a3 = calc_a3_quasi_butter(a2);
  sqrth = calc_sqrt_h(qt, ql, a1, a3);
  to_zero = sqrth / qt + 1 / (ql * sqrth) - a3;
end

function sqrt_h = calc_sqrt_h(qt, ql, a1, a3)
  sqrt_h = (a3 * ql - a1 * qt) / (ql / qt - qt / ql);
end

function alpha = calc_alpha(qt, ql, a2, hval)
  alpha = hval * (ql * qt * a2 - 1) / (ql * qt) - 1 - hval^2;
end

function a1 = calc_a1_quasi_butter(a2)
  a1 = sqrt(2 * a2);
end

function a3 = calc_a3_quasi_butter(a2)
  a3 = (a2^2 + 2) / sqrt(8 * a2);
end

function filt = G(s, a1, a2, a3, h)
  % normalize the frequency
  sn = s / sqrt(h);
  filt = sn^4 / (sn^4 + a1 * sn^3 + a2 * sn^2 + a3 * sn + 1);
end

function resp = mag_squared_of_G(w, a1, a2, a3, h)
  cplx_G = G(j * w, a1, a2, a3, h);
  resp = real(cplx_G)^2 + imag(cplx_G)^2;
end

function to_zero = cutoff_freq_of_G(w, a1, a2, a3, h)
  to_zero = mag_squared_of_G(w, a1, a2, a3, h) - 0.5;
end

function w = find_f3_norm(a1, a2, a3, h)
  guess_for_w = sqrt(h);
  w = fzero(@(w)cutoff_freq_of_G(w, a1, a2, a3, h), guess_for_w, get_fzero_opts());
end