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.
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,
ql, which correspond to QT and QL respectively. It will be shown shortly that
params is a data structure that contains the values of α, h, f3 / fS, 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_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 QTB, the value of QT required for a Butterworth alignment, is computed. This computation corresponds to (19) and (20). Based on the value of QT passed in relative to the computed QTB,
calc_vented_box_params() passes control to either
find_system_params_quasi_butter(), which do the actual "grunt work" of calculating α, h, f3 / fS, and the ripple in dB. In this implementation, if QT < QTB, 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.
When QT ≥ QTB, 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
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
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
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
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 a1, a2, a3 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
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 f3 / fS 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 f3 / fS 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, f3 / fS 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.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.
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.
When QT is less than QTB,
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 a1 and a3 can be expressed solely as functions of a2. This leads to a procedure that finds the design parameters by solving a single nonlinear equation in the single unknown a2, and finding all the other parameters from there. This equation is formally described in (40), but in practice a2 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 a2. 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 a2. The anonymous function technique holds QT and QL constant for each call to
func_to_solve_quasi_butter() such that the zero of the anonymous function is the value of a2 required for a quasi-Butterworth alignment. The second argument to
fzero() is the first guess for a2. The a2 value for a Butterworth alignment is used here to aid convergence. The third argument to
fzero() is the
options structure described earlier.
After a2 is found, a1 and a3 are computed by
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 f3 / fs is found using the
The data structure
params returned by
find_system_params_quasi_butter() consists of the components
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.
This is the function whose root gives the value of a2 consistent with a quasi-Butterworth alignment when its
ql arguments are held constant. It first computes a1 and a3 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 a2 in a more efficient way than using (40) directly, as described in the quasi-Butterworth section.
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 fs, G(s) is computed using
G(s, a1, a2, a3, h) by replacing T0 with 1 / sqrt(h) in (1). As a result,
find_f3_norm() returns the normalized frequency ω3 / ωs = f3 / fs. 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 QES value
qes are specified. Then QE is calculated and assigned to the
qe variable. In this specific case, QE has the same value as QES because the generator resistance is specified as zero. Next, QT is calculated. Its value is the same as QTS, again because the generator resistance is zero. Following this, the free-air resonant frequency
fs is entered, along with the volume compliance VAS, assigned to the variable
vas. It should be mentioned that VAS is entered in units of cubic meters. VAS is then converted to cubic feet. Next, the function
calc_vented_box_params(), described previously above, is called. It takes the
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