Several types of time varying modifications can be made to a general FDN algorithm to enhance the sound of its reverberation or to add life to an otherwise static reverberation.  For example, changing the feedback matrix of an FDN in real-time would create a reverberator whose timbre changes with time.  One could also change the output matrix (coefficients c_N in Fig. 3.7).  This thesis investigates the use of modulation to change the FDN delay lengths in real time.  Changing these delay lengths causes the resonant frequencies to change constantly, thus achieving a reverberation with a flatter perceived frequency response.

Fig. 4.1 and Fig. 4.2 show a typical modulated delay line where the output position of the delay line moves constantly around a center position.  The output of the delay line is given by:
 

		( 4.1 )
	where	( 4.2 )

Here, NOMINAL_DELAY is the center position of the modulation, MOD_DEPTH is the excursion of the modulation on each side of the center position, y_mod(n) is the output of the waveform generator, and D is a positive real number that can be split into an integer and a fractional part:
 

( 4.3 )

When d  0, the value of x(n - D) falls between two samples and thus must be interpolated.

Fig. 4.1.  A modulated delay line.

Fig. 4.2.  Detailed version of the delay line of Fig. 4.1.

4.1         Modulation Types

4.1.1        Sinusoid Oscillator

The most common modulation source is a simple sinusoid oscillator.  There are many ways to implement such oscillators.  One of them uses look-up tables where coefficients are already computed and stored in memory.  However, since the purpose of this thesis is to produce an algorithm that requires low memory requirements (in addition to being efficient), we chose another way to do it requiring only a few words of memory.  The implementation is based on the z-transform of the sin function:
 

	<=>	( 4.4 )
		( 4.5 )

where , f_osc is the desired oscillator frequency, and f_s is the sampling rate.  X(w) is an impulse signal used to start the oscillator, defined in the time domain by:
 

x(n) = 1,  for n = 0 x(n) = 0,  for n 0

( 4.6 )

Taking the inverse z-transform of ( 4.5 ), the following time domain equations are obtained:
 

			( 4.7 )
		= 0  for n > 1
	<=>
	for n > 1		( 4.8 )

where the initial conditions are given by:
 

		( 4.9 )
		( 4.10 )

where f_osc is the desired oscillator frequency and f_s is the sampling rate.  It is important to note that this algorithm can be unstable for some frequencies, due to the finite word length of the implementation platform.  More details on the implementation of this algorithm will be presented in section 5.4.1.

4.1.2        Filtered Pseudo-Random Number Sequence

       We can also use a low-pass filtered pseudo-random number sequence as a modulation source.  In practice, a simple uniform pseudo-random number generator (RNG) is sufficient and can be implemented efficiently.  One class of RNG proposed by Knuth [23] is based on the linear congruence method, which uses the following equation:
 

( 4.11 )

The initial value or seed is generally not important, unless multiple RNG are used at the same time.  In this case, it is desirable to choose a unique seed for each of them.  The choice of parameters a and c is crucial since a good generator should generate m values before repeating the output sequence, and every number between 0 and m-1 should be output once in every period of the generator.  However, if multiple generators have to be used at the same time, a different seed should be chosen for each RNG so that all the sequences appear to be different from each other.  The modulus m is usually chosen according to the word length of the implementation platform, to execute the mod function “transparently.”  For example, we chose a modulus of 32 because our implementation platform uses a 32-bit word length processor.  More implementation details will be given in section 5.4.2.  The other parameters we used (a = 1664525 and c = 32767) were chosen according to the rules of Knuth.

Once we get the raw pseudo-random number sequence, we must process it to obtain smooth oscillations.  We can achieve slow oscillations using one pseudo-random number for every f_s/N sample, where f_s is the sampling rate and N is a constant to be determined experimentally.  We could then obtain the intermediate values using linear interpolation.  However, the abrupt changes in oscillation slopes would create sudden pitch changes in the reverberation output making this method unsuitable for our application.  A better sounding (although more computationally intensive) approach is to pad zeros in the intermediate points and to low-pass filter the resulting sequence.  Zero padding the original sequence makes the entire process more efficient since we compute one random number every N samples (we found that N = 100 gave good results).  We tried several low-pass filters and found that a fourth order Chebyshev filter (obtained using the bilinear transform), with a cut-off frequency of 2 Hz and a pass-band ripple of 0.5 dB was required to obtain a smooth output with no specific prominent oscillation frequency.  The resulting signal with a seed of 2 is shown in Fig. 4.3.

Fig. 4.3.  Filtered pseudo-random number sequence.

Triangular waveform generators do not give good results when used as modulators, as they will produce reverberation tails consisting of two distinct pitches every period.  The first pitch will correspond to the rising ramp and the second pitch will correspond to the falling ramp.  To avoid this sudden change in pitch, we could low-pass filter the waveform.  However, this would require an additional amount of computations.  The resulting waveform would then be similar to a sinusoidal waveform, but would require more computations.  Since our goal was to reduce computations to a minimum, we did not use filtered triangular waveforms in our implementation.

Another common waveform is the sawtooth.  Even if sawtooth modulators produce a reverberation tail having a constant pitch, the tail will suffer from clicks every generator period caused by the abrupt transition in the waveform.  Again, appropriate low-pass filtering could remove these clicks, but it would result in an unacceptable increase in computation time.

4.2         Interpolation

In a time-varying implementation, a modulator is constantly changing the length of a specific delay line.  For example, from a minimum delay of 500 to a maximum delay of 508 samples.  Thus most of the intermediate delay values will not be exact integers (for example 501.25), so a delay of a fraction of a sample will be requested.  Since the modulator rate will be slow (no more than a few hertz) and only a few samples wide, the delay length can not simply be rounded to the nearest output sample in the delay line.  It is important to use interpolation, as rounding will result in audible noise or clicks in the reverberation tail.

4.2.1        General Concepts

The z-domain transfer function of the ideal delay of Fig. 4.1 is:
 

( 4.12 )

The frequency response of this ideal system is given by:
 

( 4.13 )

where  = 2ft is the normalized frequency and T is the sampling period.  Thus the magnitude and phase response of H() are respectively:
 

	for all	( 4.14 )
		( 4.15 )

Finally, we find the group delay by taking the derivative of the phase response and inverting its sign:
 

( 4.16 )

The implementation of a fractional delay system to interpolate values from a band-limited signal can be considered as an approximation of the ideal linear-phase all-pass filter with a unity magnitude and a constant group delay D.

The impulse response of such a filter is computed by taking the inverse Fourier transform:
 

for all n

( 4.17 )

Using H_id() found in ( 4.13 ), the ideal impulse response of the all-pass filter can be found:
 

for all n

( 4.18 )

For integer values of D, the impulse response is a single impulse centered at n = D.  For all other values of D, the impulse response corresponds to a shifted equivalent of the sinc function centered at D.  Such a filter is impossible to realize in a real-time system since it is infinitely long and also non-causal.

Fig. 4.4.  Impulse response of an ideal all-pass filter with delay a) D = 3 and b) D = 3.3 (from [25]).

4.2.2        Lagrange Interpolation (FIR)

The FIR methods all have the same z-domain transfer function:
 

( 4.19 )

Coefficients are determined such as the norm of the frequency-domain error function
 

( 4.20 )

is minimized.

The idea behind Lagrange interpolation is to make this error function maximally flat at a certain frequency, typically at _o = 0, so that the approximation is at its best close to this frequency.  This can be achieved by making the derivatives of the frequency-domain error function equal to zero at this point:
 

for n = 0, 1, 2, … N

( 4.21 )

Differentiating and substituting _o = 0 in equation ( 4.21 ), a set of L = N + 1 linear equations are obtained.  They can be expressed as an impulse response of the form:
 

for n = 0, 1, 2, … N

( 4.22 )

or, in the matrix notation form:
 

( 4.23 )

where h is a vector containing the coefficients, V is an L x L Vandermonde matrix:
 

( 4.24 )

and v is:
 

( 4.25 )

The solution to ( 4.21 ) is equal to the classical Lagrange interpolation formula and the coefficients are obtained by making the interpolation polynomial pass through a given set of data values.  This leads to the solution:
 

for n = 0, 1, 2, … N

( 4.26 )

If N equals 1, this corresponds to linear interpolation between two samples.  In this specific case, h(0) = 1 - D and h(1) = D.  The magnitude and phase delay responses of low-order Lagrange interpolation filters are shown in Fig. 4.5 and Fig. 4.6.

Fig. 4.5.  Lagrange interpolating filters of length L = 2, 3, 4, 5, and 6 with d = 0.5 (from [25]).

Fig. 4.6.  Lagrange interpolating filter of length L = 4, with d = 0 to 0.5 (from [25]).

As we can see, the even length filters (L = 2, 4, and 6) are linear phase.  However, their magnitude responses fall to zero at  = .  On the other hand, the odd length filters (L = 3 and 5) have a better magnitude response, but their phase delays are worse.  An important property that holds for all Lagrange interpolation filters is that their frequency response never exceeds 0 dB.  This is important in our case since stability must be ensured in the feedback network.

In general, recursive IIR filters are able to achieve the same frequency domain specifications as FIR filters with a smaller number of multiplications.  However, their design is more complicated and they can be unstable if some of the poles move outside the unit circle during real-time coefficient updates.  The z transfer function of an N^th-order all-pass filter is of the form:
 

( 4.27 )

where the numerator is the mirrored version of the denominator D(z).

The simplest IIR filter design is based on Thiran’s analytic solution for an all-pole low-pass filter with a maximally flat group delay response at the zero frequency.  Since the group delay of an all-pass filter is twice the corresponding all-pole filter, each delay value of the all-pole Thiran formulas must be doubled.  The solution for the all-pass filter coefficients approximating the delay D = N + d is:
 

for k = 0, 1, 2, … N

( 4.28 )

where  is a binomial coefficient [25].  Coefficient a_o is always equal to 1, so the polynomial is automatically scaled as desired.  Thiran also proved that for large enough delay D, the resulting all-pass filter would always be stable.  Phase delay curves of some maximally flat group delay designs are shown in Fig. 4.7 and Fig. 4.8.  We can see that the delay response is linear over a large frequency range, even for small order filters.

Fig. 4.7.  Phase delay curves of 1, 2, 3, 5, 10, and 20^th-order all-pass filters with d = 0.3 (from [25]).

Fig. 4.8.  Phase delay curves of 2^nd-order all-pass filter with d = -0.5 to 0.5 (from [25]).

for  k = 0, 1, 2, … N

( 4.29 )

where .  This approximation increases the phase distortion of the filter, but as he noted, “in some musical contexts the induced phase distortion is not subjectively objectionable.”  In these applications, using this approximation can save computations.

4.3         Previous Usages of Modulation in Artificial Reverberation

The idea of using delay length modulation in reverberation algorithms is not new.  It was first used in early hardware reverberation units such as the original EMT model 250 (released in 1975).  In those days, memory was limited and the use of modulation to randomize delay lengths was imperative to emulate the eigenmode density of a good sounding room.  Even today, many reverberation applications still have limited memory and computation resources (such as game platforms or multi-media applications).  Note that high-end commercial reverberation units such as the Lexicon 480L and the TC Electronics M3000 also use modulation to produce richer reverberation.  However, very little has been available in the literature on time-varying techniques.

4.3.1        Dattoro’s Plate Reverberator

As mentioned in section 2.7.2 when Dattoro’s plate reverberator was described, modulation can be use to reduce the coloration of the reverberation tail.  Dattorro suggested the use of low-frequency oscillators (LFO) to modulate some (or all) of the delay lines in the tank.  He suggested the use of linear (1^st order Lagrange) interpolation or even all-pass interpolation for better results.  He noted that the use of modulation was avoiding the build-up of resonances, especially when using high-frequency content input signals such as drums or percussions.  However, he also noted that modulation must be used with care to avoid audible vibrato (pitch modulation inherent to the modulation of the delay lengths) on instruments like piano.

For applications where processing power is available, he suggested the use of different modulation rates and depths for each delay line.  For applications where processing time is limited, he suggested the use of sinusoids oscillators in quadrature to modulate only two of the four delay lengths of the tank.

4.3.2        Smith’s Application Notes Related to Modulation

At the end of his first paper on artificial reverberation using digital waveguides [44], Julius Smith included in [45] a few notes on the use of modulation.  He mostly emphasized modulating the scattering coefficients of waveguide reverberation algorithms (which is equivalent to modulating the feedback matrix of a FDN network).  However, he also disscussed varying the algorithm’s delay lengths.  He compared this delay length variation to the movement of the walls of a room.  To avoid energy modulation, he suggested reducing the length of one delay line while simultaneously increasing the length of another delay line by the same amount.  We will use this method in our final implementation since it also helps mask the perceived pitch change during the reverberation tail.

4.3.3        Griesinger’s Time-variant Synthetic Reverberation

Griesinger presented a way to improve the acoustics of a room through time-variant synthetic reverberation [13].  He used artificial reverberation in an electro-acoustical sound reinforcement system to increase the perceived reverberation of a hall during musical performances.  To reduce feedback normally present in such a reinforcement system, he modulated the delay lengths in the reverberation algorithm.

He noted that interpolation was a necessity to avoid noise and clicks in the reverberation tail, as mentioned earlier.  He also noted that since delay modulation results in a pitch shift, the modulation sources must be selected with care.  For example, if RNG were used as modulators on all delay lines, all of them could go in the same direction at the same time and thus would create a noticeable shift in the pitch of the reverberation.