2024年7月31日发(作者：畅靖之)
Image method for efficiently simulating small-room acoustics
Jont B. Allen and David A. Berkley
Acoustics Research Department, Bell Laboratories, Murray Hill, New Jersey 07974
(Received 6 June 1978)
Image methods are commonly used for the analysis of the acoustic properties of enclosures. In this paper
we discuss the theoretical and practical use of image techniques for simulating, on a digital computer, the
impulse response between two points in a small rectangular room. The resulting impulse response, when
convolved with any desired input signal, such as speech, simulates room reverberation of the input signal.
This technique is useful in signal processing or psychoa½oustic studies. The entire process is carried out on
a digital computer so that a wide range of room parameters can be studied with accurate control over the
experimental conditions. A •oaza^• implementation of this model has been included.
PACS numbers: ,
INTRODUCTION
(1) We are most interested in the office environment,
which is usually a rectangular geometry.
In some recent experiments, which studied the per-
ceptual effects of reverberation properties of a small
room, t'2 a carefully controlled, easily changed, acou-.
stic environment was required. It was decided to utilize
a computer simulation of the acoustic space. This pa-
per describes both the general theoretical approach and
the specific implementation techniques used (the
FORTRAN program). We believe that the resulting room
model for a broad range of investigations,
from our original experiments mentioned above, to
basic studies of room acoustics.
(2) This model can be most easily realized in an ef-
ficient computer program.
(3) The image solution of a rectangular enclosure
rapidly approaches an exact solution of the wave equa-
tion as the walls of the room become rigid.
The image model is chosen because we are interested
in the point-to-point (e.g., talker-to-microphone) trans-
fer function of the room. In order to obtain a good
transient description of the response, a time domain mo-
del is required. A normal-mode solution of the enclosure
would require calculation of all modes within the fre-
The room model assumed is a rectangular enclosure
with a source-to-receiver impulse response, or trans-
fer function, calculated using a time-domain image ex-
pansion method. Frequent applications have been made
of the image method in the past as in deriving the re-
quency range of interest (i.e., 0.1-4.0 kHz), plus cor-
rections for those outside this range. The image meth-
od includes only those images contributing to the im-
pulse response. Thus the contributing images are those
within a radius given by the speed of sound times the
verberation-time equations, 3 for theoretical studies of
sound behavior in enclosures, 4'? and in the study of
architectural acoustics and perceptual properties of
reverberation time. 6 (The exact relationship between
the normal-mode solutions and the image solution, for
a lossless room, is discussed in Appendix A.) The im-
portant information used here is that in the time-do-
main, each image contributes only a pure impulse of
known strength and delay while each normal mode is
a decaying exponential which contributes to all times.
Furthermore, whereas an image has only delay and
gain as parameters, a normal mode computation re-
quires the solution of transcendental equations to find
the pole location plus the evaluation of a relatively
complex function to find the mode gain (i.e., the residue
of the pole).
rooms. 8'11 In addition, there has been a considerable
amount of important theoretical work on the approxi-
mate 12 use of images produced by a single soft-wail
(finite impedance) reflection. Several recent papers on
this subject which have good bibliographies are Refs.
13, 14, and 15. Computer methods have also recently
been applied to image computations in enclosures (see
for example Refs. 6, ?, 10, and 11). In the current
paper the computational technique is specifically aimed
at being simple, easy to use, and fast. In addition the
resulting'room responses have been used to realisticai-
1¾ mode[ speech transmission in rooms and to investi-
gate the effects of various forms of digital speech sig-
nai processors .t6, l ?
In the following we will first briefly discuss theoretic-
al aspects of the method. Then we will outline the compu-
tational approach and, finally, we will give some ex-
amples of applications.
I. IMAGE VERSUS NORMAL MODE MODELS
A. The image model
We model a talker in a room as a point source in a
rectangular cavity. A single frequency point source of
acceleration in free space emits a pressure wave of the
form
We model the rooms of interest as simple rectangular
enclosures. This choice of geometry is made for sev-
eral reasons-
P(co,X,X') exp[{o•(R/c - t)]
= 4•.R '
where
(t)
943 J. Acoust. Soc. Am. 65(4), Apr. 1979 0001-4966/79/040943-08500.80 (D1979 Acoustical Society of America 943
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
P = pressure,
IMAGE EXPANSION p(t)
•y
c• =2•rf,
f = frequency,
t = time,
X X
X X
X X
X X
0
XX
XX
X X
XX
XX XX
XX
X X
X X
XX
XX
R= Ix-x'l,
x =vector talker location (x,3,z),
(2)
XX
X' = vector microphone location (x', 3 ', z'),
i=v•/,
c =speed of so•d.
8 •o t-
p(t)=• •.
c
I 8( IRp+
p=lr=-m
When a rigid wall is present, the rigid wall (zero norm-
al velocity) boundary condition may be satisfied by
p•ci• an im•e symmetrically on the far side of the
wall. Thus,
R r=(2nLx,21Ly,2mLz)
Rp= (x- + x', y+_ y',z_ + z')
Rp-I- R r = (x_ + x'+ 2nLx,y + y'+ 21Ly,z*_zø+2rnLz )
FIG. 1. A slice through the image space showing how the
The solid box
represents the original room. The actual image space is
three dimensional.
R(•,X,X') [exp[i(•/c)R+] )
images of the source are spatially arranged.
=L 4•R + +
exp[i(•/c)R..]]exp(_{•t
4•R.. J
(3)
where we define the two distances from the microphone
to the source R. and to the im•e R. by
a z =(•- •')z + (3 -3')z + (z - z') z
' '
a z =(x+ •)z +(y -y')z +(z -z') z
+
(4)
be derived directly from the normal-mode solution as
shown in Appendix A.
The wall has been placed at x=0 in this case (note the
s•n in the x terms of R+ and R.).
In the general case of sk walls the sit•tion becomes
more complicated because each image is itself imaged.
The pressure may then be written (as shown in Appen-
dk A)
B. Case of nonrigid walls
If the room walls are not rigid, the soluHon in terms
of point images may no longer be exact. A precise
statement of the effects of finite impedance walls is
presently impossible, since the effects on even a single
image are quite complicated. la'14'15 Therefore we have
continued to assume the approximate point image model
even for nonrigid walls. In addition we have assumed
an angle independent pressure wall reflection coefficient
fl. This assumption is equivalent to assuming that the
R(•,X,X')= •
permutations over ß of
4•1• +R .I '
(6)
where • represents the eight vectors given by the eight
wall impedance is proportional to sec(e), where e is the
and
angle of incidence of a plane wave with respect to the
wall normal. We presently do not understand the exact
physical interpretation of the above assumptions. How-
ever, we believe that they do not introduce serious
problems into the final result under typical conditions.
By typical, we mean over the frequency range of 100
Hz-4. kHz, wall reflection coefficients of greater than
0.?, typical office room geometries, and where both
source and receiver are not close to the wall. Many,
•=(x•x',
•=2(.n•.
y •y', z•z')
ZL•. •n•).
r is the integer vector triplet (•,•,m),
(V)
where (Lx, L•, L•) are the room dimensions. Equation
(5) is the pressure frequency response ass•ing rigid
walls for a point source at X =(x,y,z) and receiver at
X'=(• ,3',z'). If Eq. (5) is Fourier transformed, we
find the room impulse response function (time domain
Green' s f•ction)
p(t,X ,X' )= :.- •1• + •1 '
•r,• •[t-(I R•+ • l/c)] (8)
if not all, of the above conditions are probably not cri-
tical and could be relaxed. We merely wish to carefully
point out the nonexact nature of the results.
The above assumptions result in the Sabine energy
absorption coefficient a for a uniform reflection coeffi-
cient fi on a given wall of the form
An interpretation of Eq. (8) is given in Fig. 1 where we
show a part of the image space for a two-dimensional
slice through the room. When the accelerative source
a = t - f. (O)
location (talker) X is excited, each image point is si-
multaneously excited, creating spherical pressure
waves which propagate away from each image point.
Our assumptions are similar to those of geometrical
acoustics 18 and are the same as those required for spec-
ular angle-independent ray tracing. In current imple-
mentations of the model we also do not allow frequency
variations in the reflection coefficients. Both the angle
dependence and frequency dependence could be included
944
Equation (8) is the exact solution to the wave equation
in a rectangular, rigid-wall (lossless), room and may
944 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979
J. Allen and D. Berkley: Method for simulating smalr-room acoustics
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
in our computer program, but only at the expense of
significantly complicating and slowing down the compu-
tational model.
TABLE I. Computation parameters--room size (feet) 10' x 15'
x 12.5' 8 kHz sampling rate.
Introducing the effects of finite, angle independent
Length
Impulse response
Image Computation
Convolution
rate
wall absorption into Eq. (8) leads to the modified room
impulse response
1 .•
(ms)
64
128
256
No. points
512
1024
2048
count
585
4690
37500
time (s)
i
8
60
(s/s)
12.5
13.8
15.0
•(t,X,X')=• •-• •ln'ql•lnl•ll'jl•lll•lm'kl• Irnl
•xl Px2 P•I •2 P•I •2
(10)
where • is now expressed in terms of the integer 3-
vector p = (q,j, k) as
I for our implementation on a Data General, Eclipse
•=(x-x • + 2qx', y -y' + 2jy', z-z' + 2kz') . (11)
• as given by Eq. (6) is similar to that of (11), but is
indexed differently from (11). The hera's are the pres-
sure reflection coefficients of the s• boundary planes,
with the subscript 1 referring to walls adjacent to the
S/200 computer. (On this machine the computation
time required for each image is about 1.6 ms.) The
actual FORTRAN programs used are given in Appendix B.
coordinate origin (see Fig. 1). Subscript 2 is the oppos-
ing wall. Eq. (10) has beenderiv• heuristically from
The temporal quantization in the impulse function
computation causes slight statistical errors in the com-
puted arrival times of each image pulse relative to the
geometric• consideratio• of Fig. 1. The sum • with
vector ind• p is used to indicate three sums, namely
one for each of the three components of p=(q,j,k).
r=(n,l,m) is a similar sum. Physically these s•s
are over a three-dimensional lattice of •ints. For p
there are eight points in •e lattice and for r, the lattice
is i•inite.
exact delay as given in Eq. (10). This error can be
thought of as effectively "moving" each image source by
0• ezzoz• AR/2 relative to the receiver. This effect
could be removed, in principle, by using a band-limited
source pulse. However, the error is small for most
(if not all) purposes and it greatly complicates-the com-
putation to remove this approximation. We have esti-
mated that the error due to the slight moving of the im-
ages could not be perceived even in a digital simulation
of a binaural hearing experiment.
II. IMPLEMENTATION OF THE MODEL
The primary consideration in a computer (sampled
data) implementation of Eq. (10) is the method of spatial
sampling. In addition, an apparently nonphysical be-
havior of the model at zero frequency is removed by a
The subroutine SROOM of Appendix B requires as pa-
rameters the number of impulse response points de-
sired (NPTS), the source location R0, the receiver lo-
low-frequency (0.01 of the sampling frequency) high-
pass digital filter. •9
t.0
A calculated impulse response is built up as a "histo-
gram" of image pulses received at different time de-
lays. The width of each histogram bin is equal to the
time sampling period T initially assumed, which in
turn is determined by the highest frequency to be rep-
resented. For example, all images with the range
N•R to (N + 1)AR, where •R=cT (T is the sampling
period and c the speed of sound), are added together
IMPULSE RESPONSE
with appropriate amplitude as given by Eq. (10).
The choice of sampling rate is the appli-
cation. If speech is to be studied in small rooms one
might choose T =0.1 ms. (sampling frequency of 10
kHz; highest frequency of 5 kHz). But, if reverberation
times of large enclosures are being studied (and convo-
lution with speech is not required) much lower rates
can be useful.
- 1.O
2048 POINTS
8 KHZ SAMPLING
I I I
RATE
I
The time length of the calculated impulse response is
also a consideration. For a given sampling rate the
0 TIME (MS) 256
number of points in p(t) increases linearly with its
length while the computation time (and number of im-
ages) goes up approximately as the cube of response
length. This is shown in the first four columns of Table
945 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979
FIG. 2. Plot of a typical impulse response for a room
80 x120 x100 sample lengths long. Wall reflection coefficients
were all 0.9 ceiling and floor coefficients were 0.7. X and X'
were at (30,100, 40) and (50, 10, 60) sample periods.
J. Allen and D. Berkley: Method for simulating small-room acoustics 945
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
cation R, the room dimensions, all specified in terms
use a modification of the integrated tone-burst method 22
of the sample length (AR), and the reflection coefficients
of each of the six wall surfaces (•). Figure 2 shows an
example of the impulse response obtained for a room of
dimensions 80 x 120 x 100 sample lengths with equal wall
reflection coefficients of 0.9 (c• =0.19) and with floor
and ceiling reflection coefficients (•) of 0.7 (a -- 0.51). X
and X' were (30,100,40)and (50,10,60)sample
tervals, respectively.
where E(t) is the average energy decay, k is a propor-
tionality constant, and p(•) is the calculated pressure
impulse response from Eq. (10). For cases where the
impulse response has been truncated before most of the
in-
It is usually convenient to interpret the model param-
eters as a true distance rather than as multiples of
AR. This requires the choice of a sampling rate and
then conversions may be performed in the users main
program which calls the subroutines of Appendix B.
Figure 2 is labeled assuming an 8 kHz sampling rate.
decay has taken place, (12) may lead to errors. These
errors are usually obvious in the E(t) plots.
Another, approximate, approach is to simply mea-
sure the short-time average energy decay of the im-
pulse itself (e.g., using a simulated level recorder).
For exponential or near-exponential decays, both
methods should give approximately the same value of
For this assumption (and assuming a sound speed of 1
ft/ms) the room dimensions are 10' x 15' x12.5'.
III. APPLICATIONS
Our room image model has been applied to several
physical evaluation of room reverberation effects 1'2 and
a study of critical distance measurements using spec-
problems. We will discuss two examples: a psycho-
(o)
tral response variance. zø We have also used the model
to test a signal processor intended to reduce perceived
reverberationl• 6 and to study problems associated with
mathematical inversion (inverse filtering) of room
transfer functions.• ?
A. Psychophysics of room reverberation
Once a simulated room impulse response has been
calculated using the image model, the psychophysical
effects of this simulated reverberation on speech may
be directly studied. A reverberant sample of speech
DECAY CURVE
2048 POINTS
8 KHZ SAMPLING RATE
was produced by convolving an anechoic (unreverberant)
speech sample with the calculated impulse response
[P(t)l. This can be done efficiently using a Fast Fourier
Transform (FFT) method (overlap-add) to perform the
convolution. zi The last column of Table I shows the
measured convolution rate, for various length impulse
responses. The convolution rate only increases as
TIME (MS)
256
logz(N), (where N is the room response length in time
sample periods T) so even large impulse responses can
be convolved with speech quite efficiently. For exam-
ple, to convolve (filter) one second of speech, sampled
at 8 kHz, with a 256 ms long impulse response (2048
points) requires a 15 s computation.
ENERGY DECAY CURVE
The speed of processing makes multivariate psycho-
physical studies quite practical. Ease of modification
and perfect control of room parameters avoids the
problems which have made such experiments so difficult in
the past. The actual experiments used 16 different simu-
lated" rooms" (impulse res ponses) convo lved with ten dif-
ferent sentences spokenby four different speakers. Itwas
discovered that the experimental rooms were perceptually
8 KHZ. SAMPLING
I ,
RATE
, I, , I ,,
256
well characterized by their spectral variance [Eq. (14)]
andby the reverberation time. This latter measure,
reverberation time, deserves some discussion.
, i , i •l, ,, I ..... 1
o
TIME (MS)
Given the impulse responses, reverberation time may
be estimated in a number of ways. One method is to
946 J. Acoust. $oc. Am., Vol. 65, No. 4, April 1979
Fig. 2 using the Schroeder integration method. 22 (b) Impulse
energy decay curve for a simulated level recorder.
946
FIG. 3. (a) Energy decay curve for the impulse response of
J. Allen and D. Berkley' Method for simulating small-room acoustics
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
reverberation time. Example plots of E(t) for both
methods are shown in Figs. 3(a) and 3(b) using the im-
pulse response of Fig. 2. Experience indicates thatEq.
(12) gives the most satisfactory results. We have found
empirically that calculated reverberation times, for a
number of simulated enclosures, agree well with
Eyring's formula s over a wide range of Beta values. 23
In the experiments discussed above 1' •-4 we discovered a
monotonic relationship between A/A c( Fig. 4)the micro-
(a) STANDARD DEVIATION OF SPECTRAL RESPONSE
IN COMPUTER-SIMULATED ROOM: 17' x !$' x
7
phone-talker distance when normalized by the room c ri-
tical distance (the distance at which reverberant energy
equals direct sound energy), the reverberation time, and
psychophysical preference for the resulting speech.
B. Critical distance measurement
z
o
n- 5
. A new method has been proposed •'ø for measurement
of critical distances (or reverberation radius) in rooms.
In this technique a measurement is made of the log frequen-
cy response variance gL defined as
u•4
L(w)=20 logtiP(w) [ ] (13)
(14)
z
o
z
given the room transfer-function P(•o) [Fourier trans-
form of Eq. (10)] for several microphone-source spac-
ings. The measured values are fitted to a theoretical
curve for (rL based on the assumption of simultaneously
excited, uncorrelated, normal modes, combined with
the calculated direct sound energy. The resulting fit
was shown to give an accurate value for the room's
critical distance.
o
z
This new method was extensively studied using our
0
-$0 -24
I
-18
I
-12
I
-6
i
0
I
6
[
12:
I
18
I
image model before being applied successfully to real
24
REVERBERANT/DIRECT ENERGY RATIO (dB)
rooms. Since the direct and reverberant energy are
(•/•c)'
(•) STANDARD DEVIATION OF SPECTRAL RESPONSE
I1• COMPUTER-SIMULATED ROOM; 47' x :51' x 15'
7•
known in the computer model, a comparison can easily
be made to the theory. The model results show excell-
ent agreement with theoretical calculation as is seen in
Figs. 4(a) and 4(b). We know of no other method by
which this study could have been carried out as effect-
ively.
• 6- o
•
f•
0
&•l•• • 0
p ß o.•
.oo
ß
ß o n
IV. SUMMARY
A simulation
AND DISCUSSION
method for small rooms based on an
approximate image expansion for rectangular nonrigid-
wall enclosures has been discussed. The method is
simple, easy to implement and efficient for computer
simulation. Several examples of its use, where other
methods would be difficult, have been discussed.
APPENDIX A
o _
-
>
We wish to derive the rigid-wall image solution di-
rectly from the normal-mode expansion for a rectangu-
lar enclosure. The frequency response function
(Green's function) for the pressure P(•o) in an enclos-
ure is given by solving the Helmholtz equation driven
by a single frequency point acceleration source.
I
t2
0 I I
-18
I
-12
, I
-6
I
0
I
6
I
18 24 -E•) -24
REVERBERANT/DIRECT ENERGY RATIO (dB)
,x'] + ,x'] =- ,(x -x '),
(A1)
(A/Ac)2
FIG. 4. Figures from Jetzt 2ø which compare the theoretical
rms deviation of the pressure in dB from the mean pressure
in dB as a function of the direct to reverberent energy ratio
(a) for a room 17x13x10 ft and (b) 47x31 x15 ft.
where • is the frequency and c is the speed of sound,
The solution to this equation, assuming rigid boundar-
ies, is given by
947 947 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
P(k,X,X')=• =.• (kZr_k z) ,
•,• •,(x)%(x' )
where k=o•/c, r=(n,l,m)
al sum, V is the room volume,
(A2)
By Fourier series analysis one may show
indicates a three dimension-
-,, ß
• ( n•)-L--•exp(i2L n•x)
(A9)
Thus [with analogus equations to (A9) for y and z ]
, , LE '
kZr= Ik, I2
and
(A3)
f
• ,=-•
• exPt(i•'(Rv•z R•)] d•
(•z_ ) ,
(A10)
where • is the vector [also given by Eq. (7)]
/n,rXcos[l•ry /m•rz
]
•=2(•L,, •L•, mL•). (A•)
where the L•'s are the room dimensions.
Each triple integral is just a plane wave expansion for
a •int source in free space since
Using the exponential expansion for cosine, multi-
plyi• the terms of Eq. (A2) together and collecting, we
obtain
P(k,X,X') =•,•
•q• (•)]
1 • •exp(•-•)
• (•-•:)
exp(ikl RI ) ='•'x• (I tj I z - k •)
4•1RI '
(A12)
,
(A5)
(A6)
Finally, using Eq. (A12), Eq. (A10) becomes
4,tr I R• + l:{rl
p(,• ,X,X,)=••__ r__•.•exp[i(•/c)l Rp + R•I]
(A13)
where • represents the eight vectors [also given by
• =(x •x', y •y', • •').
Taking the inverse Fourier transform ofEq. (A13), the
echo structure becomes explicit
Using the property of the delta Mnction onk,, k•, and k,
f" •(• - a)r(•)d• =r(a),
we may rewrite Eq. (A5) in integr• form
(A7)
(AO)
1 •
APPENDIX
C PGM:
fexp(if-•) •
, (A14)
=_• 4•' 11• + R,.I
p(t,X,X')=••__ r•'• 6[t-(lRp+Rr[/c)]
which is the same as Eq. (8) as desired.
B
SROOM
C SUBROUTINE TO CALCULATE A ROOM IMPULSE RESPONSE
C R=VECTOR RADIUS TO RECEIVER IN SAMPLE PERIODS=LENGTH/(C*T)
C R0 =VECTOR RADIUS TO SOURCE IN SAMPLE PERIODS
C RL=VECTOR OF BOX DIMENSIONS IN SAMPLE PERIODS
C BETA=VECTOR OF SIX WALL REFLECTION COEFS (0 C HT =IMPULSE RESP ARRAY C N-PTS=# OF POINTS OF HT TO BE COMPUTED C ZERO DELAY IS IN HT(1) C SUBROUTINE SROOM (R, R0, RL, BETA, HT, N-PTS) DIMENSION HT (N-PTS) DIMENSION R(3), R0(3), NR(3), RL(3), DELP(8), BETA(2,3) EQUIVALENCE (NR(1), NX), (NR(2), NY), (NR(3), NZ) DO 5 I= 1, N-PTS 5 HT(I) =0 C CK FOR MIC AND SOURCE AT SAME LOCATION DIS =0 6 DO 6 I=1,3 DIS = (R(I)-R0 (I))*'2 +DIS DIS = SQRT(DIS) IF (..5)HT(1) = 1 IF (..5)RETURN RANGE OF SUM C FIND N1 = N-PTS/(RL(1)* 2) + 1 N2 = N-PTS/(RL(2)* 2) + 1 N3 = N-PTS/(RL(3)*2) + 1 DO 20 NX=-N1, N1 DO 20 NY =--N2, N2 DO 20 NZ =--N3, N3 948 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics 948 Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17 C GET EIGHT IMAGE LOCATIONS IO=O FOR MODE # NR CALL LTHIMAGE (R, R0, RL, NR, DELP) DO 10 L = 0,1 DO 10 J =0, 1 DO 10 K = O, 1 IO=IO + 1 C MAKE DELAY =ID AN INTEGER ID = DELP (I0) + .5 FDM1 ID =ID +1 IF ()GO C PUT IN LOSS FACTOR TO 10 ONCE FOR EACH WALL REFLECTION Gm= BE TA(1,1)** IABS(NX-L) 1 *BETA(2,1)**IABS(NX) 2 *BETA (1,2)** IABS(NY-J) 3 *BETA(2,2)**IABS(NY) 4 *BE TA (1,3**IABS (NZ-K) 5 *BETA(2,3)**IABS(NZ) 6 /FDM1 C CHECK FOR FLOATING POINT UNDERFLOW HERE; C IF UNDER FLOW, SKIP NEXT LINE HT(ID) = HT(ID) +GID 10 20 CONTINUE CONTINUE NOTE LINES CONTINUATION C IMPULSE RESP HAS BEEN COMPUTED C FILTER WITH HI PASS FILT OF 1% OF SAMPLING FREQ (I.E'. 100 HZ) C IF THIS STEP IS NOT DESIRED, W = 2.*4.*ATAN(1.)* 100. T = 1E-4 RETURN HERE R1 = EXP(-W*T) R2 =R1 B1 = 2.* RI* COS(W* T) B2 =-RI*R1 A1 =-(1. +R2) A2 =R2 Y1 =0 Y2 =0 Y0 =0 C FILTER HT DO 40 I=1, NPTS X0 =HT(I) HT(I) =Y0 +AI*YI+A2*Y2 Y2 =Y1 Y1 =Y0 Y0 =BI*Y1 40 CONTINUE RETURN END +B2*Y2 +X0 C PGM: C PGM C LTHIMAGE TO COMPUTE EIGHT IMAGES OF A POINT IN BOX SUBROUTINE LTHIMAGE (DR, DR0, RL, NR, DELP) C C DR IS VECTOR C DR0 IS VECTOR C RL IS VECTOR C NR IS VECTOR C DELP C C RADIUS RADIUS TO RECEIVER TO SOURCE IMAGE IN SAMPLE IN SAMPLE IN SAMPLE TO IMAGE PERIODS PERIODS PERIODS OF BOX DIMENSIONS OF MEAN IN SAMPLE NUMBER IS VECTOR OF EIGHT SOURCE PERIODS DISTANCES DIMENSION R2L(3), RL(3), NR(3),DELP(8) DIMENSION DR0(3), DR(3), RP(3,8) C LOOP OVER ALL SIGN PERMUTATIONS AND COMPUTE R +/- R0 I0= 1 DO 10 L =--l,1,2 DO 10 J =--1,1,2 DO 10 K=-l,1,2 949 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics 949 Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17 C NEAREST IMAGE IS L = J = K =-1 RP (1, I0)=DR(1)+ L'DR0(1) RP(2,I0) =DR(2) +J. DR0(2) RP(3, I0) =DR(3) + K. DR0(3) I0 = I0 + 1 10 CONTINUE C ADD IN MEAN RADIUS TO EIGHT VECTORS TO GET TOTAL DELAY R2L (1) = 2.* RL (1)* NR(1) R2L (2) = 2.*RL (2)* NR(2) R2L ($) = 2.* RL (3)* NR(3) DO 20 I=1,8 DELSQ=0 DO 25 J = 1,3 R1 = R2L (J)-RP(J,I) DELSQ =DELSQ +Rl**2 25 CONTINUE DELP (I) = SQRT(DELSQ) 20 CONTINUE RETURN END 1Barbara McDermott and Jont Allen, "Perceptual Factors of Small Room Reverberation," J. Acoust. Soc. Am. 60, S9 (A) (1976). an Impedance Boundary," J. Acoust. Soc. Am. 59, 780-785 (1976). 16R. J. Donato, '•Spherical Wave Reflections from a Boundary 2jont B. Allen and B. J. McDermott, "A New Method for Mea- suring Perception of Room Reverberation" (unpublished). 3C. F. Eyring, "Reverberation Time in 'Dead' Rooms," J. Acoust. Soe. Am. 1, 217-241 (1930). 4D. Mintzer, "Transient Sounds in Rooms," J. Acoust. Soc. 51{. H. Bolt, P. E. Doak, and P. J. Westervelt, "Pulse Sta- tistics Analysis of Room Acoustics," J. Acoust. Soc. Am. Am. 22, 341-352 (1950). 22, 328-339 (1950). 16Jont B. Allen, D. A. Berkley, and J. Bauert, '•lViultimicro- phone signal-processor technique to remove room reverbera- tion from Speech Signals," J. Acoust. SOc. Am. 62, 912-915 (1977). of Reactive Impedance using a Modification of Cagniard's Methods," J. Acoust. Soc. Am. 60, 999-1002 (1976). ?S. T. Neeley and Jont B. Allen, "Invertibility of a Room Im- pulse Response," J. Acoust. Soc. Am. (in press). 18p. M. Morse and K. U. Ingard, Theoretical Acoustics (McGraw-Hill, New York, 1968), Sec. 9.5. 6j. M. Berman, "Behavior of Sound in a Bounded Space," J. Acoust. Soc. Am. 5?, 1275-1291 (1975). 1E. K. Dunens and R. F. Lambert, "Impulsive Sound Response Statistics," J. Acoust. Soc. Am. 61, 1524-1532 (1977). 6J. R. Power, "Measurement of Absorption in Rooms with Sound Absorbing Ceilings," J. Acoust. Soc. Am. 10, 98-101 (1938). SC. G. Mayo, '•tanding Wave Patterns in Studio Acoustics," Acustica 2, 49-64 (1951). lSWithout the high-pass filter, the Fourier transform ofp(t) shows a large energy spike at zero frequency as the absorp- tion a goes to zero. This is in part a direct result of the no- pressure-release boundary conditions we have assumed and is •urther complicated by the assumption of an accelerative source, which is nonphysical for very low frequencies (less than 50 Hz). 1øB. M. Gibbs and D. K. Jones, "A Simple Image Method for Calculating the Distribution of Sound Pressure Levels within an Enclosure," Acustica 26, 24-32 (1972). 20j. j. Jetzt, "Critical Distance Measurement of Rooms from the Sound Energy Spectral Response," J. Acoust. Soc. Am. (to be published). llF. Santon, "Numerical Prediction of Echograms and the In- telligibflity of Speech in Rooms," J. Acoust. Soc. Am. 59, 1399-1405 (1976). 12By "approximate" we mean the blurring of the image as a re- sult of the finite impedance wall. Variations in wall imped- ance as a function of frequency or other frequency dispersion effects will also make the images frequency dependent. In practice the two effects usually appear simultaneously. 21Jont B. Alien, "FASTFILT-An FFT Based Filtering Pro- gram, "in IEEE Press book on Programs for Digital Signal Pro- cessing (IEEE Press, New York, 1979). 22M. R. Schroeder, "New Method of Measuring Reverberation Time," J. Acoust. SOc. Am. 37, 409-412 (1965). 2SAgreement is best for beta the same on all surfaces and variations appear when a pair of opposing walls are signifi- canfly different in reflectivity than all other surfaces, (as is 13A. R. Wentzel, "Propagation of Waves Along an Impedance Boundary," J. Acoust. Soc. Am. 55, 956-963 (1974). 14S. Thomasson, '•eflection of Waves from a Point Source by 24D. A. Berkley, "Normal Listeners in Typical Rooms:Rever- beration Perception, Simulation and Reduction" (unpublished). the case for Figs. 2 and 3). 950 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics 950 Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
2024年7月31日发(作者：畅靖之)
Image method for efficiently simulating small-room acoustics
Jont B. Allen and David A. Berkley
Acoustics Research Department, Bell Laboratories, Murray Hill, New Jersey 07974
(Received 6 June 1978)
Image methods are commonly used for the analysis of the acoustic properties of enclosures. In this paper
we discuss the theoretical and practical use of image techniques for simulating, on a digital computer, the
impulse response between two points in a small rectangular room. The resulting impulse response, when
convolved with any desired input signal, such as speech, simulates room reverberation of the input signal.
This technique is useful in signal processing or psychoa½oustic studies. The entire process is carried out on
a digital computer so that a wide range of room parameters can be studied with accurate control over the
experimental conditions. A •oaza^• implementation of this model has been included.
PACS numbers: ,
INTRODUCTION
(1) We are most interested in the office environment,
which is usually a rectangular geometry.
In some recent experiments, which studied the per-
ceptual effects of reverberation properties of a small
room, t'2 a carefully controlled, easily changed, acou-.
stic environment was required. It was decided to utilize
a computer simulation of the acoustic space. This pa-
per describes both the general theoretical approach and
the specific implementation techniques used (the
FORTRAN program). We believe that the resulting room
model for a broad range of investigations,
from our original experiments mentioned above, to
basic studies of room acoustics.
(2) This model can be most easily realized in an ef-
ficient computer program.
(3) The image solution of a rectangular enclosure
rapidly approaches an exact solution of the wave equa-
tion as the walls of the room become rigid.
The image model is chosen because we are interested
in the point-to-point (e.g., talker-to-microphone) trans-
fer function of the room. In order to obtain a good
transient description of the response, a time domain mo-
del is required. A normal-mode solution of the enclosure
would require calculation of all modes within the fre-
The room model assumed is a rectangular enclosure
with a source-to-receiver impulse response, or trans-
fer function, calculated using a time-domain image ex-
pansion method. Frequent applications have been made
of the image method in the past as in deriving the re-
quency range of interest (i.e., 0.1-4.0 kHz), plus cor-
rections for those outside this range. The image meth-
od includes only those images contributing to the im-
pulse response. Thus the contributing images are those
within a radius given by the speed of sound times the
verberation-time equations, 3 for theoretical studies of
sound behavior in enclosures, 4'? and in the study of
architectural acoustics and perceptual properties of
reverberation time. 6 (The exact relationship between
the normal-mode solutions and the image solution, for
a lossless room, is discussed in Appendix A.) The im-
portant information used here is that in the time-do-
main, each image contributes only a pure impulse of
known strength and delay while each normal mode is
a decaying exponential which contributes to all times.
Furthermore, whereas an image has only delay and
gain as parameters, a normal mode computation re-
quires the solution of transcendental equations to find
the pole location plus the evaluation of a relatively
complex function to find the mode gain (i.e., the residue
of the pole).
rooms. 8'11 In addition, there has been a considerable
amount of important theoretical work on the approxi-
mate 12 use of images produced by a single soft-wail
(finite impedance) reflection. Several recent papers on
this subject which have good bibliographies are Refs.
13, 14, and 15. Computer methods have also recently
been applied to image computations in enclosures (see
for example Refs. 6, ?, 10, and 11). In the current
paper the computational technique is specifically aimed
at being simple, easy to use, and fast. In addition the
resulting'room responses have been used to realisticai-
1¾ mode[ speech transmission in rooms and to investi-
gate the effects of various forms of digital speech sig-
nai processors .t6, l ?
In the following we will first briefly discuss theoretic-
al aspects of the method. Then we will outline the compu-
tational approach and, finally, we will give some ex-
amples of applications.
I. IMAGE VERSUS NORMAL MODE MODELS
A. The image model
We model a talker in a room as a point source in a
rectangular cavity. A single frequency point source of
acceleration in free space emits a pressure wave of the
form
We model the rooms of interest as simple rectangular
enclosures. This choice of geometry is made for sev-
eral reasons-
P(co,X,X') exp[{o•(R/c - t)]
= 4•.R '
where
(t)
943 J. Acoust. Soc. Am. 65(4), Apr. 1979 0001-4966/79/040943-08500.80 (D1979 Acoustical Society of America 943
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
P = pressure,
IMAGE EXPANSION p(t)
•y
c• =2•rf,
f = frequency,
t = time,
X X
X X
X X
X X
0
XX
XX
X X
XX
XX XX
XX
X X
X X
XX
XX
R= Ix-x'l,
x =vector talker location (x,3,z),
(2)
XX
X' = vector microphone location (x', 3 ', z'),
i=v•/,
c =speed of so•d.
8 •o t-
p(t)=• •.
c
I 8( IRp+
p=lr=-m
When a rigid wall is present, the rigid wall (zero norm-
al velocity) boundary condition may be satisfied by
p•ci• an im•e symmetrically on the far side of the
wall. Thus,
R r=(2nLx,21Ly,2mLz)
Rp= (x- + x', y+_ y',z_ + z')
Rp-I- R r = (x_ + x'+ 2nLx,y + y'+ 21Ly,z*_zø+2rnLz )
FIG. 1. A slice through the image space showing how the
The solid box
represents the original room. The actual image space is
three dimensional.
R(•,X,X') [exp[i(•/c)R+] )
images of the source are spatially arranged.
=L 4•R + +
exp[i(•/c)R..]]exp(_{•t
4•R.. J
(3)
where we define the two distances from the microphone
to the source R. and to the im•e R. by
a z =(•- •')z + (3 -3')z + (z - z') z
' '
a z =(x+ •)z +(y -y')z +(z -z') z
+
(4)
be derived directly from the normal-mode solution as
shown in Appendix A.
The wall has been placed at x=0 in this case (note the
s•n in the x terms of R+ and R.).
In the general case of sk walls the sit•tion becomes
more complicated because each image is itself imaged.
The pressure may then be written (as shown in Appen-
dk A)
B. Case of nonrigid walls
If the room walls are not rigid, the soluHon in terms
of point images may no longer be exact. A precise
statement of the effects of finite impedance walls is
presently impossible, since the effects on even a single
image are quite complicated. la'14'15 Therefore we have
continued to assume the approximate point image model
even for nonrigid walls. In addition we have assumed
an angle independent pressure wall reflection coefficient
fl. This assumption is equivalent to assuming that the
R(•,X,X')= •
permutations over ß of
4•1• +R .I '
(6)
where • represents the eight vectors given by the eight
wall impedance is proportional to sec(e), where e is the
and
angle of incidence of a plane wave with respect to the
wall normal. We presently do not understand the exact
physical interpretation of the above assumptions. How-
ever, we believe that they do not introduce serious
problems into the final result under typical conditions.
By typical, we mean over the frequency range of 100
Hz-4. kHz, wall reflection coefficients of greater than
0.?, typical office room geometries, and where both
source and receiver are not close to the wall. Many,
•=(x•x',
•=2(.n•.
y •y', z•z')
ZL•. •n•).
r is the integer vector triplet (•,•,m),
(V)
where (Lx, L•, L•) are the room dimensions. Equation
(5) is the pressure frequency response ass•ing rigid
walls for a point source at X =(x,y,z) and receiver at
X'=(• ,3',z'). If Eq. (5) is Fourier transformed, we
find the room impulse response function (time domain
Green' s f•ction)
p(t,X ,X' )= :.- •1• + •1 '
•r,• •[t-(I R•+ • l/c)] (8)
if not all, of the above conditions are probably not cri-
tical and could be relaxed. We merely wish to carefully
point out the nonexact nature of the results.
The above assumptions result in the Sabine energy
absorption coefficient a for a uniform reflection coeffi-
cient fi on a given wall of the form
An interpretation of Eq. (8) is given in Fig. 1 where we
show a part of the image space for a two-dimensional
slice through the room. When the accelerative source
a = t - f. (O)
location (talker) X is excited, each image point is si-
multaneously excited, creating spherical pressure
waves which propagate away from each image point.
Our assumptions are similar to those of geometrical
acoustics 18 and are the same as those required for spec-
ular angle-independent ray tracing. In current imple-
mentations of the model we also do not allow frequency
variations in the reflection coefficients. Both the angle
dependence and frequency dependence could be included
944
Equation (8) is the exact solution to the wave equation
in a rectangular, rigid-wall (lossless), room and may
944 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979
J. Allen and D. Berkley: Method for simulating smalr-room acoustics
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
in our computer program, but only at the expense of
significantly complicating and slowing down the compu-
tational model.
TABLE I. Computation parameters--room size (feet) 10' x 15'
x 12.5' 8 kHz sampling rate.
Introducing the effects of finite, angle independent
Length
Impulse response
Image Computation
Convolution
rate
wall absorption into Eq. (8) leads to the modified room
impulse response
1 .•
(ms)
64
128
256
No. points
512
1024
2048
count
585
4690
37500
time (s)
i
8
60
(s/s)
12.5
13.8
15.0
•(t,X,X')=• •-• •ln'ql•lnl•ll'jl•lll•lm'kl• Irnl
•xl Px2 P•I •2 P•I •2
(10)
where • is now expressed in terms of the integer 3-
vector p = (q,j, k) as
I for our implementation on a Data General, Eclipse
•=(x-x • + 2qx', y -y' + 2jy', z-z' + 2kz') . (11)
• as given by Eq. (6) is similar to that of (11), but is
indexed differently from (11). The hera's are the pres-
sure reflection coefficients of the s• boundary planes,
with the subscript 1 referring to walls adjacent to the
S/200 computer. (On this machine the computation
time required for each image is about 1.6 ms.) The
actual FORTRAN programs used are given in Appendix B.
coordinate origin (see Fig. 1). Subscript 2 is the oppos-
ing wall. Eq. (10) has beenderiv• heuristically from
The temporal quantization in the impulse function
computation causes slight statistical errors in the com-
puted arrival times of each image pulse relative to the
geometric• consideratio• of Fig. 1. The sum • with
vector ind• p is used to indicate three sums, namely
one for each of the three components of p=(q,j,k).
r=(n,l,m) is a similar sum. Physically these s•s
are over a three-dimensional lattice of •ints. For p
there are eight points in •e lattice and for r, the lattice
is i•inite.
exact delay as given in Eq. (10). This error can be
thought of as effectively "moving" each image source by
0• ezzoz• AR/2 relative to the receiver. This effect
could be removed, in principle, by using a band-limited
source pulse. However, the error is small for most
(if not all) purposes and it greatly complicates-the com-
putation to remove this approximation. We have esti-
mated that the error due to the slight moving of the im-
ages could not be perceived even in a digital simulation
of a binaural hearing experiment.
II. IMPLEMENTATION OF THE MODEL
The primary consideration in a computer (sampled
data) implementation of Eq. (10) is the method of spatial
sampling. In addition, an apparently nonphysical be-
havior of the model at zero frequency is removed by a
The subroutine SROOM of Appendix B requires as pa-
rameters the number of impulse response points de-
sired (NPTS), the source location R0, the receiver lo-
low-frequency (0.01 of the sampling frequency) high-
pass digital filter. •9
t.0
A calculated impulse response is built up as a "histo-
gram" of image pulses received at different time de-
lays. The width of each histogram bin is equal to the
time sampling period T initially assumed, which in
turn is determined by the highest frequency to be rep-
resented. For example, all images with the range
N•R to (N + 1)AR, where •R=cT (T is the sampling
period and c the speed of sound), are added together
IMPULSE RESPONSE
with appropriate amplitude as given by Eq. (10).
The choice of sampling rate is the appli-
cation. If speech is to be studied in small rooms one
might choose T =0.1 ms. (sampling frequency of 10
kHz; highest frequency of 5 kHz). But, if reverberation
times of large enclosures are being studied (and convo-
lution with speech is not required) much lower rates
can be useful.
- 1.O
2048 POINTS
8 KHZ SAMPLING
I I I
RATE
I
The time length of the calculated impulse response is
also a consideration. For a given sampling rate the
0 TIME (MS) 256
number of points in p(t) increases linearly with its
length while the computation time (and number of im-
ages) goes up approximately as the cube of response
length. This is shown in the first four columns of Table
945 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979
FIG. 2. Plot of a typical impulse response for a room
80 x120 x100 sample lengths long. Wall reflection coefficients
were all 0.9 ceiling and floor coefficients were 0.7. X and X'
were at (30,100, 40) and (50, 10, 60) sample periods.
J. Allen and D. Berkley: Method for simulating small-room acoustics 945
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
cation R, the room dimensions, all specified in terms
use a modification of the integrated tone-burst method 22
of the sample length (AR), and the reflection coefficients
of each of the six wall surfaces (•). Figure 2 shows an
example of the impulse response obtained for a room of
dimensions 80 x 120 x 100 sample lengths with equal wall
reflection coefficients of 0.9 (c• =0.19) and with floor
and ceiling reflection coefficients (•) of 0.7 (a -- 0.51). X
and X' were (30,100,40)and (50,10,60)sample
tervals, respectively.
where E(t) is the average energy decay, k is a propor-
tionality constant, and p(•) is the calculated pressure
impulse response from Eq. (10). For cases where the
impulse response has been truncated before most of the
in-
It is usually convenient to interpret the model param-
eters as a true distance rather than as multiples of
AR. This requires the choice of a sampling rate and
then conversions may be performed in the users main
program which calls the subroutines of Appendix B.
Figure 2 is labeled assuming an 8 kHz sampling rate.
decay has taken place, (12) may lead to errors. These
errors are usually obvious in the E(t) plots.
Another, approximate, approach is to simply mea-
sure the short-time average energy decay of the im-
pulse itself (e.g., using a simulated level recorder).
For exponential or near-exponential decays, both
methods should give approximately the same value of
For this assumption (and assuming a sound speed of 1
ft/ms) the room dimensions are 10' x 15' x12.5'.
III. APPLICATIONS
Our room image model has been applied to several
physical evaluation of room reverberation effects 1'2 and
a study of critical distance measurements using spec-
problems. We will discuss two examples: a psycho-
(o)
tral response variance. zø We have also used the model
to test a signal processor intended to reduce perceived
reverberationl• 6 and to study problems associated with
mathematical inversion (inverse filtering) of room
transfer functions.• ?
A. Psychophysics of room reverberation
Once a simulated room impulse response has been
calculated using the image model, the psychophysical
effects of this simulated reverberation on speech may
be directly studied. A reverberant sample of speech
DECAY CURVE
2048 POINTS
8 KHZ SAMPLING RATE
was produced by convolving an anechoic (unreverberant)
speech sample with the calculated impulse response
[P(t)l. This can be done efficiently using a Fast Fourier
Transform (FFT) method (overlap-add) to perform the
convolution. zi The last column of Table I shows the
measured convolution rate, for various length impulse
responses. The convolution rate only increases as
TIME (MS)
256
logz(N), (where N is the room response length in time
sample periods T) so even large impulse responses can
be convolved with speech quite efficiently. For exam-
ple, to convolve (filter) one second of speech, sampled
at 8 kHz, with a 256 ms long impulse response (2048
points) requires a 15 s computation.
ENERGY DECAY CURVE
The speed of processing makes multivariate psycho-
physical studies quite practical. Ease of modification
and perfect control of room parameters avoids the
problems which have made such experiments so difficult in
the past. The actual experiments used 16 different simu-
lated" rooms" (impulse res ponses) convo lved with ten dif-
ferent sentences spokenby four different speakers. Itwas
discovered that the experimental rooms were perceptually
8 KHZ. SAMPLING
I ,
RATE
, I, , I ,,
256
well characterized by their spectral variance [Eq. (14)]
andby the reverberation time. This latter measure,
reverberation time, deserves some discussion.
, i , i •l, ,, I ..... 1
o
TIME (MS)
Given the impulse responses, reverberation time may
be estimated in a number of ways. One method is to
946 J. Acoust. $oc. Am., Vol. 65, No. 4, April 1979
Fig. 2 using the Schroeder integration method. 22 (b) Impulse
energy decay curve for a simulated level recorder.
946
FIG. 3. (a) Energy decay curve for the impulse response of
J. Allen and D. Berkley' Method for simulating small-room acoustics
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
reverberation time. Example plots of E(t) for both
methods are shown in Figs. 3(a) and 3(b) using the im-
pulse response of Fig. 2. Experience indicates thatEq.
(12) gives the most satisfactory results. We have found
empirically that calculated reverberation times, for a
number of simulated enclosures, agree well with
Eyring's formula s over a wide range of Beta values. 23
In the experiments discussed above 1' •-4 we discovered a
monotonic relationship between A/A c( Fig. 4)the micro-
(a) STANDARD DEVIATION OF SPECTRAL RESPONSE
IN COMPUTER-SIMULATED ROOM: 17' x !$' x
7
phone-talker distance when normalized by the room c ri-
tical distance (the distance at which reverberant energy
equals direct sound energy), the reverberation time, and
psychophysical preference for the resulting speech.
B. Critical distance measurement
z
o
n- 5
. A new method has been proposed •'ø for measurement
of critical distances (or reverberation radius) in rooms.
In this technique a measurement is made of the log frequen-
cy response variance gL defined as
u•4
L(w)=20 logtiP(w) [ ] (13)
(14)
z
o
z
given the room transfer-function P(•o) [Fourier trans-
form of Eq. (10)] for several microphone-source spac-
ings. The measured values are fitted to a theoretical
curve for (rL based on the assumption of simultaneously
excited, uncorrelated, normal modes, combined with
the calculated direct sound energy. The resulting fit
was shown to give an accurate value for the room's
critical distance.
o
z
This new method was extensively studied using our
0
-$0 -24
I
-18
I
-12
I
-6
i
0
I
6
[
12:
I
18
I
image model before being applied successfully to real
24
REVERBERANT/DIRECT ENERGY RATIO (dB)
rooms. Since the direct and reverberant energy are
(•/•c)'
(•) STANDARD DEVIATION OF SPECTRAL RESPONSE
I1• COMPUTER-SIMULATED ROOM; 47' x :51' x 15'
7•
known in the computer model, a comparison can easily
be made to the theory. The model results show excell-
ent agreement with theoretical calculation as is seen in
Figs. 4(a) and 4(b). We know of no other method by
which this study could have been carried out as effect-
ively.
• 6- o
•
f•
0
&•l•• • 0
p ß o.•
.oo
ß
ß o n
IV. SUMMARY
A simulation
AND DISCUSSION
method for small rooms based on an
approximate image expansion for rectangular nonrigid-
wall enclosures has been discussed. The method is
simple, easy to implement and efficient for computer
simulation. Several examples of its use, where other
methods would be difficult, have been discussed.
APPENDIX A
o _
-
>
We wish to derive the rigid-wall image solution di-
rectly from the normal-mode expansion for a rectangu-
lar enclosure. The frequency response function
(Green's function) for the pressure P(•o) in an enclos-
ure is given by solving the Helmholtz equation driven
by a single frequency point acceleration source.
I
t2
0 I I
-18
I
-12
, I
-6
I
0
I
6
I
18 24 -E•) -24
REVERBERANT/DIRECT ENERGY RATIO (dB)
,x'] + ,x'] =- ,(x -x '),
(A1)
(A/Ac)2
FIG. 4. Figures from Jetzt 2ø which compare the theoretical
rms deviation of the pressure in dB from the mean pressure
in dB as a function of the direct to reverberent energy ratio
(a) for a room 17x13x10 ft and (b) 47x31 x15 ft.
where • is the frequency and c is the speed of sound,
The solution to this equation, assuming rigid boundar-
ies, is given by
947 947 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics
Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17
P(k,X,X')=• =.• (kZr_k z) ,
•,• •,(x)%(x' )
where k=o•/c, r=(n,l,m)
al sum, V is the room volume,
(A2)
By Fourier series analysis one may show
indicates a three dimension-
-,, ß
• ( n•)-L--•exp(i2L n•x)
(A9)
Thus [with analogus equations to (A9) for y and z ]
, , LE '
kZr= Ik, I2
and
(A3)
f
• ,=-•
• exPt(i•'(Rv•z R•)] d•
(•z_ ) ,
(A10)
where • is the vector [also given by Eq. (7)]
/n,rXcos[l•ry /m•rz
]
•=2(•L,, •L•, mL•). (A•)
where the L•'s are the room dimensions.
Each triple integral is just a plane wave expansion for
a •int source in free space since
Using the exponential expansion for cosine, multi-
plyi• the terms of Eq. (A2) together and collecting, we
obtain
P(k,X,X') =•,•
•q• (•)]
1 • •exp(•-•)
• (•-•:)
exp(ikl RI ) ='•'x• (I tj I z - k •)
4•1RI '
(A12)
,
(A5)
(A6)
Finally, using Eq. (A12), Eq. (A10) becomes
4,tr I R• + l:{rl
p(,• ,X,X,)=••__ r__•.•exp[i(•/c)l Rp + R•I]
(A13)
where • represents the eight vectors [also given by
• =(x •x', y •y', • •').
Taking the inverse Fourier transform ofEq. (A13), the
echo structure becomes explicit
Using the property of the delta Mnction onk,, k•, and k,
f" •(• - a)r(•)d• =r(a),
we may rewrite Eq. (A5) in integr• form
(A7)
(AO)
1 •
APPENDIX
C PGM:
fexp(if-•) •
, (A14)
=_• 4•' 11• + R,.I
p(t,X,X')=••__ r•'• 6[t-(lRp+Rr[/c)]
which is the same as Eq. (8) as desired.
B
SROOM
C SUBROUTINE TO CALCULATE A ROOM IMPULSE RESPONSE
C R=VECTOR RADIUS TO RECEIVER IN SAMPLE PERIODS=LENGTH/(C*T)
C R0 =VECTOR RADIUS TO SOURCE IN SAMPLE PERIODS
C RL=VECTOR OF BOX DIMENSIONS IN SAMPLE PERIODS
C BETA=VECTOR OF SIX WALL REFLECTION COEFS (0 C HT =IMPULSE RESP ARRAY C N-PTS=# OF POINTS OF HT TO BE COMPUTED C ZERO DELAY IS IN HT(1) C SUBROUTINE SROOM (R, R0, RL, BETA, HT, N-PTS) DIMENSION HT (N-PTS) DIMENSION R(3), R0(3), NR(3), RL(3), DELP(8), BETA(2,3) EQUIVALENCE (NR(1), NX), (NR(2), NY), (NR(3), NZ) DO 5 I= 1, N-PTS 5 HT(I) =0 C CK FOR MIC AND SOURCE AT SAME LOCATION DIS =0 6 DO 6 I=1,3 DIS = (R(I)-R0 (I))*'2 +DIS DIS = SQRT(DIS) IF (..5)HT(1) = 1 IF (..5)RETURN RANGE OF SUM C FIND N1 = N-PTS/(RL(1)* 2) + 1 N2 = N-PTS/(RL(2)* 2) + 1 N3 = N-PTS/(RL(3)*2) + 1 DO 20 NX=-N1, N1 DO 20 NY =--N2, N2 DO 20 NZ =--N3, N3 948 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics 948 Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17 C GET EIGHT IMAGE LOCATIONS IO=O FOR MODE # NR CALL LTHIMAGE (R, R0, RL, NR, DELP) DO 10 L = 0,1 DO 10 J =0, 1 DO 10 K = O, 1 IO=IO + 1 C MAKE DELAY =ID AN INTEGER ID = DELP (I0) + .5 FDM1 ID =ID +1 IF ()GO C PUT IN LOSS FACTOR TO 10 ONCE FOR EACH WALL REFLECTION Gm= BE TA(1,1)** IABS(NX-L) 1 *BETA(2,1)**IABS(NX) 2 *BETA (1,2)** IABS(NY-J) 3 *BETA(2,2)**IABS(NY) 4 *BE TA (1,3**IABS (NZ-K) 5 *BETA(2,3)**IABS(NZ) 6 /FDM1 C CHECK FOR FLOATING POINT UNDERFLOW HERE; C IF UNDER FLOW, SKIP NEXT LINE HT(ID) = HT(ID) +GID 10 20 CONTINUE CONTINUE NOTE LINES CONTINUATION C IMPULSE RESP HAS BEEN COMPUTED C FILTER WITH HI PASS FILT OF 1% OF SAMPLING FREQ (I.E'. 100 HZ) C IF THIS STEP IS NOT DESIRED, W = 2.*4.*ATAN(1.)* 100. T = 1E-4 RETURN HERE R1 = EXP(-W*T) R2 =R1 B1 = 2.* RI* COS(W* T) B2 =-RI*R1 A1 =-(1. +R2) A2 =R2 Y1 =0 Y2 =0 Y0 =0 C FILTER HT DO 40 I=1, NPTS X0 =HT(I) HT(I) =Y0 +AI*YI+A2*Y2 Y2 =Y1 Y1 =Y0 Y0 =BI*Y1 40 CONTINUE RETURN END +B2*Y2 +X0 C PGM: C PGM C LTHIMAGE TO COMPUTE EIGHT IMAGES OF A POINT IN BOX SUBROUTINE LTHIMAGE (DR, DR0, RL, NR, DELP) C C DR IS VECTOR C DR0 IS VECTOR C RL IS VECTOR C NR IS VECTOR C DELP C C RADIUS RADIUS TO RECEIVER TO SOURCE IMAGE IN SAMPLE IN SAMPLE IN SAMPLE TO IMAGE PERIODS PERIODS PERIODS OF BOX DIMENSIONS OF MEAN IN SAMPLE NUMBER IS VECTOR OF EIGHT SOURCE PERIODS DISTANCES DIMENSION R2L(3), RL(3), NR(3),DELP(8) DIMENSION DR0(3), DR(3), RP(3,8) C LOOP OVER ALL SIGN PERMUTATIONS AND COMPUTE R +/- R0 I0= 1 DO 10 L =--l,1,2 DO 10 J =--1,1,2 DO 10 K=-l,1,2 949 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics 949 Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17 C NEAREST IMAGE IS L = J = K =-1 RP (1, I0)=DR(1)+ L'DR0(1) RP(2,I0) =DR(2) +J. DR0(2) RP(3, I0) =DR(3) + K. DR0(3) I0 = I0 + 1 10 CONTINUE C ADD IN MEAN RADIUS TO EIGHT VECTORS TO GET TOTAL DELAY R2L (1) = 2.* RL (1)* NR(1) R2L (2) = 2.*RL (2)* NR(2) R2L ($) = 2.* RL (3)* NR(3) DO 20 I=1,8 DELSQ=0 DO 25 J = 1,3 R1 = R2L (J)-RP(J,I) DELSQ =DELSQ +Rl**2 25 CONTINUE DELP (I) = SQRT(DELSQ) 20 CONTINUE RETURN END 1Barbara McDermott and Jont Allen, "Perceptual Factors of Small Room Reverberation," J. Acoust. Soc. Am. 60, S9 (A) (1976). an Impedance Boundary," J. Acoust. Soc. Am. 59, 780-785 (1976). 16R. J. Donato, '•Spherical Wave Reflections from a Boundary 2jont B. Allen and B. J. McDermott, "A New Method for Mea- suring Perception of Room Reverberation" (unpublished). 3C. F. Eyring, "Reverberation Time in 'Dead' Rooms," J. Acoust. Soe. Am. 1, 217-241 (1930). 4D. Mintzer, "Transient Sounds in Rooms," J. Acoust. Soc. 51{. H. Bolt, P. E. Doak, and P. J. Westervelt, "Pulse Sta- tistics Analysis of Room Acoustics," J. Acoust. Soc. Am. Am. 22, 341-352 (1950). 22, 328-339 (1950). 16Jont B. Allen, D. A. Berkley, and J. Bauert, '•lViultimicro- phone signal-processor technique to remove room reverbera- tion from Speech Signals," J. Acoust. SOc. Am. 62, 912-915 (1977). of Reactive Impedance using a Modification of Cagniard's Methods," J. Acoust. Soc. Am. 60, 999-1002 (1976). ?S. T. Neeley and Jont B. Allen, "Invertibility of a Room Im- pulse Response," J. Acoust. Soc. Am. (in press). 18p. M. Morse and K. U. Ingard, Theoretical Acoustics (McGraw-Hill, New York, 1968), Sec. 9.5. 6j. M. Berman, "Behavior of Sound in a Bounded Space," J. Acoust. Soc. Am. 5?, 1275-1291 (1975). 1E. K. Dunens and R. F. Lambert, "Impulsive Sound Response Statistics," J. Acoust. Soc. Am. 61, 1524-1532 (1977). 6J. R. Power, "Measurement of Absorption in Rooms with Sound Absorbing Ceilings," J. Acoust. Soc. Am. 10, 98-101 (1938). SC. G. Mayo, '•tanding Wave Patterns in Studio Acoustics," Acustica 2, 49-64 (1951). lSWithout the high-pass filter, the Fourier transform ofp(t) shows a large energy spike at zero frequency as the absorp- tion a goes to zero. This is in part a direct result of the no- pressure-release boundary conditions we have assumed and is •urther complicated by the assumption of an accelerative source, which is nonphysical for very low frequencies (less than 50 Hz). 1øB. M. Gibbs and D. K. Jones, "A Simple Image Method for Calculating the Distribution of Sound Pressure Levels within an Enclosure," Acustica 26, 24-32 (1972). 20j. j. Jetzt, "Critical Distance Measurement of Rooms from the Sound Energy Spectral Response," J. Acoust. Soc. Am. (to be published). llF. Santon, "Numerical Prediction of Echograms and the In- telligibflity of Speech in Rooms," J. Acoust. Soc. Am. 59, 1399-1405 (1976). 12By "approximate" we mean the blurring of the image as a re- sult of the finite impedance wall. Variations in wall imped- ance as a function of frequency or other frequency dispersion effects will also make the images frequency dependent. In practice the two effects usually appear simultaneously. 21Jont B. Alien, "FASTFILT-An FFT Based Filtering Pro- gram, "in IEEE Press book on Programs for Digital Signal Pro- cessing (IEEE Press, New York, 1979). 22M. R. Schroeder, "New Method of Measuring Reverberation Time," J. Acoust. SOc. Am. 37, 409-412 (1965). 2SAgreement is best for beta the same on all surfaces and variations appear when a pair of opposing walls are signifi- canfly different in reflectivity than all other surfaces, (as is 13A. R. Wentzel, "Propagation of Waves Along an Impedance Boundary," J. Acoust. Soc. Am. 55, 956-963 (1974). 14S. Thomasson, '•eflection of Waves from a Point Source by 24D. A. Berkley, "Normal Listeners in Typical Rooms:Rever- beration Perception, Simulation and Reduction" (unpublished). the case for Figs. 2 and 3). 950 J. Acoust. Soc. Am., Vol. 65, No. 4, April 1979 J. Allen and D. Berkley: Method for simulating small-room acoustics 950 Redistribution subject to ASA license or copyright; see /content/terms. Download to IP: 202.118.235.3 On: Sun, 13 Nov 2016 11:43:17