next Digital Waveguide Networks
previous Summary and Outline
up Circulant and Elliptic Feedback Delay Networks for Artificial Reverberation   Contents   Global Contents
global_index Global Index   Index   Search

Feedback Delay Networks

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

Figure: Order 3 Feedback Delay Network.
\includegraphics[scale=0.5]{eps/FDN.eps}

As depicted in Fig. Fig. 1, an FDN is built using N delay lines, each having a length in seconds given by \(\tau_{i}=m_{i}T\), where \(T=1/F_{s}\) is the sampling period. The complete FDN is given by the following relations:

$\displaystyle y(n)$ $\textstyle =$ $\displaystyle \sum_{i=1}^N c_{i}s_{i}(n) + dx(n)$ (1)
$\displaystyle s_{i}(n+m_i)$ $\textstyle =$ $\displaystyle \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
$\displaystyle Y(z)$ $\textstyle =$ $\displaystyle {\bf c}^T {\bf S}(z) + dX(z)$ (3)
$\displaystyle {\bf S}(z)$ $\textstyle =$ $\displaystyle {\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

$\displaystyle y(n)$ $\textstyle =$ $\displaystyle {\bf\gamma}^T {\bf w}(n) + dx(n)$  
$\displaystyle {\bf w}(n+1)$ $\textstyle =$ $\displaystyle {\bf A}^\dagger{\bf w}(n) + {\bf\beta}x(n)$ (5)

where
$\displaystyle {\bf\beta}^T$ $\textstyle =$ $\displaystyle [0, \ldots, 0, {\bf b}^T]$ (6)
$\displaystyle {\bf\gamma}^T$ $\textstyle =$ $\displaystyle [0, \ldots, 0, {\bf c}^T, \underbrace{0, \ldots, 0}_{N}]$ (7)

The state-transition matrix is

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

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


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

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

We have
$\displaystyle {\bf U}_j$ $\textstyle \in$ $\displaystyle {\cal C} ^{{(m_j - 2)} \times {(m_j - 2)}}$  
$\displaystyle {\bf R}_j$ $\textstyle \in$ $\displaystyle {\cal C} ^{{(m_j - 2)} \times N}$ (12)
$\displaystyle {\bf P}_j$ $\textstyle \in$ $\displaystyle {\cal C} ^{N \times {(m_j - 2)}}$  

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

\begin{displaymath}
N^{\dagger} = \sum_{i=1}^N m_i
\end{displaymath}

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

\begin{displaymath}
H(z) = {Y(z)\over X(z)} = {\bf c}^T [{\bf D}(z^{-1}) - {\bf A}]^{-1}{\bf b} + d . \\
\end{displaymath} (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

$\displaystyle {\left[{\bf I}- {\bf D}(z){\bf A}\right]^{-1}}$
  $\textstyle =$ $\displaystyle {\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
$\displaystyle {\Vert{\bf I}+ \cdots + {\bf D}^k(z){\bf A}^k + \cdots \Vert}$
  $\textstyle \leq$ $\displaystyle \Vert{{\bf I}}\Vert + \cdots + \Vert{\bf D}^k(z){\bf A}^k\Vert + \cdots$  
  $\textstyle \leq$ $\displaystyle 1 + \cdots + \Vert{\bf D}^k(z)\Vert\cdot\Vert{\bf A}^k\Vert + \cdots$  
  $\textstyle \leq$ $\displaystyle 1 + \cdots + \Vert{\bf D}(z)\Vert^k\cdot\Vert{\bf A}\Vert^k + \cdots$  
  $\textstyle \leq$ $\displaystyle 1 + \cdots + \Vert{\bf A}\Vert^k + \cdots$  
  $\textstyle =$ $\displaystyle \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

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

or
\begin{displaymath}
\det[z{\bf I}- {\bf A}^\dagger] = 0, \quad z \neq 0
\end{displaymath} (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

\begin{displaymath}
{\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)
\end{displaymath} (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

\begin{displaymath}
\det[{\bf A} - {\bf b}{\frac {1}{d}}{\bf c}^T- {\bf D}(z^{-1})] = 0
\end{displaymath} (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]

\begin{displaymath}
g_i = \alpha ^ {m_i}
\end{displaymath} (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:

\begin{displaymath}
G_i(z) = A^{m_i}(z),
\end{displaymath} (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.


next Digital Waveguide Networks
previous Summary and Outline
up Circulant and Elliptic Feedback Delay Networks for Artificial Reverberation   Contents   Global Contents
global_index Global Index   Index   Search

``Circulant and Elliptic Feedback Delay Networks for Artificial Reverberation'', by Davide Rocchesso and Julius O. Smith III, preprint of version in IEEE Transactions on Speech and Audio, vol. 5, no. 1, pp. 51-60, Jan. 1996.

Download PDF version (cfdn.pdf)
Download compressed PostScript version (cfdn.ps.gz)

(Browser settings for best viewing results)

Copyright © 2005-03-10 by Davide Rocchesso and Julius O. Smith III
Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
CCRMA  (automatic links disclaimer)