Estimation of Modal Decay Parameters from Noisy Response Measurements
Matti Karjalainen1 , Poju Antsalo1 , Aki Mäkivirta2 , Timo Peltonen3 , and Vesa Välimäki1
1Helsinki University of Technology, Laboratory of Acoustics and Audio Signal Processing,
P.O. Box 3000, FIN-02015 HUT, Finland
2

Genelec Oy, Iisalmi, Finland, and 3 Akukon Oy, Helsinki, Finland

ABSTRACT
Estimation of modal decay parameters from noisy measurements of reverberant and resonating systems is a common
problem in audio and acoustics, e.g., in room and concert hall measurements or musical instrument modeling. In thispaper, reliable methods to estimate the initial response level, decay rate, and noise ﬂoor level from noisy measurement
data are studied and compared. A new method, based on nonlinear optimization of a model for exponential decay
plus stationary noise ﬂoor, is presented. Comparison with traditional decay parameter estimation techniques using
simulated measurement data shows that the proposed method outperforms in accuracy and robustness, especially inextreme SNR conditions. Three cases of practical applications of the method are demonstrated.
0

INTRODUCTION

Parametric analysis, modeling, and equalization (inverse
modeling) of reverberant and resonating systems ﬁnd many
applications in audio and acoustics. These include room andconcert hall acoustics, resonators in musical instruments, andresonant behavior in audio reproduction systems. Estimation of reverberation time or modal decay rate are importantmeasurement problems in room and concert hall acoustics [1],where S/N ratios of only 30-50 dB are common. The sameproblems can be found for example in the estimation of parameters in model-based sound synthesis of musical instru-

ments, such as vibrating strings or body modes of string instruments [2]. Reliable methods to estimate parameters fromnoisy measurements are thus needed.
In an ideal case of modal behavior, after a possible initialtransient, the decay is exponential until a steady state noise
ﬂoor is encountered. The parameters of primary interest tobe estimated are:
- Initial level of decay (LI )- Decay rate/time or reverberation time (TD )- Noise ﬂoor level (LN )In a more complex case there can be two or more modal fre-

KARJALAINEN ET AL.
quencies, whereby the decay is not simple anymore, but shows
additional ﬂuctuation (beating) or a two-stage (or multiplestage) decay behavior. In a diﬀuse ﬁeld (room acoustics) the
decay of a noise-like response is approximately exponential inrooms with compact geometry. The noise ﬂoor may also be
non-stationary. In this article we primarily discuss a simple
mode (i.e., a complex conjugate pole pair in transfer function)
or a dense set of modes with exponential reverberant decay,
together with a stationary noise ﬂoor.
Methods presented in literature and common-sense or ad-hoc
methods will ﬁrst be reviewed. Techniques based on energytime curve analysis of signal envelope are known as methods
where the noise ﬂoor can be found and estimated explicitly.
Backwards integration of energy, so called Schroeder integration [3]-[4], is often applied ﬁrst to obtain a smoothed envelope for decay rate estimation. AR modeling of modes by
estimating the transfer function poles and group delay analysis are examples of straightforward methods which are notparticularly robust against background noise.
The eﬀect of background noise ﬂoor is known to be problematic, and techniques have been developed to compensate theeﬀect of envelope ﬂattening when the noise ﬂoor in a measured response is reached, including limiting the period of
integration [5], subtracting an estimated noise ﬂoor energy
level from a response [6], or using two separate measurements
to reduce the eﬀect of noise [7]. The iterative method by Lundeby et al. [8] is of particular interest since it addresses withcare the case of noisy data. This technique, as most othermethods, analyzes the initial level LI , decay time TD , andnoise ﬂoor LN parameters separately, typically starting froma noise ﬂoor estimate. Iterative procedures are common in
accurate estimation.
A diﬀerent approach was taken by Xiang [9] where aparametrized signal-plus-noise model is ﬁtted to Schroederintegrated measurement data by searching for a least squares(LS) optimal solution. In this study we have elaborated asimilar method of nonlinear LS optimization further to makeit applicable to a wide range of situations, showing good convergence properties. A speciﬁc parameter and/or a weighting function can be used to further ﬁne-tune the method forspeciﬁc problems. The technique is compared with the Lundeby et al. method by applying them to simulated cases of
exponential decay plus stationary noise ﬂoor where the exact
parameters are known. The improved nonlinear optimization
technique is found to outperform traditional methods in accuracy and robustness, particularly in diﬃcult conditions withextreme signal-to-noise ratios.
Finally, the applicability of the improved method is demonstrated by three examples of real measurement data: (a) reverberation time of a concert hall, (b) low-frequency mode
analysis of a room, and (c) parametric analysis of guitar string
behavior for model-based sound synthesis. Possibilities forfurther generalization of the technique to more complex problems, such as two-stage decay, will be discussed brieﬂy.
1

DEFINITION OF PROBLEM DOMAIN

A typical property of resonant acoustic systems is that theirimpulse response is a decaying function after a possible initialdelay and the onset. In the simplest case the response of a

ESTIMATION OF MODAL DECAY PARAMETERS
single mode resonator system ish(t) = Ae−τ (t−t0 ) sin[ω(t − t0 ) + φ]u(t − t0 )(1)

where u(t − t0 ) is a step function with value 1 for t ≥ t0
and zero elsewhere, A is initial response level, t0 response
latency for example due to propagation delay of sound, τdecay rate parameter, ω = 2πf angular frequency, and φ initial phase of sinusoidal response. In practical measurements,
when there are multiple modes in the system and noise (acoustic noise plus measurement system noise), a measured impulse
response is of the form1
h(t) =NX

Ai e−τi (t−t0 ) sin[ωi (t − t0 ) + φi ] + An n(t)

(2)

i=1

where An is the rms value of background noise and n(t) is
unity level noise signal. Figure 1 illustrates a single delayedmode response corrupted by additive noise.

1

delay

(a)

exponential decay

0.5

noise floor level

0
-0.5
-1
dB
initial level (IL)

0

(b)

-5
-10
-15

LN

noise floor level (LN)

-20
-25
-30
dB
initial level
level (IL)(IL)
initial

0

(c)
(b)

-5
-10
-15

LN
L
N

noise floor
floor level
level (L
(LN))
noise
N

-20
-25
-30
0

200

400

600

800

1000

1200 time in samples

Fig. 1: (a) Single mode impulse response (sinusoidal decay)
with initial delay and additive measurement noise. (b) Absolute value of response on the dB scale to illustrate the decay
envelope. (c) Hilbert envelope, otherwise same as (b).
The problem of this study is deﬁned as a task to ﬁnd reliable estimates for the parameter set {Ai , τi , t0 , ωi , φi , An },given a noisy measured impulse response of the form (2). The
main interest here is focused on systems of (a) separable single modes of type (1) including additive noise ﬂoor or (b)dense (diﬀuse, noiselike) set of modes resulting also in exponential decay similar to Fig. 1. In both cases the parametersof primary interest are A, τ, An , and t0 .

1
In a more general case the initial delays may diﬀer and there can be simple non-modal exponential terms, but these cases
are of less importance here.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

2

KARJALAINEN ET AL.

ESTIMATION OF MODAL DECAY PARAMETERS

Often the decay time TD is of main interest, for example inroom acoustics where the reverberation time2 of 60 dB decay,
T60 , is related to τ by
T60 = −

1
6.908
ln(10−3) ≈
τ
τ

(3)

Modern measurement and analysis techniques of system response are carried out by digital signal processing whereby
discrete-time formulation for modal decay (without initial delay) with sampling rate fs and sample period Ts = 1/fs becomes(4)
h(n) = Ae−τd n sin(Ωn + φ)
where n is sample index, τd = Ts τ , and Ω = 2πTs f .
2

Line ﬁtting (linear regression to log envelope)Fitting a line in a logarithmic decay curve is a conceptuallyand computationally simple way of decay rate estimation.
The decay envelope y(t) can be computed simply as a dBscaled energy-time curvey(t) = 10 lg{x2 (t)}

Zt

min

DECAY PARAMETER ESTIMATIONIn this section an overview of known techniques for decay parameter estimation will be presented. Initial delay and level
estimation is discussed ﬁrst brieﬂy. The main problem, i.e.,
decay rate estimation, is the second topic. Methods to smooththe decay envelope from a measured impulse response are presented shortly. Noise ﬂoor estimation, an important subproblem, is discussed next. Finally, techniques for combined noise
ﬂoor and decay rate estimation are reviewed.
2.1

Initial delay and initial level estimation

In most cases these two parameters are relatively easy to estimate. The initial delay may be short, not needing any attention, or the initial bulk delay can be cut oﬀ easily up to
the edge of response onset. Only when the onset is relatively
irregular or the S/N ratio is low, the detection of onset time
moment can be diﬃcult.
A simple technique to eliminate initial delay is to computethe minimum-phase component hmphase (t) of the measured
response [10]. An impulse response can be decomposed
as a sum of minimum-phase and excess-phase components:
h(t) = hmphase (t) + hephase (t). Since the excess-phase component will have allpass properties manifested as delay, computation of the minimum-phase part will remove the initialdelay.
The initial level in the beginning of the decay can be detected
directly from the peak value of the onset. For improved robustness, however, it may be better to estimate it from thematched decay curve, particularly its value at the onset time
moment.
In the case of a room impulse response, the onset corresponds
to direct sound from the sound source. It may be of special
interest for the computation of source-to-receiver distance orin estimating the impulse response of the sound source itself
by windowing the response before the ﬁrst room reﬂection.
2.2

Decay rate estimation

Decay rate or time estimation is in practice based on ﬁtting
a line to the decay envelope, such as the energy-time curve,
mapped on a logarithmic (dB) scale. Before computerized age
this was done graphically on paper. The advantage of manual inspection is that an expert can avoid data interpretation
errors in pathological cases. However, automatic determination of the decay rate or time is, however, highly desirable inpractice.
2(5)where x(t) is the measured impulse response or a band-passﬁltered part of it, such as an octave or 1/3-octave band. It
is common to apply techniques such as Schroeder integrationand Hilbert envelope computation (to be described below) inorder to smooth the decay curve before line ﬁtting. Leastsquares line ﬁtting (linear regression) is done by ﬁnding the
optimal decay rate k ask,a2

{y(t) − [k(t) + a]}2 dt

(6)

t1

for example using the Matlab function polyfit [11].
Practical problems with line ﬁtting are related to the selectionof interval [t1 , t2 ] and cases where the decay of the measured
response is inherently nonlinear. The ﬁrst problem is avoided
by excluding onset transients in the beginning and noise ﬂoor
biasing at the end of interval [t1 , t2 ]. The second problem isrelated to such cases as a two-stage decay (initial decay rateor early reverberation and late decay rate or reverberation)or beating (ﬂuctuation) of the envelope because of two modes
very close in frequency (see Fig. 9b).

Nonlinear regression (Xiang’s method)
Xiang [9] formulated a method where a measured andSchroeder-integrated energy-time curve is ﬁtted to a parametric model of a linear decay plus a constant noise ﬂoor.
Since the model is not linear in parameters, nonlinear curve
ﬁtting (nonlinear regression) is needed. Mathematically thisis done by iterative means such as starting from a set of initial values for the model parameters and applying gradientdescent to search for a least squares optimum

Zt n
2

min

x1 ,x2 ,x3ysch (t) − [x1 e−x2 t + x3 (L − t)]

o2

dt

(7)

t1

where ysch (t) is the Schroeder-integrated energy envelope, x1
is the initial level, x2 the decay rate parameter, x3 a noise
ﬂoor related parameter, L the length of response, and (t1 , t2 )the time interval of nonlinear regression. Notice that the lastterm of noise ﬂoor eﬀect is a descending line instead of a constant due to backward integration of the noise energy [9].
Nonlinear optimization is mathematically more complex thanlinear ﬁtting and care should be taken to guarantee convergence. Even with convergence, the result may be only a local
optimum, and generally the only way to know that a global
optimum is found is to apply exhaustive search over possible
value combinations of model parameters which, in a multiparameter case, is often computationally too expensive.
Nonlinear optimization techniques will be studied in more detail later in this paper by introducing generalizations to themethod of Xiang and by comparing the performance of different techniques in decay parameter estimation.

Autoregressive modeling (linear prediction)
For a single mode of Eq. (1) the response can be modeled as animpulse response of a resonating second-order all-pole ﬁlter.

See reverberation time deﬁnition notes on page 10.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

3

KARJALAINEN ET AL.

ESTIMATION OF MODAL DECAY PARAMETERS

More generally, a combination of N modes can be modeled
as a 2N -order all-pole ﬁlter. Auto-regressive (AR) modeling [12] is a way to derive parameters for such a model. In
many technical applications this method is called linear prediction [13]. For example the function lpc in Matlab [14]processes a signal frame through autocorrelation coeﬃcientcomputation and solving the normal equations by Levinson
recursion, resulting in a N th order z-domain transfer func−i ). Poles are obtained by solving thetion 1/(1 + N
i=1 αi z
roots of the denominator polynomial. Each modal resonanceappears as a complex-conjugate pole pair (zi , zi∗ ) in the complex z-plane with angle φ = arg(zi ) = 2πf /fs and radiusr = |zi | = e−τ /fs , where f is modal frequency, fs is sampling
rate, and τ is the decay parameter of the mode in Eq. (1).

P

Decay parameter analysis by AR modeling is not robust withlong decay times (poles close to the unit circle) and background noise. Additive noise ﬂattens the power spectrumof a modal resonance and moves the poles from unit circle
towards the origin, thus resulting in shortened decay time estimates, which is contrary to the result in linear curve ﬁtting.
AR modeling works for decay parameter estimation only if it
can spectrally resolve each mode. Thus it is not applicableto high modal densities, such as typical reverberation timemeasurements.

Decay envelope smoothing techniques

In the methods of linear or nonlinear curve ﬁtting it is desirable to obtain a smooth decay envelope before the ﬁtting
operation. The following techniques are often used to improve
the regularity of the decay ramp.

Hilbert envelope computation
In this method, signal x(t) is ﬁrst converted to an analytic signal xa (t) so that x(t) is the real part of xa (t) and the complex
part of xa (t) is the Hilbert transform (90o phase shift) [10] ofx(t). For a single sinusoid this results in an entirely smoothenergy-time envelope. An example of Hilbert envelope for anoisy modal response is shown in Fig. 1c.

Schroeder integration
A monotonic and smoothed decay curve can be produced by‘backward integration’ of the impulse response h(t) over measurement interval [0,T ] and converting it to a logarithmic scale

R T h2 (τ )dτ !
RtT h2 (τ )dτ [dB](8)0This process is commonly known as the Schroeder integration [3, 4]. Based on its superior smoothing properties it isused routinely in modern reverberation time measurements.
3

0.5
(a)
0
-0.5
-1
0
dB
0

200

400

600

800

1000

(e)

(f)

800

1000

1200 time in samples

-5
-10
-15

(b)
(c)

-25
0
dB
0

200

(d)

(b)

-20

400

600

(g)
1200 time in samples

-5

-15

A complementary method to AR modeling is to use the groupdelay, i.e., phase derivative Tg (ω) = −dϕ(ω)/dω, as an estimate of the decay time for separable modes of an impulse
response. While AR modeling is sensitive to power spectrum
only, the group delay is based on phase properties only. For aminimum-phase single mode response the group delay at themodal frequency is inversely proportional to the decay parameter, i.e., Tg = 1/τ . Group delay computation is somewhat
critical due to phase unwrapping needed, and the method is
sensitive to measurement noise.

L(t) = 10 lg

1-10Group delay analysis

2.3

A known problem with it is that if the background noise ﬂoor
is included within the integration interval, the process produces a raised ramp that biases upwards the late part of decay. This is shown in Fig. 2 for the case of noisy single-mode
decay of (a) for full response integration shown as curve (d).

-20

(h)

-25
0

200

400

600

800

1000

1200 time in samples

Fig. 2: Results of Schroeder integration applied to noisy decay of a mode: (a) measured noisy response including initial delay, (b) true decay of noiseless mode (dashed straightline), (c) noise ﬂoor (-26 dB), (d) Schroeder integration of
total measured interval, (e) integration over a short interval(0,900), (f) integration over interval (0,1100), (g) integration
after subtracting noise ﬂoor from energy-time curve, and (h)a few decay curves integrated by the Hirata’s method.
The tail problem of Schroeder integration has been addressedby many authors, for example in [15, 8, 5, 6], and techniques
to reduce slope biasing have been proposed. In order to apply
these improvements, a good estimate of the noise ﬂoor levelis needed ﬁrst.
2.4

Noise ﬂoor level estimation

The limited signal-to-noise ratio inherent in practically all
acoustical measurements, and especially measurements performed under ﬁeld conditions, call for attention concerning
the upper time limit of decay curve ﬁtting or Schroeder integration. Theoretically this limit is set to inﬁnity, but inpractical measurements it is naturally limited to the lengthof the measured impulse response data. In practice, measured
impulse responses must be long enough to accommodate forlarge enough dynamic range or the whole system decay down
to the background noise level3.
Thus the measured impulse response typically contains not
only the decay curve under analysis, but also a steady level ofbackground noise, which dominates for some time at the endof the response. Fitting the decay line over this part of envelope or Schroeder integrating this steady energy level along

This is needed to avoid time aliasing with MLS and other cyclic impulse response measurement methods.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

4

KARJALAINEN ET AL.
with the exponential decay curve causes an error both in the
resulting decay rate (see Fig. 2) and in the time-windowed
energies (i.e., energy parameters).
To avoid bias by noise, analysis must be performed on the impulse response data to ﬁnd the level of background noise andthe point where the room decay meets the noise level. Thisway it is possible to eﬀectively truncate the impulse responseat the noise level, minimizing the noise energy mixed with theactual decay.
The determination of noise ﬂoor level is diﬃcult without using
iterative techniques. The method by Lundeby et al. that will
be outlined below is a good example of iterative techniques,
integrated with decay rate estimation.
A simple way to have a reasonable estimate of backgroundnoise ﬂoor is to average a selected part of the measured response tail or to ﬁt a regression line to it [16]. The level iscertainly overestimated if the noise ﬂoor is not reached, but
this is not necessarily problematic contrary to underestimating it. Another technique is to look at the background level
before the onset of main response. This works if there isenough initial latency in the system response under study.
2.5

Decay estimation with noise ﬂoor reduction

In addition to determining the response starting point, it isthus essential to ﬁnd an end point where the decay curve
meets background noise, and to truncate the noise from theend of the response. Figure 2 illustrates the eﬀect of limiting
the Schroeder integration interval. If the interval is too short
as in (e), the curve is biased downwards. Curve (f) shows acase where the bias due to noise is minimized by considering
the decay only downto 10 dB above the noise ﬂoor.

ESTIMATION OF MODAL DECAY PARAMETERS
2. A ﬁrst estimate for the background noise level is determined from a time segment containing the last 10 %of the impulse response. This gives a reasonable statistical selection without a large systematic error, if the
decay continues to the end of the response.
3. The decay slope is estimated using linear regression between the time interval containing the response 0 dB
peak, and the ﬁrst interval 5–10 dB above the background noise level.
4. A preliminary crosspoint is determined at the intersection of the decay slope and the background noise level.
5. A new time interval length is calculated according to
the calculated slope, so that there are 3–10 intervals
per 10 dB of decay.
6. The squared impulse is averaged into the new local timeintervals.
7. The background noise level is determined again. The
evaluated noise segment should start from a point corresponding to 5–10 dB of decay after the crosspoint, ora minimum of 10 % of the total response length.
8. The late decay slope is estimated for a dynamic rangeof 10–20 dB, starting from a point 5–10 dB above thenoise level.
9. A new crosspoint is found.
Steps 7–9 are iterated until the crosspoint is found to converge(max. 5 iterations).

There are no standardized exact methods for determining thelimits for Schroeder integration and decay ﬁtting or noisecompensation techniques. Such methods are discussed below.

Response analysis may be further enhanced by estimating theamount of energy under the decay curve after the truncation
point. The measured decay curve is artiﬁcially extended beyond the point of truncation by extrapolating the regression
line on the late decay curve to inﬁnity. The total compensation energy is formed as an ideal exponential decay process,the parameters of which are calculated from the late decay
slope.

Limited integration or decay matching interval

Subtraction of noise ﬂoor level

There are several recommendations about dealing with noise
ﬂoor and the point where the decay meets noise. For example,
according to the ISO 3382 standard [1] for room reverberation
determination, the noise ﬂoor must be 10 dB below the lowest
decay level used for calculation of the decay slope. Morgan[5] recommends to truncate at the knee point and then tomeasure the decay slope of the backward integrated response
down to a level about 5 dB above the noise ﬂoor.

Chu [15] proposed a subtraction method in which the meansquare value of the background noise is subtracted from theoriginal squared impulse response before the backward integration. Curve (g) in Fig. 2 illustrates this case. If the noise
ﬂoor estimate is good and the noise is stationary, the resulting
backward integrated curve is close to the ideal decay curve.

Faiget et al. [16] propose a simple but systematic method forpost-processing noisy impulse responses. The latter part of a
response is used for the estimation of the background noiselevel by means of a regression line. Another regression line isused for the decay, and the end of the useful response is determined at the crossing point of the decay and backgroundnoise regression lines. The decay parameter ﬁtting intervalends at 5dB above the noise ﬂoor.

Hirata’s method
Hirata [7] has proposed a simple method for improving thesignal-to-noise ratio by replacing the squared single impulse
response h2 (t) with the product of two impulse responses
measured separately at the same position:Z∞h2 (t)dt ⇐t

1. The squared impulse response is averaged into local
time intervals in the range of 10–50 ms, to yield asmooth curve without losing short decays.

Z t∞[h1 (t) + n1 (t)][h2 (t) + n2 (t)]dt[h1 (t)h2(t) + h1 (t)n2 (t) +=tLundeby’s method
Lundeby et al. [8] presented an algorithm for automaticallydetermining the background noise level, the decay-noise truncation point, and the late decay slope of an impulse response.
The steps of the algorithm are:Z∞

Z ∞ n2+ h (t)n1 (t) + n1 (t)n2 (t)]dt=

h1 (t)h2 (t) +t+ n1 (t)n2 (t)[1 +

≈

Z∞

(9)

o

h1 (t)h2 (t)+] dt
n1 (t)n2 (t)h2 (t)dt + K(t)t

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

5

KARJALAINEN ET AL.

ESTIMATION OF MODAL DECAY PARAMETERS

The measured impulse responses consist of the decay terms
h1 (t), h2 (t) and the noise terms n1 (t), n2 (t). The highly
correlated decay terms h1 (t) and h2 (t) yield positive values
corresponding to squared response h2 (t), whereas the mutually uncorrelated noise terms n1 (t) and n2 (t) are seen as arandom ﬂuctuation K(t) superposed on the ﬁrst term.
Curves (h) in Fig. 2 illustrate a few decay curves obtainedby backward integration with Hirata’s method. In this simulated case they correspond to the case of (g), the noise ﬂoor
subtraction technique.

Other methods
Under adverse noise conditions, a direct determination of theT30 decay curve from the squared and time-averaged impulse
response has been noted to be more robust than the backward
integration method (Satoh et al. [17]).
3 NONLINEAR OPTIMIZATION OF A DECAYPLUS-NOISE MODEL
The nonlinear regression (optimization) method proposed by
Xiang [9] was shortly described above. In the present study
we have worked along similar ideas, using nonlinear optimization for improved robustness and accuracy. Below we introduce the nonlinear decay-plus-noise model and its applicationin several cases.
Let us assume that in noiseless conditions the system understudy results in simple exponential decay of response envelope, corrupted by additive stationary background noise. We
will study two cases that ﬁt into the same model category.
In the ﬁrst case there is a single mode (i.e., a complex conjugate pole pair in transfer function) that in the time domain
corresponds to an exponential decay function
hm (t) = Am e−τm t sin(ωm t + φm )(10)

Here Am is the initial envelope amplitude of the decaying sinusoidal, τm is a coeﬃcient that deﬁnes the decay rate, ωmis the angular frequency of the mode, and φm is the initial
phase of modal oscillation.
The second case that leads to a similar formulation is where
we have a high density of modes (diﬀuse sound ﬁeld) with exponential decay, resulting in an exponentially decaying noise
hd (t) = Ad e−τd t n(t)(11)

where Ad is the initial rms level of response, τd is a decayrate parameter, and n(t) is stationary Gaussian noise withrms level of 1 (= 0 dB).
In both Eqs. (10) and (11) we assume that a practical measurement of the system impulse response is corrupted withadditive stationary noisenb (t) = An n(t)(12)

where An is the rms level of Gaussian measurement noise on
the analysis bandwidth of interest, and it is assumed to be
uncorrelated with the decaying system response. Statisticallythe rms envelope of the measured response is thena(t) =qh2 (t) + n2b (t) =q

A2 e−2τ t + A2n(13)This is a simple decay model that can be used for parametricanalysis of noise-corrupted measurements. If the amplitude
4

envelope of a speciﬁc measurement is y(t), then optimized
least squares (LS) error estimate for parameters {A, τ, An }can be achieved by minimizing the following expression overa time span [t0 , t1 ] of interest

Zt

min

A,τ,An1[a(t) − y(t)]2 dt

(14)

t0

Since the model of Eq. (13) is nonlinear in parameters{A, τ, An }, nonlinear LS optimization is needed to search forthe minimum LS error.
By numerical experimentation with real measurement data it
is easy to observe that LS ﬁtting of the model of Eq. (14)places emphasis on large magnitude values, whereby noise
ﬂoors well below the signal starting level are estimated poorly.
In order to improve the optimization, a generalized form ofmodel ﬁtting can be formulated by minimization

Zt n
1

min

A,τ,Anf (a(t), t) − f (y(t), t)

o2

dt

(15)

t0

where f (y, t) is a mapping to balance weight for diﬀerent envelope level values and time moments.
The choice of f (y, t) = 20 log(y(t)) results in ﬁtting on the
decibel scale. It turns out that low-level noise easily has
a dominating role in this formulation. A better result inmodel ﬁtting can be achieved by using a power law scaling
f (y, t) = ys (t) with factor s < 1, which is a compromise between amplitude vs. logarithmic scaling. A value of s ≈ 0.5
has been found a useful default value4.
A time-dependent part of mapping f (y, t), if needed, can beseparated as a temporal weighting function w(t). A generalized form of the entire optimization is now to ﬁnd

Zt n
1

min

A,τ,Anw(t)as (t) − w(t)ys (t)

o2

dt

(16)

t0

There is no clear physical motivation for the magnitude compression factor s. Speciﬁc temporal weighting functions w(t)can be applied case by case, based on extra knowledge on the
behavior of the system under study and goals of the analysis,
such as focusing on the early decay time (early reverberation)of a room response.
The strengths of the nonlinear optimization method are apparent especially under extreme SNR conditions where all
three parameters, {A, τ, An }, are needed with the best accuracy. This occurs both at very low SNR conditions where
the signal is practically buried in background noise and at theother extreme where the noise ﬂoor is not reached within themeasured impulse response but still an estimate of the noiselevel is desired. The necessary assumption for the method towork in such cases is that the decay model is valid, implying
an exponential decay and stationary noise ﬂoor. Experiments
show the method to work both for single mode decay andreverberant acoustic ﬁeld decay models.
Figure 3 depicts three illustrative examples of decay modelﬁtting of a single mode at an initial level of 0 dB and diﬀerent noise ﬂoor levels. Because of simulated noisy responses it
is easy to evaluate the estimation accuracy of each parameter.
White curves show the estimated behavior of the decay-plusnoise model. In Fig. 3a the SN-ratio is only 6 dB. Errors inparameters in this case are a 0.5 dB underestimation of A,

Interestingly enough, this is close to the loudness scaling in auditory perception known from psychoacoustics [18].

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

6

KARJALAINEN ET AL.

ESTIMATION OF MODAL DECAY PARAMETERS

level/dB

3.5 % underestimate in decay time related to parameter τ ,and 1.8 dB overestimate of noise ﬂoor An . In Fig. 3b a similar case is shown with a moderate 30 dB SN-ratio. Estimation
errors of parameters are +0.2 dB for A, -2.8 % for decay time,and +1.2 dB for An . In the third case of Fig. 3c the SNratio is -60 dB so that the noise ﬂoor is hardly reached within
the analysis window. In this case the estimation errors are+0.002 dB for A, -0.07 % for decay time, and -1.0 dB for
An . This shows that the noise ﬂoor is estimated with highaccuracy also in this extreme case.

0
-10
-20

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

level/dB

0
-10
-20
-30
-40

level/dB

0
-20
-40
-60

Fig. 3: Nonlinear optimization of decay-plus-noise model forthree synthetic noisy responses with initial level of 0 dB andnoise levels of (a) -6 dB, (b) -30 dB), and (c) -60 dB. Black
curves show the Hilbert envelopes of the simulated responses
and white curves depict the estimated decay behavior.
The nonlinear optimization used in this study is based on using the Matlab function curvefit5 and the functions which
implement the weighting by parameter s and weighting function w(t) can be found at:
http://www.acoustics.hut.fi/software/decayof this level is important. A rule for an optimal value of thescaling parameter s is to use s ≈ 1.0 for very low SN-ratios,such as in Fig. 3a, and let it approach a value of 0.4–0.5 when
the noise ﬂoor is low as in 3c (see also Fig. 4 below).
4 COMPARISON OF DECAY PARAMETER ESTIMATION METHODS
The accuracy and robustness of methods for decay parameter estimation can be evaluated by using synthetic decay
signals or envelope curves, computed for sets of parameters{A, τ, An }. By repeating the same for diﬀerent methods, theirrelative performances can be compared. In this section we
present results from the comparison of the proposed nonlinear optimization and the method of Lundeby et al.
The accuracy of the two methods was analyzed in the following setting. A decaying sinusoid of 1 kHz with a 60 dB decay
time (‘reverberation time’) of 1 second was contaminated withwhite noise of Gaussian distribution and zero mean. Initialsinusoidal level to background noise ratio was varied from 0
dB to 80 dB by steps of 10 dB. Each method under study wasapplied to analyze the decay parameters and the error to the‘true’ value was computed in decibels for the initial and thenoise ﬂoor levels and as a percentage for decay time.
Figure 4 depicts the results of evaluation for the nonlinear
optimization proposed in this paper. The accuracy of decay
time estimation in Fig. 4a is excellent for SN-ratios above 30
dB and useful (below 10 % typically) even for SN-ratios of 0–
10 dB. Initial level is accurate within 0.1 dB for SNR above
20 dB and about 1 dB for SNR of 0 dB. Noise ﬂoor estimate
is within approximately 1–2 dB up to SNR of 60 dB and gives
better than a guess up to 70–80 dB of SNR (notice that SNR
alone is not important here but rather if the noise ﬂoor isreached in the analysis window or not).
%Our experience in the nonlinear decay parameter ﬁtting described here is that it still needs some extra information or
top-level iteration for best results. It is advantageous to select the analysis frame so that the noise ﬂoor is reached not
too early or late. If the noise ﬂoor is reached in the very beginning of the frame, the decay may be missed. Not reaching
the noise ﬂoor in the frame is a problem only if the estimate
5
6

average: —0.70 % (s=0.5)+0.04 % (s=1.0)

-5
-10

0

10

20

30

40

50

60

SNR/dB

dB
0

80

(b)

-0.5

average: +0.03 dB (s=0.5)—0.01 dB (s=1.0)

-1
-1.5

The optimization routines are found robustly converging inmost cases, including such extreme cases as Fig. 3a and 3c,and the initial values of parameters for iteration are not critical. However, it is possible that in rare cases the optimization
diverges and no (even a local) optimum is found6. It would
be worth working out a dedicated optimization routine guaranteeing a result in minimal computation time.

(a)

0

0

10

20

30

40

50

60

dB

SNR/dB80(c)0average:+1.2 dB+1.8 dB

-20
-40
0

10

20

30

40

50

60

SNR/dB80

Fig. 4: Sine-plus-noise decay parameter estimation errors(average of 20 trials) for the nonlinear optimization method
proposed in this paper as a function of SN-ratio: (a) decay
time estimation error in percentage, (b) initial level estimation error in dB, (c) noise ﬂoor estimation error in dB. Solidline: s = 0.5, dotted line: s = 1.0.

In new versions of Matlab, function curvefit is recommended to be replaced by function lsqcurvefit.
Function curvefit also prints warnings of computational precision problems even when optimization results are excellent.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

7

KARJALAINEN ET AL.

ESTIMATION OF MODAL DECAY PARAMETERS

Figure 5 plots the same information for decay parameter estimation using the method of Lundeby et al. without noisecompensation, implemented by us in Matlab. Since this iterative technique is not developed for extreme SN-ratios, suchas 0 dB, it cannot deal with these cases without extra tricks,
and even then it may have severe problems. We have usedsafety settings whereby it did not try to yield decay time values for SNR below 20 dB, and low SNR parts of the decay
parameter estimate curves are omitted.
%
6-5 dB ... noise floor+5dB integration

(a)

4
2

Figure 6 depicts three diﬀerent cases of ﬁtting the decay-plusnoise to this case of control room with a short reverberation
time. in Fig. 6a the ﬁtting is applied to the entire decay curve
including the initial delay, and the resulting model is clearlybiased towards too long reverberation time. In Fig. 6b the
initiall delay is excluded from model ﬁtting and the result isbetter. However, after the direct sound there is a period of
only little energy due to the ﬁrst reﬂections before the rangeof dnese reﬂections and diﬀuse response. If the reverberation time estimate should describe the decay of this diﬀuse
part, the case of Fig. 6c with ﬁtting, starting from about 30
ms, yields the best match to reverberation decay, and theapproaching noise ﬂoor is also estimated well7 .

-5 dB ... noise floor+10 dB integration

0
10

20

30

40

(b)

50

60

SNR/dB80

level/dB0
dB-2

average: —3.1 dB

-20
(a)
-40

-3

-60

-4

0

-5
20

30

40

50

60

SNR/dB

0.05

0.1

time/s80(c)average: —0.4 dB

-20
(b)
-40
-60

10

20

30

40

50

60

SNR/dB

080

Fig. 5: Decay parameter estimation errors for the Lundeby
et al. method as a function of SN-ratio: (a) decay time estimation error percentage with truncated Schroeder integration
but without noise compensation, (b) initial level estimation
error in dB, and (c) noise ﬂoor estimation error in dB.

level/dB

dB
10
5
0
-5
-10
0

10

level/dB

0

0.05

0.1

time/s

-20
(c)
-40
-60
0

0.05

0.1

time/sFor moderate SN-ratios the results of the method are fairly
good and robust. Decay time shows a positive bias of few percent except. The noise ﬂoor estimate is reliable in this caseonly up to about 50 dB SNR. Notice that the method is designed for practical reverberation time measurements rather
than this test case where it could be tuned to perform better.

Fig. 6: Decay-plus-noise model ﬁtting by nonlinear optimiza-

5 EXAMPLES OF DECAY PARAMETER ESTIMATION BY NONLINEAR OPTIMIZATIONThe next case deals with the modeling of low-frequency modesof a room. Below a critical frequency (so-called Schroeder
frequency) the mode density is low and individual modes canbe decomposed from the measured room impulse response.
The task here was to ﬁnd most prominent modes and toanalyze their modal parameters fm and τm , frequency anddecay parameter, respectively. The case studied was a hardwalled, partially damped room with moderate reverberation
time (≈ 1 sec) at mid and high frequencies, but much longer
decay times at lowest modal frequencies. The following procedure was applied:• A short-time Fourier analysis of the measured impulse
response was computed to yield a time-frequency representation, shown in Fig. 7 as a waterfall plot.
• At each frequency bin (1.3 Hz spacing is used), the dBscaled energy-time decay trajectory was ﬁtted to the
decay-plus-noise model with the nonlinear optimization
technique to obtain the optimal decay parameter τ .

In this section we present examples of applying the nonlinear estimation of a decay-plus-noise model to typical acoustic
and audio applications including reverberation time estimation, analysis and modeling of low-frequency modes of a room
response, and decay rate analysis of plucked string vibration
for model-based synthesis applications.
5.1

Reverberation time estimation

Estimation of the reverberation time of a room or a hall isrelatively easy if the decay curve behaves regularly and noise
ﬂoor is low enough. Often in practice the case is quite different. Here we demonstrate the behavior of the nonlinear
optimization method inan example where the measured impulse response includes an initial delay, irregular initial part,and a relatively high measurement noise ﬂoor.

tion to a room impulse response: (a) ﬁtting range includes
initial delay, transient phase, and decay, (b) ﬁtting includes
transient phase and decay, and (c) ﬁtting includes only the
decay phase.
5.2

Modeling of low-frequency room modes

7
In this example, decay parameter analysis is applied to the entire impulse response for demonstration purposes. In reality
it should be done as a function of frequency, i.e., to octave or 1/3-octave band decay curves.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

8

KARJALAINEN ET AL.

ESTIMATION OF MODAL DECAY PARAMETERS• Based on decay parameter values and spectral levels, arule was written to pick up the most prominent modal
frequencies and the related decay parameter values.

of tones that need to be analyzed are their fundamental frequency and the decay time of all harmonic partials that are
audible. While the estimation of fundamental frequency isquite easy, the measurement of decay times of harmonics (=modes of the string) is complicated by the fact that they all
have a diﬀerent rate of decay and also the initial level canvary within a range of 20-30 dB. There may also not be information about the noise ﬂoor level for all harmonics.
One method used for measuring the decay times is based on
the short-time Fourier analysis. A recorded single guitar tone
is sliced into frames with a window function in the time domain. Each window function is then Fourier transformed
with the FFT using zero-padding to increase the spectral
resolution, and harmonic peaks are hunted from the magnitude spectrum using a peak-picking algorithm. The peakvalues from the consecutive frames are organized as tracks,
which correspond to the temporal envelopes of the harmonics. Then it becomes possible to estimate the decay rate ofeach harmonic mode. In the following, we show how thisworks with the proposed algorithm. Finally, the decay rateof each harmonic is converted into a corresponding target response, which is used for designing the magnitude responseof a digital ﬁlter that controls the decay of harmonics in thesynthesis model.

Fig. 7: Time-frequency plot of room response at low frequencies. Lowest modes (especially 40 Hz) show long decay times.

-20
-40
-60

level/dB

(a)

0
0

0.5

1

1.5

2

2.5

3

3.5

4

time/s5-40(b)level/dBlevel/dB

In this context we are interested in how well the decay parameter estimation works with noisy measurements. Application of the nonlinear optimization resulted in decay parameter estimates, some of which are illustrated in Fig. 8 by acomparison of the original decay and decay-plus-noise model
behavior. For all frequencies in the vicinity of a mode themodel ﬁts robustly and accurately.

Figure 9 plots three examples of modal decay analysis of guitar string harmonics (string 5, open string plucking). Harmonic envelope trajectories were analyzed as described above.
The decay-plus-noise model was ﬁtted in a time window thatstarted from the maximum value position of the envelope
curve. In case (a), the second harmonic shows a highly regular decay after initial transient of plucking whereby decay
ﬁtting is almost perfect. Case (b), harmonic number 24, depicts a strongly beating decay where probably the horizontal
and vertical polarizations have a frequency diﬀerence that after summation results in beating. Case (c), harmonic 54,
shows a harmonic trajectory where the noise ﬂoor is reachedwithin the analysis window. In all shown cases, the nonlinear
optimization works as perfectly as a simple decay model can
do.

-20
0
0
1

1.5

2

2.5

3

3.5

4

-40

time/s

5

(c)

-60

0

0.5

1

1.5

2

2.5

3

3.5

4

time/s5

Fig. 8: Fitting of decay-plus-noise model to low-frequency
modal data of a room (see Fig. 7): (a) at 40 Hz mode, (b) at
104 Hz, and (c) at 117 Hz (oﬀ-mode fast decay). Dotted line:
measured, dashed line: modeled.
5.3

Analysis of decay rate of plucked string tones

Model-based synthesis of string tones can produce realistic
guitar-like tones if the parameter values of the synthesis model
are calibrated based on recordings [2]. The main properties

level/dBlevel/dB

0.5level/dB

-80

-80

-10
-15

-60

0

(a)

-5

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

time/s
(b)

-5
-10
-15
-20
0
0
-10
-20
-30
-40
-50
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

time/s
(c)

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

time/sFig. 9: Examples of modal decay matching for harmonic components of a guitar string: (a) regular decay after initial transient, (b) strongly beating decay (double mode), and (c) fast
decay that reaches the noise ﬂoor. Dotted line: measured
envelope, solid line: optimized ﬁt.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

9KARJALAINEN ET AL.
As can be concluded from case (b) of Fig. 9, a string can
exhibit more complicated behavior than a simple exponential decay. Even more complex is the case of piano tones
because there are 2–3 strings slightly oﬀ-tuned, and the envelope ﬂuctuation can be more irregular. Two-stage decay isalso common where the initial decay is faster than later decay.
In all such cases a more complex decay model is needed toachieve a good match with the measured data. It remains afuture problem to investigate ﬁtting such models using linear
optimization techniques.
6

SUMMARY AND CONCLUSIONSIn this paper, an overview of modal decay analysis methods for noisy impulse response measurements of reverberant
acoustic systems is presented, and further improvements are
introduced. The problem of decay time determination is important for example in room acoustics for characterizing thereverberation time. Another application where a similar problem is encountered is when estimating string model parameters for model-based synthesis of plucked string instruments.
It is shown in this article that the further developments of the
decay-plus-noise model yield highly accurate decay parameter estimates, outperforming traditional methods especially
under extreme SN-ratio conditions. Challenges for further
studies remain to make the method (with increased numberof parameters) able to analyze complex decay characteristics,
such as double decay behavior and strongly ﬂuctuating responses due to two or more modes very close in frequency.
The Matlab code for nonlinear optimization of decay parametres, including data examples, can be found at:
http://www.acoustics.hut.fi/software/decay
ACKNOWLEDGMENTS
This study is part of the VÄRE technology program, project
TAKU (Control of Closed Space Acoustics), funded by Tekes(National Technology Agency). The work of Vesa Välimäki
has been ﬁnanced by the Academy of Finland.

Notes on reverberation time deﬁnitions
Reverberation time T60 or just T is deﬁned as the time in seconds it takes for the energy in the steady-state sound ﬁeld ina room to decay 60 dB after the source of sound excitation issuddenly turned oﬀ [19]. According to deﬁnition by Beranek[20], the reverberation time might be deﬁned in any of thefollowing ways:
1. the interval between turning oﬀ the sound source andthe time when the instantaneous value of the pressurelevel ﬁrst falls to 60 dB below its steady state value
2. the time interval until the average value of the ﬂuctuating pressure ﬁrst falls to this value
3. the time until the average value of the ﬂuctuating log
p decay curve ﬁrst falls to 60 dB below its steady statevalue, or
4. the time required for the evaluated decay slope to dropoﬀ 60 dB.
In practice, the reverberation time is often determined fromthe slope of a decay curve using only the ﬁrst 25 or 35 dBof decay and extrapolating the result to 60 dB. For the recommended practice of reverberation time determination, seestandard [1].

ESTIMATION OF MODAL DECAY PARAMETERSREFERENCES[1] ISO 3382-1997, Acoustics – Measurement of the Reverberation Time of Rooms with Reference to Other Acoustical Parameters. Geneva, Switzerland: International
Standards Organization, 1997. 21 p.
[2] V. Välimäki, J. Huopaniemi, M. Karjalainen, andZ. Jánosy, “Physical Modeling of Plucked String Instruments with Application to Real-Time Sound Synthesis,”
J. Audio Eng. Soc., vol. 44, no. 5, pp. 331–353, 1996
May.
[3] M. R. Schroeder, “New Method of Measuring Reverberation Time,” J. Acoust. Soc. Am., vol. 37, pp. 409–412,
1965.
[4] M. R. Schroeder, “Integrated-impulse Method Measuring Sound Decay Without Using Impulses,” J. Acoust.
Soc. Am., vol. 66, no. 2, pp. 497–500, 1979 Aug.
[5] D. Morgan, “A parametric error analysis of the backward
integration method for reverberation time estimation,”
J. Acoust. Soc. Am., vol. 101, no. 5, pp. 2686–2693, 1997
May.
[6] L. Faiget, C. Legros, and R. Ruiz, “Optimization of
the Impulse Response Length: Application to Noisy andHighly Reverberant Rooms,” J. Audio Eng. Soc., vol. 46,
no. 9, pp. 741–749, 1998 Sept.
[7] Y. Hirata, “A Method of Eliminating Noise in Power
Responses,” J. Sound Vib., vol. 82, no. 4, pp. 593–595,
1982.
[8] A. Lundeby, T. E. Vigran, H. Bietz, and M. Vorländer,“Uncertainties of Measurements in Room Acoustics,”
Acustica, vol. 81, pp. 344–355, 1995.
[9] N. Xiang, “Evaluation of Reverberation Times Using a
Nonlinear Regression Approach,” J. Acoust. Soc. Am.,
vol. 98, no. 4, pp. 2112–2121, 1995 Oct.
[10] A. V. Oppenheim and R. W. Schafer, Digital Signal Processing. Englewood Cliﬀs, N.J.: Prentice-Hall, 1975.
(Chapter 10, pp. 480–531).
[11] MathWorks, Inc., MATLAB User’s Guide, 2001.
[12] S. M. Kay, Fundamentals of Statistical Signal Processing: Volume I: Estimation Theory. Englewood Cliﬀs,
N.J.: Prentice-Hall, 1993.
[13] J. D. Markel and J. A. H. Gray, Linear Prediction of
Speech. Berlin, Germany: Springer-Verlag, 1976.
[14] MathWorks, Inc., MATLAB Signal Processing Toolbox,
User’s Guide, 2001.
[15] W. T. Chu, “Comparison of reverberation measurements
using Schroeder’s impulse method and decay-curve averaging method,” J. Acoust. Soc. Am., vol. 63, no. 5,
pp. 1444–1450, 1978 May.
[16] L. Faiget, R. Ruiz, and C. Legros, “Estimation of Impulse Response Length to Compute Room Acoustical
Criteria,” Acustica, vol. 82, no. Suppl. 1., p. S148, 1996
Sept.
[17] F. Satoh, Y. Hidaka, and H. Tachibana, “Reverberation
Time Directly Obtained from Squared Impulse Response
Envelope,” in Proc. Int. Congr. Acoust., vol. 4, (Seattle,
WA), pp. 2755–2756, 1998 June.
[18] E. Zwicker and H. Fastl, Psychoacoustics — Facts andModels. Berlin: Springer-Verlag, 1990.
[19] W. C. Sabine, Architectural Acoustics. 1900. Reprint by
Dover, 1964, New York.
[20] L. Beranek, Acoustic Measurements. 1949. Reprint by
Acoust. Soc. Am., 1988.

AES 110TH CONVENTION, AMSTERDAM, THE NETHERLANDS, 2001 MAY 12–15

10