# Description of Matlab/Octave Design Software for Unassisted Alignments

Since the unassisted vented-box alignments are much more commonly used than the assisted ones, this section will describe in detail the implementation of some Matlab code for design with a given driver of only unassisted alignments using the techniques previously described herein. Then, in later sections, the design equations for the external filters of nth-order assisted alignments will be derived, and Matlab code for their design discussed. It isn't necessary to know Matlab programming to understand how the code works, but it's very helpful to know some programming language. Quirks specific to Matlab programming will be briefly explained as needed.

All the code provided here has also been verified to work under *GNU Octave*, a freeware Matlab workalike.

A browser-displayable version of the source code described below can be found at The Complete Matlab/Octave Source Code for Design of Unassisted Alignments. You can also download a zipped Matlab/Octave ".m" file of this code.

`calc_vented_box_params()`

The Matlab design software consists of a main function called `calc_vented_box_params()`

and a collection of helper functions. In Matlab, functions are contained in files whose file name is the name of the first function of the file with a ".m" extension appended. The file name in this case must be "calc_vented_box_params.m". The main function, always the first function of the file, is the only one callable by code defined in other .m files or typed from the Matlab command line. The function declaration `function params = calc_vented_box_params(qt, ql)`

specifies that a variable called `params`

is implicitly returned from the function. The function takes two arguments, `qt`

and `ql`

, which correspond to Q_{T} and Q_{L} respectively. It will be shown shortly that `params`

is a data structure that contains the values of α, h, f_{3} / f_{S}, and the ripple in dB (which is only non-zero for Chebyshev alignments).

After a comment line (which always begins with the "%" character), the first calculations done by `calc_vented_box_params()`

are of `theta_a`

and `theta_b`

. These correspond to θ_{a} and θ_{b}, the Butterworth pole angles introduced earlier. They have the values π/8 and 3π/8 respectively as per (43) for n = 4. Next, `qtb`

, which corresponds to Q_{TB}, the value of Q_{T} required for a Butterworth alignment, is computed. This computation corresponds to (19) and (20). Based on the value of Q_{T} passed in relative to the computed Q_{TB}, `calc_vented_box_params()`

passes control to either `find_system_params_cheby_butter()`

or `find_system_params_quasi_butter()`

, which do the actual "grunt work" of calculating α, h, f_{3} / f_{S}, and the ripple in dB. In this implementation, if Q_{T} < Q_{TB}, a quasi-Butterworth alignment is automatically chosen. The sub-Chebyshev alignment could also be chosen in this case depending on user preference, but the purpose of this simplified code is not to illustrate user interface issues, but to show how the equations previously derived are implemented in working code.

`find_system_params_cheby_butter()`

When Q_{T} ≥ Q_{TB}, the work of computing the design parameters is passed on to `find_system_params_cheby_butter()`

, which handles the Chebyshev and Butterworth alignments. Its purpose is to first compute k of (66) and (67) using a root-finding algorithm, then compute the remaining design parameters, all of which follow from k. It does so by first computing the eight constants `a`

through `h`

. These constants correspond to A through H of (76). Then Matlab's built-in root-finding function `fzero()`

is called to solve for k. The specific form of `fzero()`

used here takes three arguments. The first argument is a function of a single variable that returns a single scalar value. The second argument is an initial guess for the zero of the function. The third argument allows for configuring certain options of `fzero()`

. In this case, the function whose root we wish to find is `func_to_solve_cheby_butter()`

, and the first guess of the root is `guess_for_k`

. Since `func_to_solve_cheby_butter()`

has many arguments, it's necessary to make a new function having only one argument and to pass that function to `fzero()`

. This is done using the Matlab anonymous function feature, invoked by `@(k)`

, which creates a temporary, unnamed function of the lone variable k to pass to `fzero()`

. This anonymous, temporary function is just `func_to_solve_cheby_butter()`

with all arguments except k held constant. This technique prevents `a`

through `h`

, which never change, from being re-computed on each iteration of the root-finding algorithm. The relationship of `func_to_solve_cheby_butter()`

to the design equations previously derived herein will be described shortly. The `options`

structure passed to `fzero()`

is obtained by calling `get_fzero_opts()`

. The implementation of `get_fzero_opts()`

first gets the default options of `fzero()`

using `optimset()`

, then sets its `TolX`

component to the number `eps`

, which is about 2.2204*10^{-16}. The `TolX`

value controls the x-axis termination tolerance used by `fzero()`

. Since Matlab and Octave use different default `TolX`

values for `fzero()`

, the tolerance is explicitly set so the two programs get similar accuracy in the solution of k.

After k is found, `calc_tf_coeffs_cheby_butter()`

is called, which uses (81), (82), (83) and (80) to compute a_{1}, a_{2}, a_{3} and ω_{1} respectively. This collection of four values is returned as a structure from `calc_tf_coeffs_cheby_butter()`

and assigned to the `coeffs`

structure. The components of this structure are subsequently assigned to the corresponding individual variables `a1`

, `a2`

, `a3`

and `w1`

. Once the transfer function coefficients have been found, the square root of h is computed by `calc_sqrt_h()`

using (24), then α is computed by `calc_alpha()`

using (28).

After computing the key parameters h and α, the ripple in dB and f_{3} / f_{S} are all that remain to be found. To obtain the ripple in dB, the ripple factor ε is first calculated using (84) and assigned to the variable `epsi`

. Then the ripple in dB is computed using (85).

Finding f_{3} / f_{S} is done in several steps. First, μ is computed using (65) and assigned to the variable `mu`

. Then, ω_{3LP} is calculated using (69) and assigned to the variable `w3LP`

. Finally, f_{3} / f_{S} is computed by combining the results of (70) and (14).

The data structure `params`

returned by `find_system_params_cheby_butter()`

consists of the components `params.h`

, `params.alpha`

, `params.dBripple`

and `params.f3_norm`

. This data structure is returned to the calling function, `calc_vented_box_params()`

, which supplies the same data structure to code provided by the user. An example of a Matlab script that uses the code described herein is provided in the next section.

`func_to_solve_cheby_butter()`

For the Chebyshev and Butterworth cases, the final function to be described is `func_to_solve_cheby_butter()`

, whose zero corresponds to the value of k needed to establish a Chebyshev or Butterworth alignment. This function is passed indirectly to `fzero()`

, which holds all its arguments constant except k using the anonymous function approach described previously. So `func_to_solve_cheby_butter()`

is in effect a function of a single variable. It first calls `calc_tf_coeffs_cheby_butter()`

, described previously, to compute the transfer function coefficients. Then `calc_sqrt_h()`

, also described previously, is called to calculate the square root of h. The function to be set to zero is the left-hand side of (25). This is done in lieu of using (26), as it improves efficiency by eliminating the redundant evaluation of common subexpressions. Earlier attempts to use (22) and/or (23) were scrapped due to the singularity caused by having the square root of h in the denominator causing `fzero()`

to fail.

`find_system_params_quasi_butter()`

When Q_{T} is less than Q_{TB}, `find_system_params_quasi_butter()`

, is called to handle the quasi-Butterworth case. The techniques for finding the design parameters are described earlier in the quasi-Butterworth section. The solution depends on the fact that a_{1} and a_{3} can be expressed solely as functions of a_{2}. This leads to a procedure that finds the design parameters by solving a single nonlinear equation in the single unknown a_{2}, and finding all the other parameters from there. This equation is formally described in (40), but in practice a_{2} is found using a more efficient approach described at the end of the quasi-Butterworth section.

The procedure begins with a call to `fzero()`

, whose purpose is to compute a_{2}. This is done by passing in `func_to_solve_quasi_butter()`

as an argument and using the anonymous function technique described previously to treat `func_to_solve_quasi_butter()`

as a function of only a_{2}. The anonymous function technique holds Q_{T} and Q_{L} constant for each call to `func_to_solve_quasi_butter()`

such that the zero of the anonymous function is the value of a_{2} required for a quasi-Butterworth alignment. The second argument to `fzero()`

is the first guess for a_{2}. The a_{2} value for a Butterworth alignment is used here to aid convergence. The third argument to `fzero()`

is the `options`

structure described earlier.

After a_{2} is found, a_{1} and a_{3} are computed by `calc_a1_quasi_butter()`

and `calc_a3_quasi_butter()`

respectively. These functions use (38) and (39) respectively to obtain their results.

Computation of h and α is done in the same way as previously described for the Chebyshev and Butterworth cases. The ripple in dB is just set to zero, since the quasi-Butterworth alignments have no ripple. The ratio f_{3} / f_{s} is found using the `find_f3_norm()`

function.

The data structure `params`

returned by `find_system_params_quasi_butter()`

consists of the components `params.h`

, `params.alpha`

, `params.dBripple`

and `params.f3_norm`

. This data structure is returned to `calc_vented_box_params()`

, which passes it on to the user code. An example of such user code in the form of a Matlab script is provided in the next section.

`func_to_solve_quasi_butter()`

This is the function whose root gives the value of a_{2} consistent with a quasi-Butterworth alignment when its `qt`

and `ql`

arguments are held constant. It first computes a_{1} and a_{3} using (38) and (39) respectively. Then the square root of h is calculated using (24). The return value of `func_to_solve_quasi_butter()`

is the left-hand side of (25). This allows for finding a_{2} in a more efficient way than using (40) directly, as described in the quasi-Butterworth section.

`find_f3_norm()`

This function uses `fzero()`

to find the normalized frequency for which |G( jω)|^{2} = 0.5. Computation of |G(jω)|^{2} is done using `mag_squared_of_G()`

, which calls `G(s, a1, a2, a3, h)`

, passing in s = jω and computing the squared magnitude of the result. The `cutoff_freq_of_G()`

function just subtracts 0.5 from `mag_squared_of_G()`

, so it is the zero of this function that solves for the system -3dB frequency. Since this frequency must be normalized to f_{s}, G(s) is computed using `G(s, a1, a2, a3, h)`

by replacing T_{0} with 1 / sqrt(h) in (1). As a result, `find_f3_norm()`

returns the normalized frequency ω_{3} / ω_{s} = f_{3} / f_{s}. Software testing has shown that the technique of computing |G(jω)|^{2} directly and setting its value to 0.5 has better numerical stability than solving for d in (42).

## Matlab/Octave Script for Sample Vented-Box Design

A browser-displayable version of a script that exercises the source code described above can be found at Usage Example of Matlab/Octave Code for Unassisted Alignments. You can also download a zipped Matlab/Octave ".m" file of this script.

For those not familiar with Matlab or Octave, it should be mentioned that a statement ending in a semicolon is evaluated but its result is not displayed. A statement not ending in a semicolon is evaluated and its result displayed along with the variable to which the result is assigned.

The first line of the script is just a comment stating that the Thiele-Small parameters used are for the Dayton Audio RSS390HF subwoofer driver sold by Parts Express. The box Q is assumed to be 5. All parameters are from data supplied by Parts Express. Next, the values of any previously-defined variables are cleared, and the display format set to the "long" mode, done here mainly for testing purposes to track down potential errors that might show up as only small discrepancies from ideal values. Following this, the generator resistance `Rg`

, dc voice coil resistance `Re`

, and Q_{ES} value `qes`

are specified. Then Q_{E} is calculated and assigned to the `qe`

variable. In this specific case, Q_{E} has the same value as Q_{ES} because the generator resistance is specified as zero. Next, Q_{T} is calculated. Its value is the same as Q_{TS}, again because the generator resistance is zero. Following this, the free-air resonant frequency `fs`

is entered, along with the volume compliance V_{AS}, assigned to the variable `vas`

. It should be mentioned that V_{AS} is entered in units of cubic meters. V_{AS} is then converted to cubic feet. Next, the function `calc_vented_box_params()`

, described previously above, is called. It takes the `qt`

and `ql`

variables as arguments and returns a data structure, assigned here to the `params`

variable. This single function call finds all the needed design parameters. Once `params`

is found, the box tuning frequency `fb`

, the box volume in cubic feet `vb_cuft`

, the ripple in dB `dBripple`

, and the -3 dB frequency `f3`

are very straightforwardly calculated. Finally, the display format is set back to the default with the `format`

call.