WaveWarp Case Study: last updated 25 April 2000

Random Reverb Case Study

Sounds Logical, April 2000

[Note: this Case Study makes use of the MATLAB® technical computing environment. MATLAB is required in order to recompute the design parameters, but is not required in order to run the associated DrawingBoards in WaveWarp. These DrawingBoards will not run in the demo version of WaveWarp - the full version is required. All files used in the Case Study (WaveWarp DrawingBoards, ASCII data files, and parameter files, MATLAB m-files and MAT files, and Microsoft Word documents plus MATLAB Word notebook (this document), are contained in the downloadable self-extracting archive "RandomReverb.exe". After downloading, double-click on the archive executable, and the setup procedure will commence. Follow the on-screen instructions, and the necessary project files will be copied into the appropriate sub-directory of the WaveWarp root directory. Note that you must have a purchased version of WaveWarp already installed, otherwise the archive cannot be restored. The included "ReadMe.txt" (automatically displayed at the end of the setup) summarises all the project files associated with the Case Study. Before using the MATLAB Word notebook or m-files, make sure that the wwmatlab sub-directory of the WaveWarp root directory is included on your MATLAB path in order that the m-files bundled with WaveWarp are accessible. To activate the MATLAB notebook, set the MATLAB working directory to the directory containing the Case Study project files then type notebook('RandomReverb.doc') at the MATLAB prompt. All MATLAB calculations included in the 'RandomReverb.doc' notebook are also included in the 'randomreverb.m' script file which can be used as an alternate reference with the advantage of being easier to modify ].

Introduction

This Case Study illustrates the use of WaveWarp in combination with MATLAB to perform designs and real-time protoyping of digital reverberation algorithms based on random FIR filters, as presented in [1]. The main idea is to construct an artificial reverberator (room simulator) based on an exponentially-decaying pseudo-random FIR filter to represent the "early reflections segment", with a feedback delay path ("around" the FIR filter) to create the dense reverberant field. The main emphasis here is on the creation of a realistic reverberant field, so the details of the early reflections - in particular the stereo field which is important in spatial perception - is not covered. Consequently, the designs presented are single channel (mono), the conventional approach when modelling the reverberant "diffuse" sound field, which, by definition, is void of any directional information.

All the designs in the Case Study closely follow those presented in [1].

Basic "room modelling" parameters

Consider a generic "large concert hall" with length, width, and height represented by RL, RW, and RH respectively. Consider nominal values given by (in MATLAB):

RL=50.0; %(m) RW=25.0; %(m) RH=18.0; %(m)

(Note: the script m-file 'randomreverb.m' included with this Case Study contains two examples: the "large concert hall" discussed here, and a small "concrete room".)

The volume V and wall area S, are given by (in MATLAB):

V=RL*RW*RH %Volume (m^3) S=2*(RL*RW+RH*RL+RH*RW) % Area (m^2)

V = 22500 S = 5200

The effective absorbing wall area, A, can be can be defined in terms of the individual absorption coefficients for each partial surface, as follows (see [2]) (in MATLAB):

%Absorption coefficients aWall=0.1; %Light drapery aCeil=0.7; %Suspended accoustic tiles aFloor=0.2; %Carpet on pad aPeople=0.5; % A=RL*RW*(aCeil+aFloor)+2*(RH*RL+RH*RW)*aWall+1000*aPeople% Effective Absorption area (m^2)

A = 1895

A useful interpretation of the effective effective absorbing wall area is to consider a hypothetical room with walls which are perfectly reflecting (i.e. an absorption coefficient of zero) except for an open window of area A through which the sounds escapes (i.e. an absorption coefficient of 1).

From these basic parameters, the Reverberation Time, T60, defined as the time taken for the sound pressure level to decrease by 60 dB due to absorption in the walls, is given (approximately) by the Sabine equation (see [2] and references therein) (in MATLAB):

T60=0.161*V/A %(sec)

T60 = 1.9116

This well-known relation, verified in practice, is founded on the simple physical reasoning that the sound energy in the room, which is proportional to the room volume and power of the sound source, is depleted at a rate proportional to the effective absorbing area.

Case I: Reverberation using a pseudo-random FIR filter alone

To start with, consider a pure FIR filter for simulating the reverberant field. In order to achieve a realistic reverb, the length of the filter should be on the order of the reverberation time, i.e. (in MATLAB):

Te=T60; %Set the FIR length equal to the reverb time (in seconds) Fs=44100; %Sample rate (Hz) Le=ceil(Fs*Te) % Length of FIR filter (in samples)

Le = 84302

The next step is to design the pseudo-random FIR filter of length Le. Following [1], an effective design is to use a pseudo-random sequence of 1's, 0's, and -1's. A key parameter in determining the quality of the reverberation is the echo density, denoted De, defined as the number of reflections per second. In terms of the FIR filter, this quantity corresponds to the density of non-zero FIR filter tap values. From psychoacoustical evaluations (see [1]), it is found that De should be no less than 2000 (reflections per second). Values of 4000 or so are preferable, and even higher values, of up to 10000, have been investigated. Furthermore, precise modelling of the reverberant field would suggest an echo density which increases in time. However, the perceptual differences observed for a quadratically increasing echo density as opposed to a constant echo density were found to be subtle and inconclusive (see [1]). Hence, for present purposes, a constant echo density will be used, as follows (in MATLAB):

De=4400; % Density (per second) of non-zero filter taps

The pseudo-random FIR filter with a non-zero filter tap density specified by De, is obtained using following MATLAB code (adapted from [1]):

p=De/Fs; % Probability of non-zero filter taps h=floor((1+p)*rand(1,Le)-p/2); % Good approximation for p<=0.2 The resulting sequence, h, is identically the impulse response of the filter. The first 10000 values are plotted in Figure 1 (via the following MATLAB code):

plot(h(1:10000));

xlabel('filter coefficient index');

ylabel('filter coefficient');

Figure 1 Initial portion of the un-weighted pseudo-random FIR filter

(Only the first 10000 filter taps have been included in the plot, since the complete filter is too large to be plotted with clarity.)

The next step in the design is to account for the wall absorption by scaling the uniform (+1/-1) filter taps by an appropriate decay function. The simplest model is an exponential decay, motivated by the fact that the sound energy is depleted by an (approximately) constant factor per reflection (see [1], [2], and references therein). This weighting envelope is computed from consideration of the reverberation time, as follows (in MATLAB):

alpha=log(1000)/T60; %Time-constant (since -60dB corresponds to 1/1000) tvec=(0:length(h)-1)/Fs; %Time vector (sec) env=exp(-tvec*alpha);

Using this envelope, the weighted FIR filter is computed as follows (in MATLAB):

hdecay=h.*env;

The first 10000 values of the impulse response of the weighted filter is plotted in Figure 2 (via the following MATLAB code):

plot(hdecay(1:10000));

xlabel('filter coefficient index');

ylabel('filter coefficient');

Figure 2 Initial portion of the exponentially-weighted pseudo-random FIR filter

Implementing the FIR filter in WaveWarp

The filter can be exported to an ASCII file ('pureFIR.txt') readable by any of WaveWarp's FIR filter modules using the wwfir command (supplied with WaveWarp), as follows:

wwfirw('pureFIR.txt',hdecay);

The WaveWarp DrawingBoard 'CaseStudyRandomReverb1.dwb' implements this pure FIR reverb design using the 'Large FIR (FFT)' module which utilises the method of fast convolution for highly-efficient real-time computation of long FIR filters. The screenshot in Figure 3 illustrates the key features of the DrawingBoard. By playing the DrawingBoard and listening to the output, it is readily apparent that the pseudo-random FIR filter produces a realistic reverberant effect, as discussed in [1]. In the next section, a truncated FIR filter in combination with feedback is used to achieve a similar result, but with significantly shorter FIR filter length.

Figure 3 Screenshot of the WaveWarp DrawingBoard 'CaseStudyRandomReverb1.dwb' showing the 'Large FIR (FFT)' module and its associated File I/O dialogue box for importing the 'pureFIR.txt' filter designed in MATLAB. Despite the long filter length (84302), it runs in real-time (on a Pentium©-II- class PC) at 44.1 kHz, processing the data from the input WAV file 'wavewarp.wav' and routing it to the soundcard ouput (in this case, a Soundblaster-compatible Yamaha DS-XG soundcard). A live input could also be processed (instead of the wav file), but the output will be delayed by the length of the FIR filter since the 'Large FIR (FFT)' module performs straightforward fast convolution without any latency-reduction segmentation techniques.

Case II: Reverberation using a short(er) pseudo-random FIR filter in combination with feedback

Following [1], a similarly natural reverberant effect can be achieved using a much shorter FIR filter (the "early segment") in combination with a feedback path to provide the slow decay. This can lead to a significantly more efficient implementation with respect to processor and memory requirements. The duration, Te, of the "early segment" should correspond (roughly) to the start time of the dense reverberations. From the room dimensions and "mean-free-path" considerations, the "start time" can be computed as follows (see [1]) (in MATLAB):

RM=4*V/S; % mean-free path (average distance between wall reflections) c=343; % speed of sound (m/s) at 20 celsius Te=1000*(4*RM/c) %starting-time for stochastic reverberation (ms) %(typically in the range from 10-200 ms)

Te = 201.8390

Expressed in samples, this corresponds to an FIR filter length, Le, given by (in MATLAB):

Le=ceil(Fs*Te/1000) % Length of early segment

Le = 8902

which is considerably shorter - by about an order of magnitude - than the full length used earlier. The filter can be taken directly from the initial segment of the previous design, as follows (in MATLAB):

hdecay=hdecay(1:Le); The filter can be exported to an ASCII file ('earlyFIRsegment.txt') readable by any of WaveWarp's FIR filter modules using the wwfir command (supplied with WaveWarp), as follows:

wwfirw('earlyFIRsegment.txt',hdecay);

Feedback

Following [1], the dense reverberations can be achieved using a simple delayed feedback loop around the FIR "early segment", as illustrated in Figure 4. The reverberation tail is thus generated by periodic repitition of the "early segment" in successively scaled-down versions.

Figure 4 Schematic for creating dense reverberations using a short FIR segment in combination with a delayed feedback path (following [1])

The feedback gain, k, and feedback delay length Lf, are design parameters. Both affect the temporal diffusion. Following [1], the feedback delay length should equal the FIR filter length (i.e. Lf = Le) such that the repitition frequency of the repeated scaled-down "early segments" is given by 1/Te. Furthermore, based on the earlier findings by the authors of [1], this sets a requirement on the "early segment" FIR filter length to be greater than (approximately) 10% of the reverberation time such that the repitition frequency is less than (approximately) 10 Hz for a reverberation time of 1 second (and scaled accordingly for different reverberation times). This contraint is roughly in accordance with the choice of FIR filter length based on the start time of the dense reverberations, discussed earlier. Likewise, the feedback gain (damping factor) should be chosen to yield the appropriate decay rate corresponding to the reverberation time. Assuming a simple "comb filter" model for the delayed feedback system, then the required gain is computed as follows (see [4 ]) (in MATLAB):

g=10^(-(3/T60)*(Le/Fs)) % For comb filter. See [4] p 382.

g = 0.4822

Implementing the FIR+feedback scheme in WaveWarp

It is straightforward to implement the feedback scheme in WaveWarp, as demonstrated in the DrawingBoard 'CaseStudyRandomReverb2.dwb'.The screenshot in Figure 5 illustrates the key features of the DrawingBoard. By playing the DrawingBoard and listening to the output, it is readily apparent that the short "early segment" pseudo-random FIR filter in combination with a feedback path, produces a realistic reverberant effect, as discussed in [1]. Figure 5 Screenshot of the WaveWarp DrawingBoard 'CaseStudyRandomReverb2.dwb' which implements the FIR+feedback scheme (from Figure 4) for creating the reverberation. Note that there is no delay explicitly included in the feedback path (as in Figure 4) since the 'Large FIR (FFT)' module performs straightforward fast convolution and thereby imposes a delay equal to the filter length. The net result is that, in comparison with the scheme in Figure 4, the overall system output is delayed by the length of the FIR filter. The feedback gain is adjusted via the 'Gain' slider (of the 'Simple gain' module) to the required value computed earlier. This example clearly illustrates how a signal flow diagram as depicted in Figure 4 can be easily and intuitively implemented by virtue of WaveWarp's sample-by-sample modular architecture.

Low-pass filtering

Following [1], it is proposed that the "(subjective) quality" of the reverb in the previous section can be significantly improved by incorporating a lowpass filter in the feedback path. Specifically, a single-pole filter with a pole location between 0.4 and 0.7 (in the z-domain) is suggested as preferable, based on a limited set of listening experiments (and depending on the reverberation times and on the actual type of signal, i.e., speech, trumpet, etc). This filter can be designed as follows (in MATLAB, using the default pole location of 0.4):

P=0.4; % desired pole location (default value from [1]) num=(1-P); % discrete numerator polynomial for unity dc gain den=[1 -P];% discrete denominator polynomial for desired pole location

The frequency response of this filter is computed and displayed as follows (in MATLAB):

freqz(num,den,100);% Compute and plot the frequency response

Figure 6 Frequency response of single-pole lowpass filter to be used in feedback path of modified reverberation scheme

The filter can be exported to an ASCII file ('lowpassIIR.txt') readable by WaveWarp's 'IIR (from file)' filter module (for IIR filtering in the 'direct' structure) using the wwiirdw command (supplied with WaveWarp), as follows:

wwiirdw('lowpassIIR.txt',num,den);

Implementating the FIR+feedback+low-pass filtering scheme in WaveWarp

In WaveWarp, it is straightforward to incorporate the lowpass filter in the feedback path, as demonstrated in the DrawingBoard 'CaseStudyRandomReverb3.dwb'.The screenshot in Figure 7 illustrates the key features of the DrawingBoard. By playing the DrawingBoard and listening to the output, it is readily apparent that the lowpass filter improves the results compared with the previous case (i.e. without the lowpass filter). This confirms the findings discussed in [1]. Moreover, as also mentioned in [1], the output of the FIR filter can also be processed with a lowpass filter (e.g. to simulate absoprtion by the air in the room). This is an almost trivial modification to implement in WaveWarp, as demonstrated in the DrawingBoard 'CaseStudyRandomReverb4.dwb' (and illustrated in the screenshot in Figure 8) which includes an additional lowpass filter downstream of the FIR filter (with the same default settings as the lowpass filter in the feedback path, though this is simple to change). By listening to the output of this DrawingBoard, it is observed that the quality of the reverb is improved by the inclusion of this additional lowpass filter.

Figure 7 Screenshot of the WaveWarp DrawingBoard 'CaseStudyRandomReverb3.dwb' which implements the FIR+feedback scheme (from Figure 4 and Figure 5) modified to include a single-pole lowpass filter in the feedback path. The lowpass filter is implemented by the 'Filter IIR (direct) from file' module which has been confgured to import the 'lowpassIIR.txt' filter designed in MATLAB. Also shown is WaveWarp's frequency response plot for the filter (which is the same as that generated by MATLAB in Figure 6).

Figure 8 Screenshot of the WaveWarp DrawingBoard 'CaseStudyRandomReverb4.dwb' which implements the FIR+feedback+lowpass scheme (from Figure 7) modified to include an additional single-pole lowpass filter after the FIR filter. Following[1], the same filter design is used in both lowpass filters (though this can easily be changed).

Measurement of the overall impulse response using WaveWarp

It is instructive to investigate the overall impulse response of the FIR+feedback+lowpass reverberation scheme. This is straightforward in WaveWarp, as demonstrated in the DrawingBoard 'CaseStudyRandomReverb5.dwb' (and illustrated in the screenshot in Figure 9).

Figure 9 Screenshot of the WaveWarp DrawingBoard 'CaseStudyRandomReverb5.dwb' which implements a measurement of the impulse response of the FIR+feedback+lowpass scheme (from Figure 8). The measured response is saved in the ASCII file 'measurement.txt' for off-line processing. The response is also is also captured in real-time via the WaveWarp oscilloscope module (shown).

The output of the impulse response measurement is stored in the ASCII file 'measurement.txt' which can be imported and plotted in MATLAB as follows:

load -ascii 'measurement.txt'; %read the ouput file created by WaveWarp plot((0:length(measurement)-1)/Fs,measurement); xlabel('Time (s)'); ylabel('Measured impulse response');

Figure 10 WaveWarp measurement of the impulse response of the pseudo-random-FIR+feedback+lowpass reverberation scheme. The exponential decay of the diffuse field is readily apparent.

From the plot of the measured impulse response, the general behavior of the algorithm can be observed: namely, the reponse starts with the FIR "early segment", followed by successive scaled-down repititions. The overall exponential decay lasts for a duration characterised by the reverberation time (approximately 2 seconds using the default values in this example).

Conclusions

The random FIR model with feedback, as presented in [1], is confirmed to be a simple and effective approach to generating realistic artificial reverberation. The Case Study demonstrates the effective combination of WaveWarp - for real-time prototyping of signal processing algorithms - and MATLAB - for rapid design and analysis.

References

[1] Artificial reverberation based on a Pseudo-random Impulse Response II, Per Rubak, Lars G. Johansen, Paper No. 4900(G6), Presented at the AES 106th Convention, May 8-11 1999, Munich, Germany. [2] The Science Of Sound (second edition), Thomas D. Rossing, Addison Wesley, 1990 (reprinted) [click here to purchase]. [3] Digital Audio Signal Processing, Udo Zolzer, Wiley, 1997[click here to purchase].

[4 ] Elements of Computer Music, F. Richard Moore, PTR Prentice Hall, 1990.