This chapter will present an analysis of feedback delay networks.  Section 3.1 reviews the properties of the unit comb filter introduced in section 2.5.  Section 3.2 covers the properties of the parallel comb filter network and explains how to achieve a frequency dependent reverberation time with this design.  Section 3.3 derives a general feedback delay network based on the parallel comb filter network.

3.1         Unit Comb Filter Pole Study

The transfer function of the unit comb filter of Fig. 2.5 was given by ( 2.10 ) and can be rearranged as:
 

( 3.1 )

In the equation above, the poles z_k (where 0 km-1) are defined by , where  and .  As we already saw in Fig. 2.6c, the poles have equal magnitudes and are located around the unit circle.  If the feedback gain is less than unity, this pole pattern corresponds to exponentially decaying sinusoids of equal weight, as illustrated in Fig. 2.6d.  The inverse transform of ( 3.1 ) gives the following filter impulse response:
 

(for n  0)

( 3.2 )

3.2         Parallel Comb Filter Analysis

When combining N unit comb filters in parallel, the transfer function of the system becomes:
 

( 3.3 )

By choosing “incommensurate” delay lengths m_p, all the eigenmodes (frequency peaks) of the unit comb filters are distinct (except at  = 0 and = ).  The total number of resonating frequency peaks is equal to the half the sum of all unit comb filters’ delay lengths expressed in samples.

3.2.1        Equal Magnitude of the Poles

As Jot and Chaigne [16] pointed out, if the magnitudes of the poles are not the same, some resonating frequencies will stand out.  This is because the decay time of each comb filter will be different.  To make every comb filter decay at the same rate, the magnitude of all the poles must be equal.  Therefore, we have to follow the rule:
 

for any p

( 3.4 )

This relates the length m_p of every comb filter p with its feedback gain g_p.  Fig. 3.1 shows the smooth decay of a parallel comb filter’s impulse response (with equal pole magnitude).

Fig. 3.1.  Impulse response of a parallel comb filter with equal pole magnitude (with f_s = 1 kHz, m₁ = 8 in red, m₂ = 11 in green, and m₃ = 14 in blue).

Respecting condition ( 3.4 ) ensures that all the comb filters will have an equal decay time.  However, comb filters that have longer delay lengths  m_p will produce frequency peaks (also called eigenmodes) with weaker gains, as we can see in ( 3.3 ).  To reduce this effect, we can try to keep the range from the shorter delay length to the longer delay length as small as possible.  That is why Schroeder [39] suggested a maximum spread of 1:1.5 between the comb filter’s delay lengths. The other way to get around this problem is to weight the gains of each filter proportionally to their delay lengths before summing their outputs.  This ensures that every comb filter’s response will be heard equally.  The resulting frequency response is shown in Fig. 3.2.

Fig. 3.2.  Frequency response of the parallel comb filter of Fig. 3.1 with weighted gains.

	Frequency density:		( 3.5 )
	Time density:		( 3.6 )

where  is the average delay length of a comb (in seconds) and N is the number of comb filters.  Thus, given a frequency density and a time density, we can find how many combs should be used and what their average delay lengths would be using the following equations:
 

and

( 3.7 )

       For example, to achieve a time density D_t = 1000 (as suggested by Schroeder [39]) and a frequency density D_f  = 0.45, we would need 21 comb filters and an average delay of 21 ms.  Griesinger suggested that high-bandwidth reverberation units should have a higher time density, that could even exceed 10,000 echoes per second.  This would require more than 67 comb filters!  Luckily, feedback delay networks will take care of this issue as we will see in section 3.3.

3.2.3        Reverberation Time

In chapter 2, the reverberation time was defined as the time for the decay to decrease 60 dB.  For a parallel comb filter with equal magnitude, all comb units are decaying at the same rate.  The following equation gives their decay rate in dB per sample period :
 

for any p

( 3.8 )

This equation is equivalent to ( 3.4 ) but here,  and G_p represent  and g_p in dB ( = 20 log() and G_p = 20 log(g_p)).  Using Schroeder’s definition of reverberation time and ( 3.8 ), we can find the reverberation time for each comb filter:
 

( 3.9 )

where T_r is the reverberation time, T is the sampling period and _p = m_p T.  Modifying either the comb filters’ delay lengths or their feedback gains can change their reverberation time.  The physical analogy to an increase of the delay lengths would be an increase in the virtual room dimensions.  On the other hand, the feedback gain attenuation could be thought of as the amount of absorption occurring during sound propagation in the air.  However, we must be careful when changing these gains since they also account for a change in the amount of wall absorption (that should not change when attempting to modify only the room size).  Multiplying the delay lengths by  or dividing the feedback gains (in dB) by  would cause the reverberation to be increased by .  It should be noted that in both cases, the time and frequency density of the system would be modified.

3.2.4        Achieving Frequency Dependent Reverberation Time

As suggested by Schroeder [39] and implemented later by Moorer [27], inserting a low-pass filter in the feedback loop of a comb filter introduces some frequency dependent decay rates.  This can be used to simulate wall absorption, since a wall tends to absorb more high frequencies than low frequencies.

By replacing the feedback gain g_p with a filter having a transfer function h_p(z), we can achieve a frequency dependent reverberation time.  Care must be taken when inserting these absorbent filters to ensure that no unnatural coloration is introduced in the reverberation.  We would like every comb filter to decay with the same relative frequency absorption.  This can be achieved by respecting a condition that Jot and Chaigne [16] called continuity of the pole locus in the z-plane:

[…] in any sufficiently narrow frequency band (where the reverberation time can be considered constant) all eigenmodes must have the same decay time.  Equivalently, system poles corresponding to neighboring eigenfrequencies must have the same magnitude.

h_p() |

for any p

( 3.10 )

The reverberation time is also frequency dependent and is given by:
 

0

( 3.11 )

Fig. 3.3 shows a parallel comb filter with frequency dependent decay time (T_r(0) = 3 seconds and T_r() = 0.15 s) respecting the continuity of the pole locus condition.

Finally, we must note that changing the conjugate pole locations of some, or all, of the comb filters in the unit circle will cause some of the eigenmodes to have more energy than others.  This changes the general frequency response of the output, as we can see in Fig. 3.3b.  However, with the addition of a simple tone correction filter (in series with the parallel comb filter) having a transfer function t() as:
 

( 3.12 )

an overall frequency response that is independent of the frequency-specific reverberation time can be achieved, as shown in Fig. 3.4.

Fig. 3.3.  Parallel comb filter with frequency dependent decay time (with f_s = 1 kHz, m₁ = 8 in red, m₂ = 11 in green, and m₃ = 14 in blue).

Fig. 3.4.  Frequency response of the filter network of Fig. 3.3 with the addition of a tone corrector.

where

( 3.13 )

where  and .   The corresponding tone corrector he used was an FIR filter with the following transfer function:
 

( 3.14 )

where .  The resulting signal flow using both filters is shown in Fig. 3.5.

Fig. 3.5.  Flow diagram using first order IIR absorbent filters and one first order FIR tone correction filter.

3.3         General Feedback Delay Network

As we saw in the previous section, a considerable number of comb filters are required to achieve a good time density in a parallel comb filter design, given a reasonable frequency density.  Some other filter designs achieving a greater time density have been mentioned in the second chapter.  Some of those designs, such as the series all-pass filter, produce a build-up of echo density with time as can be observed in real rooms.  Ideally, we would like to find a single network that would be general enough to take advantage of the properties of both parallel comb and series all-pass filters.  We would also like to achieve the same frequency dependent characteristics we were able to model with the parallel comb filter.  This section describes the steps that were taken to achieve this goal.

Gerzon [11] was the first to generalize the notion of a unitary multi-channel network, which is basically an N-dimensional equivalent of the unit all-pass filter.  Stautner and Puckette [47] then proposed a general network based on N delay lines and a feedback matrix A, as shown in Fig. 3.6.  In this matrix, each coefficient a_mn corresponds to the amount of signal coming out of delay line n sent to the input of delay line m.  The stability of this system relies on the feedback matrix A.  The authors found that stability is ensured if A is the product of a unitary matrix and a gain coefficient g, where |g| < 1.  Another property of this system is that the output will be mutually incoherent and thus can be used without additional processing in multi-channel systems.

Fig. 3.6.  Stautner and Puckette’s four channel feedback delay network.

Jot and Chaigne further generalized this design by using the system shown in Fig. 3.7.  Using the vector notation and the z-transform, the equations corresponding to this system are:
 

		( 3.15 )
		( 3.16 )

where:
 

		( 3.17 )
		( 3.18 )

For multiple-input and multiple-output systems, the vectors b and c becomes matrices.  Eliminating s(z) in equation ( 3.17 ) and ( 3.18 ) gives the following system transfer function:
 

( 3.19 )

The system’s zeros are given by:
 

( 3.20 )

Thus, the poles of this system are the solution of the characteristic equation:
 

( 3.21 )

The analytical solution of the equation is not trivial.  However, for specific feedback matrices, this equation is easy to solve.

Jot and Chaigne pointed out that we could represent any combination of unit filters (either comb or all-pass) by using the appropriate matrix.  For example, if A is a diagonal matrix, the system represents the parallel comb filter described in section 3.2.  Triangular matrices have also been studied because in that case, equation ( 3.21 ) reduces to:
 

( 3.22 )

The series all-pass filter is itself a network with a triangular feedback matrix (having diagonal elements equal to the feedback gains g_p).  Thus, a series all-pass filter that has the same delay lengths and feedback gains as a parallel comb filter also has the same eigenmodes (resonant frequencies and decay rates) as the parallel comb filter.

Fig. 3.7.  Jot’s general feedback delay network.

3.3.2        Achieving Frequency Dependent Reverberation Time

To control the decay time of a network, Jot and Chaigne proposed the use of a lossless feedback matrix combined with absorption filters right after each delay line.  To make all the comb filters decay at the same rate and to respect the continuity of the pole locus, the following equation has to be respected:
 

( 3.23 )

This will move the poles from the unit circle to a curve defined by the reverberation time T_r.  Again, this will introduce an overall high frequency attenuation that can be rectified by the use of a general tone correction filter described in section 3.2.4.  The resulting FDN including absorption and tone correction filter is shown in Fig. 3.8.

Fig. 3.8.  Jot’s general feedback delay network with absorption and tone correction filters.

The choice of the feedback matrix is critical to the overall sound and computational requirement of the reverberator.  As we discussed in section 3.3.1, many matrices will create a lossless prototype, such as the diagonal matrix equivalent to the parallel comb filter network.  For example, Stautner and Puckette [47] used the matrix:
 

( 3.24 )

where A is lossless when |g| = 1, because of its interesting properties when used in a multi-channel output system.  However, significant improvement can be obtained by using a unitary matrix with no null coefficients.  This will produce a maximum density while requiring the same amount of delay lines.

One such matrix is taken from the class of Householder matrices [17] and can be implemented efficiently as we will see:
 

( 3.25 )

where J_N is an N x N permutation matrix, and u_N is an N x 1 column vector of ones.  The resulting matrix contains two different non-zero values, thus maximizing the echo density.  Also, A _N · x(z) can be computed efficiently by permuting the elements of x(z) according to J_N, and adding to these the sum of the elements of the input times the factor .  This requires around 2N operations, compared to the N² operations usually required for an N x N matrix multiplication.  For the special case where J_N is the identity matrix I_N, the resulting system behaves like a parallel comb filter with a maximum echo density as shown in Fig. 3.9.  Jot noted that this system produces periodic clicks at a time interval equivalent to the sum of all delay lengths.  However, these clicks can be avoided in the output of the system by inverting the sign of every other coefficient c_N.

Choosing J_N to be a circular permutation matrix is another interesting possibility, because in this configuration, the delay lines feed each other serially.  This simplifies memory management in the final implementation of the system.

Fig. 3.9.  Parallel comb filter modified to maximize the echo density.

( 3.26 )

We can use the Discrete Fourier Transform (DFT) matrix T to diagonalize the circulant matrices:
 

( 3.27 )

This implies that the eigenvalues of A can be computed by finding the DFT of the first row:
 

( 3.28 )

where {(A)} are the eigenvalues of A (the diagonal elements of D).  The circulant matrix will be lossless if these eigenvalues have unit modulus.  There are two main advantages in the use of circulant feedback matrices: a) specific eigenvalues can be specified during the design of the system, and b) the matrix multiplication only takes about N log(N) operations to compute since the Fast Fourier Transform can be used.

3.4         Comparing FDN with Waveguide Networks

Fig. 3.10 represents an N-branch Digital Waveguide Network (DWN) with a structure equivalent to an  N^th order FDN.  Each branch of the structure represents an FDN delay line and the length of each waveguide is half the length of the equivalent delay line.  This is because a traveling wave must travel from the center junction to the termination of the branch and then back to the center junction.  Odd delay line lengths can be converted to an even length waveguide with a reflecting termination made of a single delay element.
 

Fig. 3.10.  Waveguide network consisting of a single scattering junction to which N branches are connected.  Each branch is terminated with a perfect non-inverting reflection, indicated by a black dot [46].

( 3.29 )

where p⁺ is the vector of the incoming traveling wave samples arriving at the junction (the delay-line outputs in an FDN) at time n, p¯ is the vector of outgoing traveling wave samples (the delay-line inputs in an FDN) leaving the junction at time n, and A is the scattering matrix (the feedback matrix in an FDN) associated with the waveguide junction.

Rocchesso and Smith [31] have shown that the scattering matrix is said to be lossless if the total complex power at the junction is scattering invariant, i.e.,
 

		( 3.30 )
		( 3.31 )

where  is a Hermitian, positive-definite matrix which can be interpreted as a generalized junction admittance.  In the case of Fig. 3.10, each waveguide has a characteristic admittance _i and the matrix  = diag(_1,2, …, _N).  In the case where A is unitary,  = I.  Thus, a DWN where each characteristic admittance ?i is unity is equivalent to an FDN with a unitary feedback matrix.  We can see that the DWN leads to a more general class of lossless scattering matrices, because each waveguide may have a characteristic admittance different than unity.  However, not all lossless scattering matrices can be interpreted as a physical junction of N waveguides (for example, a permutation matrix).