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

function params = calc_vented_box_params_n(qt, ql, n, class)
  %CALC_VENTED_BOX_PARAMS_N Compute h, alpha and filter params given QT, QL, n and class
  angles = get_butterworth_pole_angles(n, class);
  theta_a = angles.theta_a;
  theta_b = angles.theta_b;
  theta_i = angles.theta_i;
  params = find_system_params(qt, ql, n, theta_a, theta_b, theta_i);
end

function angles = get_butterworth_pole_angles(n, class)
  if n < 5
    error('calc_vented_box_params_n is only for assisted alignments (n >= 5)');
  end
  is_odd = isodd(n);
  % Find the number of second-quadrant complex poles
  if is_odd
    ang_count = (n - 1) / 2;
  else
    ang_count = n / 2;
  end
  % Make sure the class spec has the right number of indices
  class_size = size(class, 2);
  if class_size < (ang_count - 2)
    error('class array has too few elements');
  end
  if class_size > (ang_count - 2)
    error('class array has too many elements');
  end
  if n == 5
    theta_a = pi / 5;
    theta_b = 2 * pi / 5;
    theta_i = [];
  else
    % Calculate the Butterworth pole angles
    theta_all = zeros(1, ang_count);
    for i = 1:ang_count
      if i == 1
        if is_odd
          theta_all(1) = pi / n;
        else
          theta_all(1) = pi / (2 * n);
        end
      else
        theta_all(i) = theta_all(1) + (i - 1) * pi / n;
      end
    end
    % Sort the Butterworth pole angles in descending order
    theta_all = sort(theta_all, 'descend');
    % Partition the Butterworth pole angles between loudspeaker and filter transfer functions
    spk_pole_count = 0;
    filt_pole_count = 0;
    for i = 1:ang_count
      index = find(class == i);
      if isempty(index)
        % Not found - must be a speaker pole
        spk_pole_count = spk_pole_count + 1;
        theta_ab(spk_pole_count) = theta_all(i);
      else
        % Found - must be a filter pole
        filt_pole_count = filt_pole_count + 1;
        theta_i(filt_pole_count) = theta_all(i);
      end
    end
    if size(theta_ab, 2) ~= 2 || size(theta_i, 2) ~= class_size
      error('Error in class array: may be some out-of-range, invalid or repeated indices');
    end
    theta_ab = sort(theta_ab);
    theta_a = theta_ab(1);
    theta_b = theta_ab(2);
  end
  angles.theta_a = theta_a;
  angles.theta_b = theta_b;
  angles.theta_i = theta_i;
end

function params = find_system_params(qt, ql, n, theta_a, theta_b, theta_i)
  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;
  qtb = calc_qtb_n(theta_a, theta_b, ql);
  % guess_for_k = 1.0;
  if qt < qtb
    guess_for_k = 1.5;
  elseif qt > qtb
    guess_for_k = 0.5;
  else
    guess_for_k = 1.0;
  end
  k = fzero(@(k)func_to_solve_k(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;
  w1 = coeffs.w1;
  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.k = k;
  if k < 1
    % Chebyshev case
    epsi = 1 / sinh(n * atanh(k));
    params.dBripple = 10 * log10(1 + epsi^2);
    mu = 1 / cosh(1 / n * asinh(1 / epsi));
    if isodd(n)
      w3LP = mu * cosh(1 / n * acosh(1 / epsi));
    else
      w3LP = mu * cosh(1 / n * acosh(sqrt(2 + 1 / epsi^2)));
    end
    params.f3_norm = w1 * sqrt_h / w3LP;
  elseif k == 1
    % Butterworth case
    params.dBripple = 0;
    params.f3_norm = 1;
  else
    % sub-Chebyshev case: no closed-form expression for f3
    params.dBripple = 0;
    params.f3_norm = find_f3_norm(a1, a2, a3, hval, n, k, w1, theta_i);
  end
  filt_params = find_filter_params(n, k, hval, w1, theta_i);
  params.qf = filt_params.qf;
  params.ff_norm = filt_params.ff_norm;
  params.fp_norm = filt_params.fp_norm;
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 qtb = calc_qtb_n(theta_a, theta_b, ql)
%calc_qtsb_n Compute the Butterworth QTS for the box Q given by QL
  qtb = 1 / (2 * (cos(theta_a) + cos(theta_b)) - 1 / ql);
end

function to_zero = func_to_solve_k(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 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 filt_params = find_filter_params(n, k, h, w1, theta_i)
  if isodd(n)
    filt_params.fp_norm = calc_fp_norm(k, h, w1);
  else
    filt_params.fp_norm = []; % No first-order filter
  end
  filt_count = size(theta_i, 2);
  if filt_count > 0
    filt_params.qf = zeros(1, filt_count);
    filt_params.ff_norm = zeros(1, filt_count);
    for i = 1:filt_count
      filt_params.qf(i) = calc_qf(k, theta_i(i));
      filt_params.ff_norm(i) = calc_ff_norm(k, h, w1, theta_i(i));
    end
  else
    filt_params.qf = [];
    filt_params.ff_norm = [];
  end
end

function fp_norm = calc_fp_norm(k, h, w1)
  fp_norm = sqrt(h) * w1 / k;
end

function qf = calc_qf(k, theta)
  qf = 0.5 * sqrt(1 + (tan(theta) / k)^2);
end

function ff_norm = calc_ff_norm(k, h, w1, theta)
  ff_norm = sqrt(h) * w1 / sqrt((k * cos(theta))^2 + sin(theta)^2);
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 filt = H(s, a1, a2, a3, h, n, k, w1, theta_i)
  filt = G(s, a1, a2, a3, h);
  if isodd(n)
    wpn = calc_fp_norm(k, h, w1);
    filt = filt * s / (s + wpn);
  end
  filt_count = size(theta_i, 2);
  for i = 1:filt_count
    qf = calc_qf(k, theta_i(i));
    wf = calc_ff_norm(k, h, w1, theta_i(i));
    filt = filt * hpf_second_order(s, qf, wf);
  end
end

function filt = hpf_second_order(s, qf, wf)
  sn = s / wf;
  filt = sn^2 / (sn^2 + 1/qf * sn + 1);
end

function resp = mag_squared_of_H(w, a1, a2, a3, h, n, k, w1, theta_i)
  cplx_H = H(j * w, a1, a2, a3, h, n, k, w1, theta_i);
  resp = real(cplx_H)^2 + imag(cplx_H)^2;
end

function to_zero = cutoff_freq_of_H(w, a1, a2, a3, h, n, k, w1, theta_i)
  to_zero = mag_squared_of_H(w, a1, a2, a3, h, n, k, w1, theta_i) - 0.5;
end

function w = find_f3_norm(a1, a2, a3, h, n, k, w1, theta_i)
  guess_for_w = sqrt(h);
  w = fzero(@(w)cutoff_freq_of_H(w, a1, a2, a3, h, n, k, w1, theta_i), guess_for_w, get_fzero_opts());
end