Feedback Delay Networks

This section reviews FDNs along the lines indicated by Jot [9,8] with some modifications.

Figure: Order 3 Feedback Delay Network.

As depicted in Fig. Fig. 1, an FDN is built using N delay lines, each having a length in seconds given by , where is the sampling period. The complete FDN is given by the following relations:

$$y(n) = \sum_{i=1}^N c_{i}s_{i}(n) + dx(n)$$ (1)

$$s_{i}(n+m_i) = \sum_{j=1}^N a_{i,j}s_{j}(n) + b_{i}x(n)$$ (2)

where $s_{i}(n), \: 1 \leq i\leq N$, are the delay-line outputs at time sample $n$. If $m_i=1$ for each $i$, we obtain the conventional state-variable description of a discrete-time linear system [10]. In the present case, $m_i$ are normally large integers, so the variables $s_{i}(n)$ form only a small subset of the system state at time $n$, with the remaining state variables being the samples contained within the delay lines at time $n$. Using the $z$ transform, assuming zero initial conditions, we can rewrite (2) in the frequency domain as

$$Y(z) = {\bf c}^T {\bf S}(z) + dX(z)$$ (3)

$${\bf S}(z) = {\bf D}(z)\left[{\bf A S}(z) + {\bf b} X(z)\right]$$ (4)

where ${\bf s}^T(z)=[s_{1}(z), \ldots, s_{N}(z) ]$, ${\bf b}^T=[b_{1}, \ldots, b_{N}]$, and ${\bf c}^T = [ c_{1}, \ldots, c_{N}]$. The diagonal matrix ${\bf D}(z)~=~{\bf\mbox{diag}} \left( z^{-m_{1}}, z^{-m_{2}}, \dots z^{-m_{N}} \right)$ is called the ``delay matrix'' and ${\bf A} = [a_{i,j}]_{N \times N}$ is called the ``feedback matrix.''

The state variables of the FDN can be collected into a vector $\bf w$ as follows: List the variables contained in the first delay line from the $(m_1 - 1)th$ cell to the second cell, then those contained in the second delay line from the $(m_2 - 1)th$ cell to the second cell, and so on for the other delay lines; then attach the first cell of all the delay lines in increasing order, and finally the output variables $s_1$ to $s_N$.

By assuming that each delay line is longer than two samples, the state-variable description corresponding to this variable format for the system (2) can be found to be

$$y(n) = {\bf\gamma}^T {\bf w}(n) + dx(n)$$

$${\bf w}(n+1) = {\bf A}^\dagger{\bf w}(n) + {\bf\beta}x(n)$$ (5)

where

$${\bf\beta}^T = [0, \ldots, 0, {\bf b}^T]$$ (6)

$${\bf\gamma}^T = [0, \ldots, 0, {\bf c}^T, \underbrace{0, \ldots, 0}_{N}]$$ (7)

The state-transition matrix is

$${\bf A}^\dagger = \left[ \begin{array}{lllllll} {\bf U}_1... ...& {\bf0}& \dots & {\bf0}& {\bf A}& {\bf0} \end{array}\right]$$ (8)

where

$${\bf U}_j = \left[ \begin{array}{lllll} 0 & 1 & 0 & \dots ... ...dots & 1 \\ 0 & 0 & 0 & \dots & 0 \\ \end{array} \right],$$ (9)

$${\bf R}_j = \left[ \begin{array}{llllll} 0 & \dots & 0 & 0 ... ...{0 \dots 0}_{j-1} & 1 & 0 & \dots & 0 \\ \end{array}\right],$$ (10)

and

$${\bf P}_j = \left[ \stackrel {j-1 \left\{ \begin{array}{llll... ... 0 & 0 & \dots & 0 \\ [-.5in] \end{array} \right. } \right] .$$ (11)

We have

$${\bf U}_j \in {\cal C} ^{{(m_j - 2)} \times {(m_j - 2)}}$$

$${\bf R}_j \in {\cal C} ^{{(m_j - 2)} \times N}$$ (12)

$${\bf P}_j \in {\cal C} ^{N \times {(m_j - 2)}}$$

The order of the system (5) is equal to the sum of the delay-line lengths:

$$N^{\dagger} = \sum_{i=1}^N m_i$$

From (4), the transfer function is easily found to be

$$H(z) = {Y(z)\over X(z)} = {\bf c}^T [{\bf D}(z^{-1}) - {\bf A}]^{-1}{\bf b} + d . \\$$ (13)

Note that when $m_i=1$ for all $i$, the FDN specializes to a fully general state-space description [10]. This implies any linear, time-invariant, discrete-time system can be formulated as a special case of an FDN since every state-space description is a special case. This suggests that a wide variety of stable FDNs can be generated by starting with any stable LTI system whatsoever and performing the substitution $z^{-1} \leftarrow z^{-m_i}$ on each delay element, or any other conformal mapping which takes the unit circle to itself (another example being the Schroeder-allpass transformation $z^{-1} \leftarrow (\rho + z^{-m_i})/(1+\rho z^{-m_i})$). Stability is preserved even when the unit-sample delays of the original state space description are mapped using different conformal maps. This can be seen from the matrix power series expansion

$${\left[{\bf I}- {\bf D}(z){\bf A}\right]^{-1}}$$

$$= {\bf I}+ {\bf D}(z){\bf A} + \cdots + {\bf D}^k(z){\bf A}^k + \cdots$$ (14)

As long as $\Vert{\bf D}(e^{j\omega})\Vert\leq{\bf I}$, by making use of the triangle inequality and the Cauchy-Schwarz (or mutual consistency) inequality [20], we can write

$${\Vert{\bf I}+ \cdots + {\bf D}^k(z){\bf A}^k + \cdots \Vert}$$

$$\leq \Vert{{\bf I}}\Vert + \cdots + \Vert{\bf D}^k(z){\bf A}^k\Vert + \cdots$$

$$\leq 1 + \cdots + \Vert{\bf D}^k(z)\Vert\cdot\Vert{\bf A}^k\Vert + \cdots$$

$$\leq 1 + \cdots + \Vert{\bf D}(z)\Vert^k\cdot\Vert{\bf A}\Vert^k + \cdots$$

$$\leq 1 + \cdots + \Vert{\bf A}\Vert^k + \cdots$$

$$= \frac{1}{1-\Vert{\bf A}\Vert}$$ (15)

Therefore, as long as $\Vert A\Vert^n$ decays exponentially with $n$, stability is assured. The above derivation extends immediately to time-varying feedback matrices ${\bf A}_k$, provided $\Vert{\bf A}_k\Vert \leq 1-\epsilon$ for some worst-case $\epsilon > 0$.

The poles of the FDN are the solutions of either

$$\det[{\bf A} - {\bf D}(z^{-1})] = 0$$ (16)

or

$$\det[z{\bf I}- {\bf A}^\dagger] = 0, \quad z \neq 0$$ (17)

The matrix ${\bf A}^\dagger$ is not uniquely determined by ${\bf A}$. In fact, our ordering of the state variables differs from that used by Jot [8]. Our ordering gives

$${\bf A}^\dagger {{\bf A}^\dagger}^T =\mbox{diag}\left( {\... ...}_{m_1 - 1}, \dots {\bf I}_{m_N - 1}, {\bf A}{\bf A}^T \right)$$ (18)

From (18), we see that ${\bf A}$ is unitary if and only if ${\bf A}^\dagger$ is unitary. Since a unitary matrix has eigenvalues on the unit circle, from (17) we see that it is sufficient to choose a unitary matrix in order to have all the system poles on the unit circle. This yields a lossless FDN prototype.

By application of the matrix inversion lemma [10] to the transfer function (13), the system zeros are found as the solutions of

$$\det[{\bf A} - {\bf b}{\frac {1}{d}}{\bf c}^T- {\bf D}(z^{-1})] = 0$$ (19)

The formulation of (2) represents a prototype structure in the sense that, with the appropriate choice of feedback matrix, it is a lossless structure. In practice, we must insert attenuation coefficients and filters in the feedback loop. For example, one may insert a gain [9]

$$g_i = \alpha ^ {m_i}$$ (20)

at the output of each delay line in the FDN. This corresponds to replacing $D(z)$ with $D(z/\alpha)$ in (4). With this choice of the attenuation coefficients, all the system poles are uniformly contracted by a factor $\alpha$, thus ensuring a uniform decay of all the modes.

In a practical realization, we normally need to introduce frequency-dependent losses such that higher frequencies decay faster. We can do this by introducing lowpass filters after each delay line in place of the gains $g_i$. In this case, local uniformity of mode decay is still achieved by condition (20), where $g_i$ and $\alpha$ are made frequency dependent:

$$G_i(z) = A^{m_i}(z),$$ (21)

where $A(z)$ can be interpreted as the per-sample filtering [7,9,29].

Notice that uniform decay of all the modes, albeit arguably desirable in artificial reverberators for a smooth late time response, is not found in actual rooms. Normal modes are associated with standing waves which have an absorption that depends on their orientation. For example, in a rectangular enclosure, waves traveling in a direction normal to a wall are less absorbed than oblique waves [17, p. 392], so that the corresponding standing waves (expressible as the superposition of traveling waves in opposite directions) have different reverberation times. The room-acoustics interpretation of FDNs provided in Section V points to ways of modeling such uneven decays.