## 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)
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
``````