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