The Mass-Spring-Damper Solution

• As previously indicated, the flow through the reed channel is approximated “quasi-statically” using the Bernoulli equation and given by
    $$u = w \, (x+H) \left(\frac{2 \, \vert p_{\Delta} \vert}{\rho} \right)^{1/2} \mbox{sgn($p_{\Delta}$)},$$ (31)
  where $w$ is the reed channel width, $x$ is the time-varying reed position, calculated from Eq. (29), and $H$ is the equilibrium tip opening.

• For single-reed geometries, the pressure and flow in the reed channel can be approximated as equivalent to the pressure and flow at the entrance to the instrument air column. This approximation is based on continuity and detachment of volume flow at the end of the reed channel such that pressure is not recovered in the mouthpiece.

• Thus, the acoustic interaction at the interface of the reed and air column can be solved using Eqs. (29) and (31), together with a description of the input impedance of the attached air column.

• For a cylindrical pipe, we can compute the impedance using the digital waveguide structure of Fig. 16.

  Figure 16: A digital waveguide cylindrical pipe impedance model (from Scavone (1997)).

  

• From Fig. 16, it is clear that
    $$p_{0} = 2 p_{0}^{-} + Z_{c} u_{0},$$ (32)
  where $p_{0}^{-}$ is the traveling-wave pressure entering the reed junction from the downstream air column.

• Because of mutual dependencies, however, an explicit solution of these equations can be problematic. In a discrete-time computational context, these mutual dependencies can be understood to result in delay-free loops.

• In Guillemain et al. (2005), the reed system is discretized using a centered finite difference approximation that avoids a direct feedforward path through the reed transfer function. The resulting system equations can then be expressed in terms of a second-order polynomial equation and an explicit solution found.

• The centered finite-difference approximation of Eq. (29) results in a digital filter structure of the form
    $$\frac{X(z)}{P_{\Delta}(z)} = \frac{-z^{-1}/\mu_{r}}{(f_{s}^2 + \f... ...(\omega_{r}^2 - 2 f_{s}^2) z^{-1} + (f_{s}^2 - \frac{g_{r} f_{s}}{2}) z^{-2}},$$ (33)
  where $f_{s}$ is the computational sample rate.

• As noted in Guillemain (2004), however, this filter structure is only stable for $\omega_{r} < f_{s}\sqrt{4 - (g_{r}/\omega_{r})^2}$, limiting its use at low sample rates and/or with high reed resonance frequencies.

• A direct application of the bilinear transform to the system of Eq. (29) results in a digital filter structure given by
    $$\frac{X(z)}{P_{\Delta}(z)} = \frac{-1/\mu_{r} \left[ 1 + 2 z^{-1}... ...ga_{r}^2 - \alpha^2) z^{-1} + (\alpha^2 - g_{r}\alpha + \omega_{r}^2) z^{-2}},$$ (34)
  where $a_{0} = \alpha^2 + g_{r}\alpha + \omega_{r}^2$ and $\alpha$ is the bilinear transform constant that controls frequency warping.

• Note that we can achieve an exact continuous- to discrete-time frequency match at the resonance frequency of the reed by setting $\alpha = \omega_{r} / \tan( \omega_{r} / 2 f_s )$.

• In this case, the use of the bilinear transform guarantees a stable digital filter at any sample rate. The presence of the direct feedforward path in Eq. (34), however, prohibits the explicit reed interface solution mentioned above.

• We therefore seek an alternative form of Eq. (34) that preserves stability and avoids an undelayed feedforward coefficient in the transfer function numerator.

• By default, the bilinear transform substitution produces a system with “zeroes” at $z = \pm 1$ (or at frequencies of 0 and $f_{s}/2$ Hz). While this result is often desirable for digital resonators, we can modify the numerator terms without affecting the essential behavior and stability of the resonator.

• In fact, it is the numerator terms that control the phase offset of the decaying oscillation. Thus, we can modify and renormalize the numerator to produce a filter structure of the form
    $$\frac{X(z)}{P_{\Delta}(z)} = \frac{-4 z^{-1} /\mu_{r}}{a_{0} + 2 (\omega_{r}^2 - \alpha^2) z^{-1} + (\alpha^2 - g\alpha + \omega_{r}^2) z^{-2}}.$$ (35)

• The frequency- and time-domain responses of the centered finite-difference and “modified” bilinear transform filter structures are shown in Fig. 17 for a reed resonance frequency $f_{r} = 2500$ Hz and $f_{s} = 22050$ Hz.

  Figure 17: Reed filter frequency and impulse responses for centered finite-difference and modified bilinear transform structures with $f_{r} = 2500$ Hz and $f_{s} = 22050$ Hz: magnitude frequency response (top) and impulse response (bottom).

  

• The complete clarinet model involves the calculation of the reed displacement using this stable reed model, the volume flow through the reed channel as given by Eq. (31), and the relationship between flow and pressure at the entrance to the air column as given by Eq. (32).

• Because the reed displacement given by Eq. (35) does not have an immediate dependence on $p_{\Delta}$, it is possible to explicitly solve Eqs. (32) and (21), as noted in Guillemain et al. (2005), by an expression of the form
    $$u_{0} = 0.5 \, \left( B \sqrt{(Z_{c}B)^2 + 4 A} - Z_{c} B^2 \right) \mbox{sgn($A$)},$$ (36)
  where $A = p_{m} - 2 p_{0}^{-}$ and $B = w \, (x+H) ( 2 / \rho )^{1/2}$ can be determined at the beginning of each iteration from constant and past known values.

• Whenever the reed channel height $x+H < 0$, $u_{0}$ is set to zero and $p_{0} = 2 p_{0}^{-}$.

• In Fig. 18, the normalized pressure response of the complete DW synthesis model is plotted using both reed models with $f_{r} = 2500$ Hz and $f_{s} = 22050$ Hz.

• The behaviors are indistinguishable for these system parameters, though as indicated above it is possible to run the modified bilinear transform model at significantly lower sample rates (and with higher reed resonance frequencies).

  Figure 18: Normalized pressure response from complete DW synthesis model using the centered finite-difference and modified bilinear transform structures with $f_{r} = 2500$ Hz and $f_{s} = 22050$ Hz.