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₃ / 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₃ / 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₁, a₂, a₃ and ω₁ 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₃ / 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₃ / 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₃ / 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₁ and a₃ can be expressed solely as functions of a₂. This leads to a procedure that finds the design parameters by solving a single nonlinear equation in the single unknown a₂, and finding all the other parameters from there. This equation is formally described in (40), but in practice a₂ 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₂. 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₂. 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₂ required for a quasi-Butterworth alignment. The second argument to fzero() is the first guess for a₂. The a₂ value for a Butterworth alignment is used here to aid convergence. The third argument to fzero() is the options structure described earlier.

After a₂ is found, a₁ and a₃ 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₃ / 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₂ consistent with a quasi-Butterworth alignment when its qt and ql arguments are held constant. It first computes a₁ and a₃ 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₂ 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ω)|² = 0.5. Computation of |G(jω)|² 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₀ with 1 / sqrt(h) in (1). As a result, find_f3_norm() returns the normalized frequency ω₃ / ω_s = f₃ / f_s. Software testing has shown that the technique of computing |G(jω)|² 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.