Audio Engineering Society

Convention Paper 5590
Presented at the 112th Convention
2002 May 10–13
Munich, Germany

This convention paper has been reproduced from the author’s advance manuscript, without editing, corrections, or
consideration by the Review Board. The AES takes no responsibility for the contents. Additional papers may beobtained by sending request and remittance to Audio Engineering Society, 60 East 42nd Street, New York, New York
10165-2520, USA; also see www.aes.org. All rights reserved. Reproduction of this paper, or any portion thereof, is
not permitted without direct permission from the Journal of the Audio Engineering Society.

AR/ARMA Analysis and Modeling of Modes
in Resonant and Reverberant Systems
Matti Karjalainen1 , Paulo A. A. Esquef1 , Poju Antsalo1 , Aki Mäkivirta2 , and Vesa Välimäki3,1
1 Helsinki University of Technology, Lab. of Acoustics and Audio Signal Processing, P.O. Box 3000, FIN-02015 HUT, Finland
2 Genelec Oy, Olvitie 5, FIN-74100 Iisalmi, Finland
3 Tampere University of Technology, Pori School of Technology and Economics, P.O. Box 300, FIN-28101 Pori, Finland

Correspondence should be addressed to Matti Karjalainen (matti.karjalainen@hut.fi)

ABSTRACT
Discrete-time analysis and modeling of reverberant and resonating systems has many applications in audio andacoustics. In a recent paper (AES110, Preprint 5290) we formulated techniques for the estimation of modal decay
parameters from noisy response measurements, targeting to systems such as room reverberation and modal decay
as well as musical instrument modeling. In this paper we extend the methodology to AR and ARMA modeling of
measured responses by all-pole and pole-zero ﬁlters. In addition to an overview of standard techniques we propose aspectral zooming technique that is useful for resolving very closely positioned modes and high-density modal clusters.
Sensitivity to background noise is also studied. Application cases are taken from analysis and modeling of room
responses, loudspeaker-room equalization, and estimation of parameters for musical instrument modeling.

1

INTRODUCTION

Parametric analysis and modeling is an increasingly common
task in acoustics and audio. In this paper we focus on suchaudio related problems where a target system response canbe measured and the task is to model it for computational

simulation/synthesis or to inverse model it for equalization,i.e., compensation for nonidealities in the observed response.
Representative examples of the ﬁrst group, forward simulation
for synthesis, are room response modeling including artiﬁcial
reverb design, or just estimation of eigenmodes at low fre-

KARJALAINEN ET AL.

AR/ARMA MODELING

quencies in room acoustical studies, and modeling of musicalinstruments. The second category of inverse modeling is common in audio, where equalization of non-ideal response properties is a frequent task in high-quality sound reproduction, as
well as prevention of acoustic feedback in sound reinforcement
systems.
The behavior of acoustic/audio systems at low frequencies canoften be modeled analytically and parametrized, at least withlumped models, for example low-frequency modal behavior ofa rectangular room [1], musical instrument body [2], or loudspeaker enclosure [3]. For irregular structures and/or higher
frequencies it is much more diﬃcult or impossible to ﬁnd analytical or numerical models that are useful in practice. In
such cases it is still possible to measure system responses andto apply signal modeling techniques to analyze, simulate, or
real-time synthesize a given response.
Particular interest of this study is focused on the resonant
and reverberating characteristics of complex acoustic systems.
The modal behavior, i.e., decomposition of eigenmodes assuming a linear and time-invariant (LTI) system, can be extremely complicated. In one-dimensional resonators, such as
strings and tubes in musical instruments, the density of eigenmodes is not necessarily high, but these modes can exhibit
complicated details, for example two-stage decay and strong
beating in a decaying envelope. In 2-D and 3-D resonators,
such as membranes, plates, and enclosed spaces, the modal
density increases towards higher frequences, resulting in dense
modal patterns and reverberation type behavior when neighboring eigenmodes overlap essentially. Also the temporal behavior becomes complex, especially in rooms where the direct
sound and early reﬂections are followed by an increasing density of reﬂections resulting in late reverberation. Thus the
problem of modeling measured responses may beneﬁt of atime-frequency modeling viewpoint where the properties ofauditory perception are taken into account.
A rich literature exists on signal modeling of LTI systems [4,5, 6, 7] from many branches of engineering, systems sciences,
and applications. Software tools for modeling are available,
for example in MATLAB1 . Here we assume that the target
system and the desired model are, in addition to linearity
and time invariance, also stable and causal. We also assume
that the measured response is not heavily noise-contaminated
so that estimating the system model parameters is practical
for the applications at hand. Then the measured impulse
response h(n) can be approximated by a rational expression
in the z-transform domain by
H(z) =

N−kb0 + b1 z−1 + . . . + bN z−N
k=0 bk z=
P−1−P
1 + a1 z+ . . . + aP z
1 + k=1 ak z−k(1)which makes it possible to eﬃciently simulate or synthesize
the system by various digital ﬁlter implementations [5, 10, 11]of the given transfer function.
In this paper we only deal with discrete-time representations
for digital signal processing. Thus the easiest way “to model”
a measured response h(n) or its truncated/windowed version−k .
is to take it directly as an FIR ﬁlter H(z) = N
k=0 h(k)z
For complex systems, the length of a ﬁnite impulse response1
In this study we utilize particularly the functions found in
the Signal Processing Toolbox [8]. Another MATLAB toolbox ofinterest is the System Identiﬁcation Toolbox [9].

required for suitable representation may be too long, preventing real-time implementations. On the other hand, shortening
the ﬁlter length reduces the capabilities of identifying the inherent resonant properties of the system under study. IIR
ﬁlters for inﬁnite impulse response can come in two forms:(a) all-pole models where the numerator of Eq. (1) is reduced
to a single gain coeﬃcient b0 or (b) pole-zero models withboth the numerator and denominator being non-trivial polynomials of z.
In systems science and engineering, for example in control
theory for estimation and identiﬁcation tasks, the terms autoregressive (AR), moving average (MA), and autoregressive
moving-average (ARMA) have been used for modeling processes similar to all-pole, FIR, and pole-zero ﬁlter modeling,respectively. For the sake of convenient use of abbreviations
AR and ARMA, as well as to draw attention to the rich knowledge from various other ﬁelds than digital audio signal processing, we apply the terms MA, AR, and ARMA here whenreferring to speciﬁc types of models.
Our aim does not stop in obtaining a useful approximation ofa measured target system by a transfer function of Eq. (1).
We are interested in decomposing it into a parametric description of its constituent components, particularly the complexconjugate pole pairs, i.e., the complex-valued roots of thedenominator polynomial, which represent the eigenmodes ofthe system and result in the resonant and reverberant behavior. In theory they are common in spatially distributed
systems such as a room [12] or an instrument body, whileN−k ) are essentiallyzeros (roots of the numerator
k=0 bk z
position-dependent. We are interested in accurate estimationof modal parameters, such as pole angle and radius or, equivalently, mode frequency and decay coeﬃcient.

In a recent paper [13] we studied this problem of modal parameter estimation using traditional time-frequency analysis
techniques by ﬁrst trying to isolate potential mode frequencies and then estimating the modal decay rate from a spectrogram, such as short-time Fourier analysis or cumulative
decay spectrum [14]. Decay rate estimation was also applied
to wide-band signals, for example to robust estimation of reverberation time T60 . The problem with such methods is tomodel overlapping modes that result in non-exponential decay
in any reasonable frequency span. AR and ARMA approaches
try to model the target response as a whole by minimizing agiven modeling error criterion, typically a least squares (LS)error. Thus the interactions of overlapping modes are taken
into account simultaneously and systematically.
One problem of straightforward optimization is the inﬂexibility of the global optimization criteria, for example to take
into account the properties of diﬀerent modes. Also the selection of proper values for the order parameters N and P of
Eq. (1) is not easy. A practical problem is that solving the(complex-valued) roots of a high-order polynomial is a diﬃcult numerical task. To avoid problems with high-order models we use methods where a part of the given audio frequency
range is modeled at a time to obtain an accurate descriptionof the modes within this frequency span. Frequency zooming
ARMA modeling (FZ-ARMA) is shown to be a powerful way
to decompose highly complex resonant responses into modal
representations, and related IIR ﬁlter implementations canbe used for simulation and synthesis of such systems. Our research is particularly inﬂuenced by earlier studies on selective
linear prediction [15], multiband modeling of musical signals

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

2

KARJALAINEN ET AL.

AR/ARMA MODELING

2

AR MODELINGThe impulse response of a resonant system shows one or more
exponentially decaying sinusoids. Each such ‘mode’ can be
inherently modeled by a complex conjugate pole pair, whichsuggests AR modeling with corresponding inﬁnite impulse response (IIR) ﬁlters. There is a long tradition for ﬁnding leastsquares optimal ﬁt of such models to measured LTI system
responses, either to a given impulse response or to inputoutput signal pairs. Here we shortly refer to the theory of
linear prediction (LP) which has found application especiallyas a powerful spectral modeling technique in speech processing [18, 19].
2.1

1e- t/τ sin (ω t)e- t/τ cos(ω t)

(a)

0

0

-0.5

-0.5

-1

0

0.2

0.4
0.6
normalized time

0.8

1

40

-1

0

0.2

0.4
0.6
normalized time

0.8

(d)
20

0
-20
0

1

0
-20
0.02
0.04
0.06
0.08
normalized frequency f/fs

0.1

-40

0

2

0.02
0.04
0.06
0.08
normalized frequency f/fs

(e)
1

0

0

-0.5

-1
0

0.2

0.4
0.6
normalized time

0.8

0.1
(f)

0.5

-1

1

40
(c)

20

-40

(b)

0.5

level [dB]

This paper is structured as follows. Sections 2 and 3 present
an overview of AR and ARMA modeling methods and techniques. Examples are given to illustrate the modeling ability
and limitations of these basic approaches. Section 4 introduces the Frequency-Zooming ARMA (FZ-ARMA) methodthat is able to analyze high-order systems with overlapping
modes and dense modal distributions. The eﬀects of nonidealities, including noise and non-LTI behavior, are discussed.
Three cases of audio applications that use AR/ARMA/FZARMA methods are described in Section 5, including modeling of measured room responses, inverse modeling and equalization of loudspeaker-room responses, and modeling as well
as sound synthesis of musical instruments. Finally a discussion and concluding remarks are given.

1
0.5

level [dB][16] or loudspeaker responses [17], and other high-resolution
system modeling techniques.

1

-2

0

0.2

0.4
0.6
normalized time

0.8

1

Fig. 1: AR modeling of single mode decay with diﬀerent initial phase: (a) damped sinusoid target response, (b) damped
cosine target function, (c) magnitude spectrum of sine response, (d) magnitude spectrum of cosine response, (e) AR
model response for sinusoidal response with order of 2, and(f) AR model response for cosine case with order of 2 (solidline) and order of 20 (dashed line).

Linear prediction

In linear prediction a signal sample x(n) is assumed to be predictable as a linear combination of previous samples x̂(n) =P
i=1 ai x(n − i). When the least squares prediction error
between x(n) and x̂(n) is minimized, the (auto)correlationcoeﬃcients x(i) x(i + k)

P −1

rx (k) =(2)

i=0

play a central role. The most frequently used version of
LP analysis is the autocorrelation method, where the optimalmodel parameters ai are solved from a linear matrix equation(normal equations)r
 r
 .r. .
0
1
2

r1
r0
r1
...

r2
r1
r0
...

rP −1

rP −2

rP −3

...
...
...
...
...

rP −1
rP −2
rP −3
...
r0

  a  r   a   r  .a. . = .r. . (3)
1
2
3

1
2
3

aP

rP

to do time domain ﬁtting to given target response unless thetarget process is strictly of AR (all-pole) type. Figure 1 illustrates this clearly in a simple case of one idealized mode. Fora damped sinusoid in Fig. 1(a) the model response in Fig. 1(e)is a very accurate replica of the target already with model order of 2 (one complex conjugate pole pair). If the initial phaseis changed 90 degrees to a damped cosine function, the second
order AR model response (solid line) in Fig. 1(f) severely falls
oﬀ from the given target in Fig. 1(b).
Further insight into the behavior of AR modeling is gained bynoticing that the model is based entirely on the autocorrelation coeﬃcients, which in the frequency domain corresponds
to the power spectrum. This means pure spectral modeling,whereby the spectrum of the cosine case in Fig. 1(d) clearly
deviates from the spectrum of sine case response in Fig. 1(c).
Irregular onsets are common in acoustic system responses,
thus indicating that simple AR modeling will have diﬃculties, and that more powerful methods are needed in accurate
temporal modeling.

Here, rk are the autocorrelation coeﬃcients rx (k) from Eq. (2)for a signal frame under study and P is the order of LP analysis (order of the all-pole model ﬁlter). The coeﬃcients ai
are the estimated polynomial coeﬃcients in the denominator−i , the numerator being only aof Eq. (1), i.e., 1 + P
i=1 ai z
gain coeﬃcient2.

A somewhat better match to the decaying cosine tail above
can be achieved by the covariance method of AR analysis3
but the phase matching problem remains and requires ARMA
modeling.

2.2

3



Limitations of AR modeling

A problem with AR modeling of real-world systems, in theformulation described above, is that the method is not able2
Linear predictive analysis is computed in MATLAB by function lpc, which ﬁrst solves the autocorrelation coeﬃcients rk andthen inverts the Toeplitz type correlation matrix to yield predictor coeﬃcients ai through levinson recursion function [8].

ARMA MODELING

ARMA modeling, yielding a pole-zero ﬁlter, has more modeling power than the AR method. It is however more diﬃcult
because no closed-form solutions are available, thus requiring
3
The covariance method function ar, option ls, can be foundin the System Identiﬁcation Toolbox of Matlab.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

3

KARJALAINEN ET AL.

Two methods are applied below: Prony’s method andSteiglitz-McBride method4. A brief discussion is given on
the problem of model order selection followed by motivationof the need for improved frequency resolution, before introducing the frequency-zooming ARMA technique.
3.1

Prony’s method

Prony’s method [20, pp. 226–228] is a stepwise algorithm
that ﬁts N + 1 ﬁrst samples of given response exactly, while
P poles of the denominator in Eq. (1) take care of tail decay
ﬁtting. Because the AR part estimation is of covariance type,the resulting ﬁlter can become unstable, even in cases where
the target system to be modeled is stable.
3.2

Steiglitz-McBride iteration

Steiglitz-McBride method [7, pp. 174–177] is an algorithmwith iterative preﬁltering for least squares ﬁt of an ARMA
model of Eq. (1) to a given impulse response or given inputoutput pair (system identiﬁcation problem). An initial estimate for the denominator can be obtained for example by
Prony’s method.
As with Prony’s method, the resulting ﬁlter from the SteiglitzMcBride iteration can be unstable especially with high ﬁlter
orders even for stable target systems. Often the model response starts in good match with the given time domain response (since this is LS ﬁtting) but after some time it starts
to explode due to pole(s) outside the unit circle.
3.3

Model order selection

Both AR and ARMA models need careful selection of ﬁlter orders N and P (N = 0 for AR models). There is no general andautomatic way to select optimal ﬁlters order(s), rather they
can be searched for by various rules to obtain good enoughmatch to a given target response [21] or the order(s) can beapproximated using a priori information about the target system to be modeled.
An illustrative example on how the order of AR modeling (linear prediction) aﬀects the estimated modal frequencies (poleangles) is shown in Fig. 2. The magnitude spectrum of a measured room response is plotted in comparison with a relatedpole frequency bifurcation map. For lowest orders only thespectral peak ranges become roughly approximated, and forincreasing model order these poles split into new pole pairs
and groups of poles.
In this study we are interested in modeling resonant and reverberant systems by methods where the poles and relatedparameters, angles and radii, can be resolved explicitly. Thisis needed in applications such as those discussed below in
the application section. The selection of model order is thenmore demanding than in cases where ﬁnding the numerator
and denominator coeﬃcients is enough. If the model order
is too low, not all modes are represented properly by complex conjugate pole pairs, or the radii of poles found remain
underestimated. If a model order is too high, single modes
4
These methods are available as Matlab Signal Processing
Toolbox functions prony and stmcb.

0level / dB

-10
-20
-30
-40
0
100
model ordernonlinear optimization. ARMA estimation algorithms are iterative, starting from a good AR model and then iterating
towards optimal parameter values. As with any nonlinear optimization methods, problems of divergence or convergence toa local optimum may occur, and plain computational problems due to the lack of numerical precision are found.

AR/ARMA MODELING

20

40

60

80

100

120

140

160

180

200

f/Hz

20

40

60

80

100

120

140

160

180

200

f/Hz

80
60
40
20
0

Fig. 2: Magnitude spectrum of a measured room response(top) and frequencies corresponding to pole angles obtained
from linear prediction of varying order (bottom).

often become overmodeled, i.e., more than one pole pair will
become allocated per mode.
3.4

Limitations of ARMA modeling

While powerful in simple cases for low model orders, the methods described above exhibit diﬃculties with high orders andcomplex target systems. Often these problems originate fromlimited computational precision. In Prony’s method and theSteiglitz-McBride iteration the potential instability is often aproblem. Although poles outside the unit circle can be mirrored inside the unit circle, yielding an equivalent magnitude
spectrum, the temporal structure of the impulse response ischanged.
Linear prediction (autocorrelation method) may yield stable
and accurate results with model orders of hundreds or thousands, particularly when the poles are not very close to theunit circle or to each others. However, if the poles (and zeros)must be solved explicitly, which is of interest in this study,
solving for roots of the denominator (and numerator) is typically more problematic than obtaining these polynomials per
se. Frequency-selective AR and ARMA modeling solves someof these problems.
4

FREQUENCY-ZOOMING ARMA (FZ-ARMA)

The problems in resolving very closely positioned modes andmode groups was the reason in this study to experiment withmethods that have better control over frequency resolution.
Several ideas are available for improvement, including frequency warping [22] and frequency selective modeling such as
selective linear prediction [15], multiband AR/ARMA techniques [16], and many other high-resolution signal modelingmethods.
Frequency warping is a convenient technique when either lowest or highest frequencies require enhanced frequency resolution. This approach can be extended to Kautz ﬁlters thatexhibit interesting properties of generalized frequency resolution control [23]. These methods are, however, out of thescope of this study.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

4

KARJALAINEN ET AL.

AR/ARMA MODELING

Frequency selective modeling has been applied for example inlinear prediction of speech. In a simple case a target response
can be low-pass ﬁltered and decimated in order to model thelow-frequency part of the response. A range of higher frequencies can be modulated down and decimated before similar
modeling. Actually any subband of a given frequency range
can be modeled this way, and ﬁnally the resulting parameters(poles and zeros) can be mapped back to the corresponding
original frequency range. This approach is called here modeling by frequency zooming. It resembles also the multiband
techniques used in [16, 17, 24].
4.1

Formulation of frequency zooming

The FZ-ARMA (or FZ-AR) analysis starts by modulating(heterodyning) the desired frequency range down to the neighborhood of zero frequency [25, 26, 27] byhm (n) = ejΩm n h(n)(4)

where Ωm = 2πfm /fs , fm is the modulation frequency, andfs is the sample rate. In the z-domain this means clockwise
rotation of poles (and zeros) zi by angle Ωm , i.e.,Ωi,rot = Ωi − Ωm = arg(zi ) − Ωm

(6)

Together mappings (5) and (6) yield new poles (and zeros)5ẑi = |zi | Kzoom ej(arg(zi )−Ωm )Kzoom(7)

Now it is possible to apply any AR or ARMA modeling to
the modulated and decimated response obtained from h(n).
Notice that this new signal is complex-valued due to the onesided modulation operation. This approach resembles multirate and subband techniques.
The advantage gained by frequency zooming is that in thezoomed subband the order of (ARMA) analysis can be reduced by increasing the zooming factor Kzoom and, consequently, the solution of poles and zeros as roots of denominator and numerator polynomials of the model function Eq. (1)is simpliﬁed. Additionally this means that a diﬀerent resolution can be used in each subband, for example based on
knowledge about the modal complexity in each subband.
After solving the poles and zeros within a zoomed subband,they must be remapped to the full band. This means inverse
scaling the radii of poles (and zeros) as well as rotating themcounter-clockwise, i.e.,(1/Kzoom ) jΩm

zi = ẑie

In the investigations of FZ-ARMA below the frequencyzooming method of solving the ARMA coeﬃcients is the
Steiglitz-McBride method. Notice also that the ﬁlter orders N
and P refer to real-valued ﬁlters with complex conjugate pairs
constructed from one-sided zeros and poles obtained from theabove method. Thus the orders are twice the numbers of zeros
and poles from the above procedures, respectively.
4.2

Modeling of higher-order modal decays

In this section, the performance of the FZ-ARMA analysis
method is illustrated through synthetic signals. In particular we are interested in investigating the modeling capability
when dealing with signals exhibiting beating and/or two-stage
decay in their envelopes. Simple signals featuring these characteristics can be obtained by(5)

but retaining the pole (or zero) radius. The next step toincrease frequency resolution is to limit the frequency range
by decimating, i.e., lowpass ﬁltering and down-sampling therotated response by a zooming factor Kzoom to obtain a new
sampling rate fs,zoom = fs /Kzoom . This implies mapping to
a new z-domain where poles (and zeros) are scaled by the rulezi,zoom = ziKzoom

task but on the other hand, partitioning of the whole problem and recombining again brings increased ﬂexibility. It isadvantageous to pick poles and zeros only within the central
part of each (overlapping) subband to avoid problems at theboundaries of subbands due to band-limitation.

(8)

Due to the one-sided down-modulation in (4), each pole (andzero) zi must be used as a complex conjugate pair in order toobtain real-valued ﬁlters.
The ﬁnal step is to combine poles and zeros obtained from
diﬀerent subbands into a full model. This is a non-trivial
5Eqs. (5) and (6) characterize how the z-domain properties of agiven response are changed through modulation and decimation,
but the estimated pole-zero pattern of an ARMA model will beobtained only in the next step.

a eMs(n) =

k

k=1

−n/τk

sin(2πnfk+ θk )
fs(9)

where M is the number of modal frequencies present in s(n),τk are the decay time constants, fk the modal frequencies,
fs the sampling frequency, and θk the initial phases of themodes.
Let us start with case (A) in which the amplitude envelope ofthe signal consisting of two modes shows beating. The parameters used to generate the signal as well as those adopted in
the FZ-ARMA modeling are given in Table 1. The target responses in sine and cosine phase, their FZ-ARMA envelopes,
and resynthesized versions are shown in Fig. 3. The envelopes
are obtained from the complex decimated signals.
Each resynthesized response is computed as the impulse response of the ﬁlter that is obtained by combining the complex
conjugate poles and zeros from the FZ-ARMA analysis. Notice that the band-limitation in frequency focusing producescompensating zeros (and poles) that show artifacts in the reconstructed impulse response. In a typical case of focusing on
a narrow resonance band the reconstructed impulse response
has a high-amplitude impulse as its ﬁrst sample (see also [27]).
This has been removed the in resynthesized response simulations of this paper. This removal has only marginal eﬀect on
the spectral properties of the response.
In the simulation results of case (A) in Fig. 3, an ARMA(4,4)
model suﬃces to represent properly the envelope decays insubplots (c) and (d), while the initial phase characteristics ofthe resynthesized signal in (f) deviate form (b). Note that it
is almost impossible to distinguish between the dashed andsolid lines in the subplots (c) and (d) of Fig. 3.
In case (B) we verify the FZ-ARMA modeling of a two-mode
response for which the amplitude envelope exhibits a twostage decay. The parameters used to generate this signal,
as well as those of FZ-ARMA modeling, are summarized inTable 2, and the results of modeling are shown in Fig. 4. The
slower decaying mode is modeled properly although its initiallevel is 10 dB below the stronger one. This ability of two-stage
decay analysis can work down to -30 dB in a clean synthetic
case.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

5

KARJALAINEN ET AL.

AR/ARMA MODELING
Table 2: Parameters in simulation case (B).

Table 1: Parameters in simulation case (A).

(c)

−20
−40
0

0.5
normalized time

1

1

0
−1
0

0.5
normalized time1

0.5
normalized time
(d)

−20
0.5
normalized time

1

1

(f)

0
−1
0

0.5
normalized time

1

1

(a)

0
−1
0

1

Synthetic Signal
a1
a2
τ1
τ2
0.2 0.8 0.3 0.02
FZ-ARMAfm /HzN
44.87
4amplitude

0

−40
0

f2 /Hz
100

Kzoom
200

(b)

0

amplitude

(e)

P
4

1

−1
0

1
envelope dB

0.5
normalized time0

amplitudeenvelope dB−10

f1 /Hz100

0.5
normalized time

0

(c)

−10
−20
0

0.5
normalized time

1

1

(e)

0
−1
0

0.5
normalized time

1

θ1
0

θ2
0
P
4

1

(b)

0
−1
0

1
envelope dB

0

θ2
0

0.5
normalized time

0

1
(d)

−10
−20
0

amplitude(a)amplitude

amplitude1

θ10

amplitude

Kzoom
200

Synthetic Signal
a1
a2
τ1
τ2
0.5 0.5 0.07 0.07
FZ-ARMAfm /HzN
52.37
4envelope dB

f2 /Hz115

amplitude

f1 /Hz100

0.5
normalized time

1

1

(f)

0
−1
0

0.5
normalized time1

Fig. 3: Case (A), amplitude beating due to two modes with

Fig. 4: Case (B), two-stage decay due to two modes withnearby frequencies: (a) synthetic signal generated according
to Eq. (9) with parameters given in Table 1; (c) original(dashed line) and modeled (solid line) amplitude envelopes(curves overlap almost perfectly); and (e) resynthesized signal based on the estimated model. Subplots (b), (d), and (f)refer to a modiﬁed signal generated with parameters shownin Table 1 but with phases replaced by θ1 = θ2 = π/2.

equal frequency: (a) synthetic signal generated according to
Eq. (9) with parameters given in Table 2; (c) original (dashedline) and modeled (solid line) amplitude envelopes; and (e)resynthesized signal based on the estimated model. Subplots (b), (d), and (f) refer to another signal generated withparameters shown in Table 2 but with phases replaced byθ1 = θ2 = π/2.

4.3Modeling of noisy responses

In simulation case (C), in order to verify the FZ-ARMA modeling when dealing with noisy signals, we contaminate the impulse responses shown in the plots (a) and (b) of Fig. 3 withzero-mean additive white Gaussian noise. In this example thevariance of the noise is chosen to produce a signal-to-noise ratio (SNR) of -5 dB in the beginning of the signal. Of course,
the local SNR decreases towards the end of the signal.
The results are displayed in Fig. 5, which follows the samestructure as the previous ﬁgures. Now, by looking at subplots(c) and (d) of Fig. 5 it can be seen that the envelopes of themodeled signals (solid lines) diﬀer substantially from those of
the noisy signals (dashed lines). Moreover, the resynthesized
signals based on the computed models, shown in subplots (e)and (f), are free of visible noise and follow closely their corresponding clean versions, which are depicted in the subplots(a) and (b) of Fig. 3.
The highly successful result of reducing the additive noise insimulation case (C) can be understood when considering thefrequency zooming to a narrow band around the modal frequencies of interest, whereby SNR is improved by the zooming
ratio, i.e., by 10 log10 200 = 23 dB in this case. The low-order
ARMA(4,4) modeling further reduces the inﬂuence of noisedue to good correlation with the modal signals only.
4.4

Envelope modeling of non-LTI system

A primary assumption when applying FZ-ARMA or any LTI
system modeling is that the frequencies of modes do notchange within the duration of the analyzed segment. If this requirement cannot be satisﬁed, e.g., in strongly plucked string
instrument tones having initial pitch shifting [28], the envelope behavior of the target signal can still be modeled. A
straightforward way, if the frequency trajectory of the pitch
shift is know, is to resample the signal so that the shift iseliminated. Another way is to apply FZ-ARMA modeling but
adopting higher orders for the numerator and denominator sothat this can capture the eﬀect of frequency shift.
Alternatively, one can compute an ARMA model for the envelope of a modulated and decimated signal (FZ-ENV-ARMA).
In that way the envelope behavior can be approximated with

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

6AR/ARMA MODELING

0.5
normalized time

1

1

(f)

0
−1
0

0.5
normalized time

−40
0

1

amplitude(d)0.5
normalized time

1

1

(e)

0
−1
0

0.5
normalized time

−1
0

1

−20

amplitude

−40
0

0

1

Fig. 5: Case (C), amplitude beating of two modes in noise:(a) synthetic signal generated according to Eq. (9) with parameters given in Table 1 immersed in white noise; (c) original(dashed line) and modeled (solid line) amplitude envelopes;and (e) resynthesized signal based on the estimated model.
Subplots (b), (d), and (f) refer to a modiﬁed signal generated
with parameters shown in Table 1 but with phases replacedby θ1 = θ2 = π/2.

a lower model order. Simulation case (D), an example thatcompares the standard FZ-ARMA modeling against the FZENV-ARMA, is shown in Fig. 6. The test signal plotted insubplot (a) is a version of that plotted in Fig. 3, but now thefrequencies of the modes start 50 Hz above the values indicated in Table 1 and then converge exponentially with a time
constant of 100 ms to the nominal values. The plots on theleft column show the original and modeled envelopes for different FZ-ENV-ARMA model orders. The plots on the rightcolumn do the same but employing FZ-ARMA models.
To resynthesize a changing pitch signal based on the FZ-ENVARMA computed model, it is necessary to estimate its pitch
behavior. Then, after obtaining a model for the amplitudeenvelope, a frequency modulation corresponding to the original frequency shift should be employed during synthesis. Fora direct FZ-ARMA modeling this is not needed if the estimation is capable of capturing the given behavior of the shifting
modal frequencies.
It can be veriﬁed from Fig. 6 that, in constrast to what happens to the FZ-ARMA modeling, increasing the model orderin the FZ-ENV-ARMA does not essentially help to improvethe model ﬁt, since the inherent phase relations of the original
signal have been lost in the computation of the envelope thatis used as a target. Nevertheless, for low-order modeling, FZENV-ARMA yields better envelope ﬁt than the equal order
FZ-ARMA modeling.
If the response of a target system is highly complex in modedensity, such as a room response at medium to high frequences, a detailed modal description may not be feasible

0

0

0.5
normalized time0(b)envelope dB

−20

0.5
normalized time

(a)

−20
−40
0

0.5
normalized time01envelope dB

(c)

−5
0

1

(d)

−20
−40
0

0.5
normalized time01envelope dB

1envelope dB

0.5
normalized timeenvelope dB

−5
0

0

envelope dB

0

(b)

5

envelope dB

amplitude

(a)

5

0

amplitudeenvelope dB

amplitude

KARJALAINEN ET AL.

(f)

−20
−40
0

0.5
normalized time

1

1
(c)

−20
−40
0

0.5
normalized time

0

1
(e)

−20
−40
0

0.5
normalized time

0

1
(g)

−20
−40
0

0.5
normalized time1

Fig. 6: Case (D), amplitude beating with pitch-shift: (a)synthetic signal generated according to Eq. (9) with parameters given in Table 1 but with a pitch shift. The plots on
the left column show the original (dashed-line) and modeled (solid-line) amplitude envelopes for the FZ-ENV-ARMAmodels of orders: (b) ARMA(4,4), (d) ARMA(8,8), and (f)ARMA(12,12). Subplots (c), (e), and (g) on the right column
show the same cases but for the FZ-ARMA modeling withcorresponding model orders.

or even desired. In such cases the envelope behavior can be
represented simply by ﬁtting a lower order model to the decaying envelope in a desired frequency range by FZ-ENV-ARMA
techniques. This can be useful in decay time estimation6.
Simulation case (E) in Fig. 7 depicts the decay envelope of an
example room response for the octave band 1–2 kHz and a
related envelope curve ﬁtting by low-order FZ-ARMA modeling.
Another form of non-LTI behavior are nonlinearities. A small
degree of nonlinearity in a system can be accepted, and evenquite severe deviation from linear can be tolerated if we accept
the fact that the parameters are then signal-dependent, forexample dependent on the level of a signal.
4.5

Selection of zooming parameters

The choice of the FZ-ARMA parameters, i.e., Ωm , Kzoom ,and the ARMA orders N and P , depends on several factors.
Considering ﬁrst the zoom factor, it can be said that thelarger Kzoom is, the higher the frequency resolution. This
6
There exist many methods, however, that are better motivated for example in reverberation time RT60 estimation [13].

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

7

KARJALAINEN ET AL.

AR/ARMA MODELING

0
(a)

-40
-60

measured

-80

0

0.2

0.4

0

0

modeledT60 ≈ 0.35 s

noise floor
0.6time/s-101level/dBlevel/dB-20

Fig. 7: Case (E), estimation of reverberant response decay

0.2

-20
-30

0.4

-40

rate and T60 by modeling the decay envelope using low-order
FZ-ENV-ARMA model. Measured response was bandpass
ﬁltered (1–2 kHz), absolute value of Hilbert envelope taken
and decimated by 500, modeled by ﬁlter orders N = 6, P = 6,
and ﬁnally the largest positive (real) pole corresponding to themain decay component was identiﬁed for slope estimation.

-50
0

0.6
50
frequency/Hz

100

150

time/s200(b)favors cases in which the modes are densely distributed infrequency. On the other hand, high values for Kzoom imply a more demanding signal decimation procedure and lesssamples available for modeling in the decimated signal.

The order of an ARMA model will be dependent on the number of modes associated with each resonance group. Experiments on two-mode resonances reveal that adopting an
ARMA(4,4) model in general yields satisfactory results forsuch cases. Better modeling accuracy can be achieved by increasing the order, although the result may not be physicallyinterpretable for a two-mode case. High-order analysis also
rises the probability of ending up with an unstable model.
Note that in FZ-ARMA the modulation frequency fm mustcorrespond to the lower edge of the focusing range and thezooming factor Kzoom , in relation to sampling rate, determines the zooming bandwith.

5

CASE STUDIES IN AUDIO APPLICATIONS

AR/ARMA modeling has many applications in modern audio signal processing. Linear and time-invariant models canbe applied for example in room acoustics, sound synthesis,
and audio reproduction. Based on the previous theoretical
overview and examples, in this section we will take examples
from these domains to study the feasibility of the methods,
particularly of the FZ-ARMA technique, in several audio applications.

0

0
-10
level/dB

A possible strategy is to deﬁne a minimum fs,zoom beforehand
and then derive Kzoom . For instance, the criterion may bebased on the number of samples that remain in the decimatedsignal. Another natural choice when there are relatively isolated modes or mode groups is to select the frequency rangeof focusing to cover such a group and its vicinity until neighboring modes or groups start to have inﬂuence. It is recommended to to choose the range of focus such that resonancepeaks are not placed at the ends of the subband. As a rule ofthumb, a suitable choice is to set fm = fr − fs,zoom /4, whichbrings the resonance frequency fr in the middlepoint of thedecimated frequency range.

0.2

-20
-30

0.4

-40
-50
0

0.6
50
frequency/Hz

100

150

time/s200

Fig. 8: Case (F): (a) cumulative decay spectrum of a measured room response below 220 Hz, for time span 0.0 . . . 1.0 s;(b) cumulative decay spectrum computed from the impulse
response of an AR model (P = 80) for the room response of
Fig. 8(a). The original sampling rate 44100 Hz was decimatedto 400 Hz before AR modeling. Level is limited to -50 dB.
5.1

Modeling of room impulse responses

A challenging application for AR/ARMA modeling is to ﬁnd
compact but perceptually valid approximations for measured(or computed) room impulse responses [12]. This is neededfor example in modal analysis of rooms at low frequencies,
artiﬁcial reverb designs, or equalization of loudspeaker-room
responses.
As case study (F), an analysis of low-frequency modal behavior of a room impulse response is carried out using diﬀerent
AR and FZ-ARMA methods. The room has approximate dimensions of 5.5 x 6.5 x 2.7 m3 . Figure 8(a) describes thetime-frequency behavior (cumulative decay spectrum) for frequencies below 220 Hz as computed from a measured room
impulse response. The room shows particularly intense andlong modal decays around 45 Hz.
A straightforward AR modeling of the room impulse response
below 220 Hz using linear prediction yields fairly accurate results when the all-pole ﬁlter order P is about 100 or above

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

8

KARJALAINEN ET AL.

1

-40

46Hz

0.5

55Hz

-0.5

0.5

1

-1

time/s
(c)

-40

0.5

0

-1

-40
-60

0.5

1
time/s

0.15

0.2

(b)

0.5
0
-0.5

0
Real part

-40

0

0.05

0.1

0.15

0.2

0

(f)

-1

0.25

0.25

time/s

-1
0

0.1

-20

-60

1

level/dB-20

0Real part

1
Imaginary partlevel/dB

(e)

0.05

0

-0.5

time/s
0

0

(d)

0.5

1

-40

time/s

-1
0

-20

-60

1

level/dB-20

0Real part

1
Imaginary partlevel/dB

0

-60

37Hz

0

-1
0

(a)

(b)

level/dB

-20

-60

0

(a)
Imaginary partlevel/dB

0

AR/ARMA MODELING

Fig. 9: Case (G), envelope match for the decay of modalgroup between 33 and 58 Hz in the room response of Fig.
8 with diﬀerent FZ-ARMA orders: (a) ARMA(4,6), (c)ARMA(4,20), (e)ARMA(6,60). Solid line = modeled, dotted line = original. Poles (x) and zeros (o) in the decimatedz-domain are plotted in (b), (d), and (f), respectively. (Somezeros fall outside the plotting range.)

for the low-frequency range. Fig. 8(b) shows the model response decay plot when P = 80. The original sample rate of
44100 Hz was decimated by a factor of 110 before AR modeling. A comparison with Fig. 8(a) reveals that the decay times
of prominent modes are quite well modeled but many weakermodes are too short or too damped due to insuﬃcient modelorder.
Direct ARMA modeling by Steiglitz-McBride method yields
a better time domain ﬁt with a given denominator order thanthe corresponding AR model. For example using numerator
order of N = 30 and denominator P = 100 worked fairly
well for the room response above, although in many cases the
Steiglitz-McBride algorithm gives an unstable result alreadywith such moderate ﬁlter orders.
FZ-ARMA is a powerful method for accurate modeling ofmodal behavior in a limited frequency range. Figure 9 depicts
modeling results of the prominent modal region 33 . . . 58 Hz
in the response of Fig. 8(a). The region includes three major
modes at frequencies of 37, 46, and 55 Hz. Figure 9 shows
the decay envelope of the modal region for the original signal(dashed line) and as a result of applying the Steiglitz-McBride
method of diﬀerent orders (solid line). Increasing the ﬁlter
order improves the envelope ﬁt, but ﬁnally it may start tomodel the background noise envelope. The pole-zero plots on
the right hand side indicate that for an order of P = 6 thepoles correspond to the three modes, while for higher orders
there are extra poles and it is not easy to associate them withthe modes.

-20
-40
-60

1

(c)

0

0.05

0.1

0.15

0.2

0.25

time/sFig. 10: Case (H), ﬁtting to room response within criticalband 920 . . . 1080 Hz by FZ-ARMA of orders (a) ARMA(6,6),(b) ARMA(6,20), (c) ARMA(6,60), described by amplitudeenvelopes of original (dashed line) and modeled (solid line)
responses.
At higher frequencies, above the critical frequency (Schroeder
frequency) [29] of the room, the modal behavior is diﬀuse, i.e.,the modal density is high and modes overlap in frequency.
Full audio range AR and ARMA modeling is diﬃcult, if not
impossible. However, it is possible to apply the FZ-ARMA
analysis to a narrow frequency band of a reverberant response.
Figure 10 describes a ﬁtting to the room response studied
within a critical band at 1 kHz (920 . . . 1080 Hz) by diﬀerent
model orders. With the highest model order P = 60, envelope
ﬁtting is good for the ﬁrst 250 ms and for about a 40 dB
dynamic range.
5.2

Loudspeaker/room equalization

Equalization of a louspeaker response or a loudspeakerroom reproduction chain means correcting the system response closer to desired perceptual or technical criteria. MA
and ARMA modeling have been reported in the literature
in several forms for loudspeaker and in-situ frequency response equalization, both in on-line and oﬀ-line formulations[30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 17, 40, 41, 42, 43, 23].
Equalization of the free-ﬁeld magnitude response (possibly
including the phase response) of a loudspeaker by DSP can be
carried out using many known techniques. For highest quality
loudspeakers there is hardly any need to improve its free-ﬁeld
response, but the loudspeaker-room combination may beneﬁt
greatly from proper equalization.
The combined task of loudspeaker and room equalization is
also demanding since it is essentially a problem of ﬁnding
a perceptually optimal time-frequency equalization, insteadof simple ﬂattening of the magnitude spectrum and/or lin-

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

9KARJALAINEN ET AL.

AR/ARMA MODELING

(a)

(b)

0

level/dB

-10

0.2

-20
-30

0.4

-40

0

0
-10
level/dB

0

0.2

-20
-30

0.4

-40

-50

0.6

0

50
frequency/Hz

100

150

time/s

-50

0.6

0

50
frequency/Hz

200

100

150

time/s

200

(c)
0.9
(d)
0.8
0.7
0.6

0

0
level/dB

0.5
0.4

-10

0.2

-20

0.3

-30

0.4

0.2

-40
-50
0

0.6
50
frequency/Hz

100

150

time/s

200

0.1
0
0

50

100

150

200

frequency/Hz

Fig. 11: Case (I), cumulative decay spectra for equalization of loudspeaker plus synthetic modes: (a) original response plus
ﬁve modes at frequencies 50, 55, 100, 130, and 180 Hz, (b) modes fully damped, (c) modal decay time (60 dB) equalized to 250
ms, and (d) mode detection function G(Ω) for the original response. Decay levels are limited to -50 dB.

earization of phase. There seems to be quite a common misunderstanding that just ﬂattening the response, at least at
low frequencies where it might be technically possible, would
be an ideal solution. A better strategy is to improve thebalance of overall acoustical parameters, particularly of thereverberation time. As discussed in [44], this can be done bycontrolling the decay times of individual modes at low frequencies, typically below 200 Hz, to match the reverberation
time at mid frequencies. This is called modal equalization. It
may be followed by a traditional correction of the envelopeof magnitude response. The need for such active correctionof room acoustics is particularly prominent around 100 Hz
even in spaces designed for listening purposes, such as audio
monitoring rooms [45].

In [44] we proposed a method for modal equalization. In thepresent paper we suggest another technique to realize modal
equalization, optionally combined with magnitude envelope
correction. The general framework of modal equalization has
been discussed in detail in the previous paper. A brief description of the procedure is:

1. Measure the combined loudspeaker plus room impulse
response in the listening position of interest. Any modern technique for reliable response measurement can beapplied.
2. Analyze the average reverberation time RT60 at mid frequencies, for example between 500 Hz . . . 2 kHz.
3. Determine an upper limit of modal decay time as a function of frequency for the low-frequency range, typicallybelow 200 Hz. This value can be allowed to grow slightly
toward lowest frequencies [46, 47], for example linearlyby 0.2 s when the frequency decreases from 300 Hz to50 Hz.
4. Find the modes that need equalization, i.e., those thathave a longer decay time than the upper limit deﬁned
above. If the magnitude level of a mode is so low that its
tail remains below a given level, it does not need modal
equalization even when its decay time is longer than theupper limit. Estimate modal parameter values for thesemodes, particularly the modal frequency and the decay
time constant, and compute the angles and radii of thecorresponding poles.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

10KARJALAINEN ET AL.

AR/ARMA MODELING

5. Design a correction ﬁlter for each mode requiring equalization so that the ﬁlter shortens the decay time to meetthe upper limit criteria speciﬁed in step 3. This means
canceling the estimated pole pair, which represents amode with a long decay time, by a zero pair, and replacing it with a new pole pair having the desired decay
time. This can be done with an IIR ﬁlter [44](1 − rejθ z−1 )(1 − re−jθ z−1 )(1 − rc ejθc z−1 )(1 − rc e−jθc z−1 )

0

0

(10)

-10

where r and rc are the (complex conjugate) pole radii ofthe original decay and the corrected decay, respectively,and θ and θc are the corresponding pole angles.

level/dBHc (z) =

(a)

0.2

-20
-30

0.4

-40

6. Compute steps 4 and 5 either in a batch mode, i.e., inparallel for each mode to be equalized, or iteratively sothat modes are equalized one by one, starting from themost prominent one and returning to step 4, to be applied to the result of the previous equalization. The process is teminated when all remaining modes meet the
decay time criteria or when a preset upper limit of correctable modes has been reached.

-50
0

0.6
10

20

30
frequency/Hz

40

50

60

70

time/s80(b)

7. Traditional magnitude equalization can be applied to theresult of modal equalization, if needed, by any method
or technique appropriate.

In the equalization cases below, mode ﬁnding and parameter
estimation are carried out iteratively in the following way:
1. Compute a function that can be used robustly to ﬁnd
the most prominent modes and their frequencies. Thiscan be done in diﬀerent ways, for example directly by
AR or ARMA analysis and ﬁnding the poles with largest
radii. Because the selection of proper model order can beproblematic, we have ﬁrst applied here a separate mode
detection function
G(Ω) =|H(Ω)| max(0, D(arg(H(Ω)))(11)

where H(Ω) is the Fourier transform of the measured
response, Ω is the normalized angular frequency (angle
in the z-plane), and D is a diﬀerentiation operator (inthe frequency domain). An example of G(Ω) functionis plotted in Fig. 11(d). Positive peaks indicate strong
modes that may need decay time equalization. Note that
this function combines both the magnitude level and the
decay time (through phase derivative) information.
2. Find the highest peak position Ωp of the detection function G(Ω) as the best candidate for modal equalization.
Run AR analysis of predeﬁned order (here we appliedorders 50 . . . 70 for the frequency range below 220 Hz)
on the minimum phase version of the target response toﬁnd poles and select the pole closest to the point ejΩp

0

0
-10
level/dB

In this context we are only interested in step 4 as a part ofbatch or iterative analysis. All other steps follow the general scheme described in [44], where the mode search andthe decay time estimation were based on a time-frequency
representation and ﬁtting of a logarithmic decay plus background noise model using nonlinear optimization. While thepreviously proposed method is found robust for modes thatare separated well enough, strongly overlapping or multiple
modes with closely similar frequency are an inherent diﬃculty of that method. Since AR/ARMA models search fora global optimum and don’t try to separate modes, they arepotentially a better alternative in such cases.

0.2

-20
-30

0.4

-40
-50
0

0.6
10

20

30
frequency/Hz

40

50

60

70

time/s80

Fig. 12: Case (J), equalization of a single mode at 46 Hz of
the room analyzed in Fig. 8: (a) original response, (b) aftermode equalization.

on the unit circle. This pole and its complex conjugate
now represent the most prominent mode.
3. If the decay time of the mode is below the upper limit
allowed and the value of G(Ωp ) is below a threshold experimentally determined, go to ﬁnalize the process instep 4. If not, design a second-order modal correction
ﬁlter of type Eq. (10) to change the modal decay time toa desired value below the upper limit. Apply this to the
response to be equalized and use the result when goingback to iterate from step 1.
4. Finally, collect the correction ﬁlters into a cascaded ﬁlter
which is now the modal equalizer for the system.
A simulated modal equalization, case (I), is illustrated inFig. 11. A loudspeaker impulse response is ﬁltered to add
ﬁve simulated modes at frequencies 50, 55, 100, 130, and 180
Hz, with 60 dB modal decay times of 1.4, 0.8, 1.0, 0.8, and 0.7
seconds, respectively. The cumulative decay spectrum of this
synthetic response is shown in Fig. 11(a). Figure 11(b) provesthat the eﬀect of the modes can be cancelled out almost perfectly, leaving the loudspeaker response only, by moving the

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

11

KARJALAINEN ET AL.

AR/ARMA MODELING

pole radii to correspond to a very short decay time (about 60ms) using the procedure described above. In Figure 11(c) the
result of modal equalization depicted is more appropriate forreal room conditions. The decay time of each mode is equalized to 250 ms. The two nearby modes partly overlapping,at 50 and 55 Hz, do not cause any diﬃculties, and the modal
equalization works almost perfectly.

0

0
-10
level/dB

In case (J) of Fig. 12 the most prominent single mode at 46 Hz
is equalized by shortening the decay time from a value above
1 second to about 300 ms using the algorithm described above
and limiting the search for modes to only one. In Fig. 12(b)the originally problematic mode decays now clearly faster.
Furthermore, the equalized response up to 80 Hz has muchsmoother shape since the modal equalization also aﬀects themagnitude spectrum. However, the decay times of some othermodes remain long.

(a)

0.2

-20
-30

0.4

-40
-50

0.6

0

50
frequency/Hz

Multi-mode equalization of the same room, case (K), is shownin Fig. 13. The room is the same as the one analyzed incase (F), Fig. 8. The procedure described above is iterated
100 times, yelding 100 second-order correction ﬁlter sections,to shorten the mode decay times. The cumulative decay
spectrum of the resulting equalized response is illustrated inFig. 13(b). The target value for equalized modal decay time
(60 dB) has been 150 ms.

0

-10

0.2

-20
-30

0.4

-40
-50

The experiments with string and bell sounds here have been
inﬂuenced by techniques published in references [50, 16, 51,
52, 53, 54]. Particularly the application of Prony’s method
in [16] is aiming to the same goal. We have found that the
Steiglitz-McBride iteration as a part of FZ-ARMA yields typically better results than Prony’s method.
As a simulated case to study the required model orders wegenerate a synthetic guitar tone using the dual polarization
model in [55]. Each partial has two modes with known parameters, i.e., resonance frequencies and time-constants of the
exponentially decaying envelope. When dealing with real instrument tones, for example guitar tones recorded in an anechoic chamber, it is necessary to estimate the frequencies of

0.6

0

50
frequency/Hz

100

150

time/s

200

5
(c)
0

-5

level/dB

Modeling of string instruments

An appealing application of the FZ-ARMA scheme is themodeling of musical instruments. For instance in sound synthesis of string instruments through any parametric approach,such as digital waveguide modeling [48], one may need to obtain information related to the partial frequencies and theirrespective decay rates from measured signals [49, 26]. The
FZ-ARMA modeling turns out to be a suitable tool for suchmodal analysis purposes.

200

0

The ﬁnal step, i.e., smooth envelope equalization of the magnitude response, is not discussed here since many known techniques could be applied to equalize the magnitude response.
An interesting choice is, however, to integrate the magnitude equalization phase together with the AR/ARMA modal
equalization process.
5.3

150

time/s(b)level/dB

In this case the result is not as perfect as in the synthetic orsingle mode case. There is about 10 dB of fast decay in the beginning, as shown by backward integrated plots in Fig. 13(c),and thereafter the decay rate follows the original one. Although the ideal shortening of the decay time is not achieved
precisely, it already makes sound reproduction in the room
more balanced in the terms of reverberation. Furthermore,
the equalization procedure can be improved by careful adjustment of the details.

100

-10

-15

-20

-25

-30

0

0.1

0.2

0.3
time/s0.4

0.5

Fig. 13: Case (K), equalization of the room response analyzed in Fig. 8: (a) original response, (b) after mode equalization, and (c) backward integrated energy decay for the original (dashed line) and equalized (solid line) responses.
the resonance peaks beforehand and select subbands aroundthem in order to proceed with the FZ-ARMA modeling.
In simulations it was found that already low-order frequencyfocusing models of ARMA(2,4), i.e., with two pole pairs, suf-

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

12

KARJALAINEN ET AL.

AR/ARMA MODELING

magnitude [dB]

−20
−40
0

−60

magnitude [dB]

−20

0

−40
−60
−80

2
2

4

0

−100

8

4

4

4
6

2

2
6

6
8

6

time [s]1010partial index

partial index

8time [s]−20
magnitude [dB]magnitude [dB]

0
−20
−40
0

−60
2
2

4

−40
−60
−80
0

−100

8

6
10

4

4
6

4
6

2

2
6
8

time [s]partial index

108time [s]partial index

Fig. 14: Case (L), analysis of guitar tone partials. The plot
on the top shows the magnitude envelope of the ﬁrst 10 partials of a A2 guitar tone while the plot on the bottom depicts
the corresponding envelope estimates using the FZ-ARMAmodeling.

Fig. 15: Case (M), analysis of piano tone partials. The plot
on the top shows the magnitude envelope of the ﬁrst 10 partials of a D4 piano tone while the plot on the bottom depicts
the corresponding envelope estimates using the FZ-ARMAmodeling.

harmonics is approximated well in this case.
ﬁce to estimate the mode frequencies properly. To matchaccurately the decay times and forms of the modes, a little
higher model order is required. ARMA(4,6) worked alreadywell.
For harmonic signals with non-varying pitch, a pitch detector
could be used to select both the down-modulation frequencies fm and the sub-band bandwidth fs,zoom . A more generalapproach should also work for inharmonic signals and in thiscase, running a peak-picking algorithm over the magnitude
spectrum of the initial part of the tone seems to be more
appropriate. In addition, the estimation of the number ofresonant modes per partial may be based on a priori information, e.g., the number of strings per note vibrating in a
piano tone.
In an example, case (L), FZ-ARMA modeling is performed
on a guitar tone. The signal was recorded in an anechoic
chamber and corresponds to plucking the ﬁfth string (A2,f0 ≈ 110 Hz). The results of modeling the ﬁrst 10 harmonics
fi , i = 1, . . . , 9, of this tone are shown in Fig. 14. Parameters used in FZ-ARMA modeling were fm = fi − 0.5fs,zoom ,Kzoom = 401, P = 6, and N = 6. The temporal structure of

As another example, case (M), Figure 15 plots the original
and estimated amplitude envelope curves for the ﬁrst 10 harmonics fi of a D4 (f0 ≈ 298 Hz) piano tone. The analysis
was performed using the FZ-ARMA modeling with parameters fm = fi − 0.5fs,zoom , Kzoom = 247, P = 6, and N = 6.
This order is able to capture the main structure of temporal
evolution in partial decays, but higher orders are required forhigh-precision modeling because for a three-string piano notethe degree of freedoms corresponds theoretically to 3·2·2 = 12
poles.
The string models obtained above can be used directly forsource-ﬁlter (= subtractive) synthesis, whereby the total ﬁlter order for single string synthesis is N = P = 50 . . . 500
depending on fundamental frequency and desired quality ofsynthesis. Since digital waveguide models [48] are attractivefrom a computational eﬃciency point of view, it would be useful to map the modal parameters of partials to a loop ﬁlterof a digital waveguide. To gain advantage by the waveguideformulation, a relatively low-order loop ﬁlter is needed thatapproximates the common properties of the separate modes,
thus leading to a nonlinear optimization task.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

13

KARJALAINEN ET AL.

AR/ARMA MODELING

level/dB

(a)

0
-1
0

0.5

1

1.5

2

time/s

100

60

40
0

0.5

1

1.5

2

time/s

2.5

60

(b)

80

50

30

2.5

(b)
level/dBlevel/dB

60

(a)

1

50
40

40
0

1000

2000

3000

4000

5000

6000

7000

8000

level/dBlevel/dB

60level/dB

40
1270
20

1280

1290

1300

1310

1320

1330

1340

0

0

0.5

1

1.5

2

time/s1

1.5

2

time/s

2.5
(c)

40
0

0.5

1

1.5

2

time/s2.5

Fig. 17: Case (N), analysis and modeling of beating envelope in a bell sound mode with diﬀerent model orders: (a)ARMA(0, 4), (b) ARMA(40, 4), and (c) ARMA(2, 6).

2.5

6

Fig. 16: Acoustical behavior of a small bell: (a) recorded
time domain signal, (b) magnitude spectrum, (c) magnitude
spectrum in the modal region around 1310 Hz, (d) decay envelope of the 1310 Hz modal group.

5.4

0.5

50

30

1350 freq/Hz
(d)

0

60

(c)

80

-20

30

9000 freq/Hz

100

Modeling of bell sounds

As another musical instrument, case (N), the analysis andmodeling of bell sounds is presented. A characteristic featureof bell sounds is that they are composed of an inharmonicset of partials [2], such as the one described by magnitude
spectrum in Fig. 16(b). Each partial is a decaying sinusoid
that, in a closer inspection, Fig. 16(c), turns out to be apair or a group of modes very closely located in frequency.
This leads to perceptually noticeable beating. In this case themodal group consists primarily of two modes with a frequency
diﬀerence of about 2.5 Hz.
FZ-ARMA is an excellent method for analyzing the modalgroups of bell sounds. Figure 17 shows the envelope match
obtained with three diﬀerent FZ-ARMA orders for the 1310
Hz modal group. The zooming factor Kzoom is 200 in eachcase. In case (a) the orders are N = 0 and P = 4. Twopole pairs should in principle be suﬃcient for a double mode,
but this all-pole (AR) case with N = 0 does not allow proper
phase matching, and thus the overall match remains poor.
For ARMA orders N = 40 and P = 4 in (b) the relatively
high number of zeros allows for good match with just two polepairs. The same can be achieved in case (c) with orders N = 2
to P = 6, i.e., by adding an extra pole pair and keeping thenumber of zeros minimal. For all resonances up to 10 kHz forthis bell sound, ﬁlter orders of N = 2 . . . 4 and P = 4 . . . 6 are
suﬃcient for good modal decay matching so that a parallelﬁlter, composed of modal group ﬁlters with a total order ofabout N = 40 and P = 50 can implement an eﬃcient andhigh-quality synthesizer for the bell sound for the sampling
rate of 22050 Hz.

DISCUSSION AND CONCLUSIONSIn this paper we have studied the modeling of acoustic andaudio system responses that exhibit resonant and reverberant
properties. Particularly the AR and ARMA modeling techniques are investigated to obtain eﬃcient all-pole or pole-zeroﬁlters. Such modeling, if accurate enough and computationally inexpensive, ﬁnds applications in solving many audiooriented problems.
The ﬁrst part of the paper is a non-theoretical overview of AR
and ARMA modeling methods to demonstrate their inherent
properties and limitations.
A speciﬁc interest of this study has been the modeling methods that can yield good temporal match to a given target
response and high frequency resolution, often at the sametime. Based on earlier studies, primarily on applying Prony’s
method to subbands, we show that frequency-zooming ARMA(FZ-ARMA) based on the Steiglitz-McBride iteration is a
powerful technique for high-resolution modeling in subbands.
Simulation examples demonstrate the ability of this approach
to model complex modal and reverberant behaviors.
These methods are then used to solve problems of practicalinterest in audio applications, covering room impulse responsemodeling, inverse modeling for equalization of loudspeakerroom responses, and modeling of musical instrument sounds.
Many general and problem-speciﬁc questions remain for further studies. Robust selection of model orders, stability andnumerical sensitivity, eﬃcient and robust implementation of
the obtained ﬁlters for applications, adaptive formulation andon-line calibration of the models, as well as ﬁnding new applications, call for more investigation.
Matlab code, examples, and sound demonstrations are available at: www.acoustics.hut.fi/demos/AR ARMA/

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

14

KARJALAINEN ET AL.
7

ACKNOWLEDGMENTS

AR/ARMA MODELING[17] J. Valakas, M. Paraskevas, and J. Mourjopoulos, “PoleZero Loudspeaker Response Analysis by Bandsplitting
and Prony Modeling,” in Preprint 4214, AES 100th
Convention, (Copenhagen), 1996 May.

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 Paulo Esquef has
been ﬁnanced by the Brazilian National Council for Scientiﬁc
and Technological Development (CNPq) and the Academy ofFinland, projects “Sound Source Modeling” and “Technology
for Audio and Speech Processing”.

[18] J. D. Markel and J. H. G. Gray, Linear Prediction of
Speech Signals. Berlin: Springer-Verlag, 1976.

8[20] T. W. Parks and C. S. Burrus, Digital Filter Design.
New York: John Wiley & Sons, Inc., 1987.

REFERENCES[19] L. R. Rabiner and R. W. Schafer, Digital Processing of
Speech Signals. Englewood Cliﬀs, New Jersey: PrenticeHall, Inc., 1978.

[1] A. D. Pierce, Acoustics, An Introduction to Its Physical Principles and Applications. Woodbury, New York:
Acoustical Society of America, 1989.

[21] L. Ljung, System Identiﬁcation — Theory for the User.
Upper Saddle River, N. J.: Prentice-Hall, 1999.

[2] N. H. Fletcher and T. D. Rossing, The Physics of MusicalInstruments. New York: Springer-Verlag, 1991.

[22] A. Härmä, M. Karjalainen, L. Savioja, V. Välimäki,
U. K. Laine, and J. Huopaniemi, “Frequency-Warped
Signal Processing for Audio Applications,” J. Audio
Eng. Soc., vol. 48, no. 11, pp. 1011–1031, 2000 Nov.

[3] M. Colloms, High Performance Loudspeakers. 5th Edition. New York: Wiley, 1997.
[4] S. M. Kay, Fundamentals of Statistical Signal Processing: Volume I: Estimation Theory. Englewood Cliﬀs,
N.J.: Prentice-Hall, 1993.
[5] A. V. Oppenheim, A. Willsky, and I. Young, Signals andSystems. Englewood Cliﬀs, New Jersey: Prentice-Hall,
1983.
[6] S. Haykin, Adaptive Filter Theory. Upper Saddle River,
New Jersey: Prentice Hall, 1996.
[7] M. H. Hayes, Statistical Digital Signal Processing andModeling. Upper Saddle River, N.J.: Prentice Hall, 1996.
[8] MathWorks, Inc., MATLAB Signal Processing Toolbox,
User’s Guide, 2001.
[9] MathWorks, Inc., MATLAB System Identiﬁcation Toolbox, User’s Guide, 2001.
[10] L. B. Jackson, Digital Filters and Signal Processing.
Boston: Kluwer Academic Publ., 1989.
[11] K. Steiglitz, A Digital Signal Processing Primer. Menlo
Park: Addison-Wesley, 1996.
[12] Y. Haneda, S. Makino, and Y. Kaneda, “Common
Acoustical Pole and Zero Modeling of Room Transfer
Functions,” IEEE Trans. Speech and Audio Processing(SAP), vol. 2, no. 2, pp. 320–328, 1992 April.
[13] M. Karjalainen, P. Antsalo, A. Mäkivirta, T. Peltonen,and V. Välimäki, “Estimation of Modal Decay Parameters from Noisy Response Measurements,” in Preprint
5290 AES 110th Convention, (Amsterdam, The Netherlands), 2001 May.
[14] J. D. Bunton and R. H. Small, “Cumulative Spectra,
Tone Bursts and Applications,” J. Audio Eng. Soc.,
vol. 30, no. 6, pp. 386–395, 1982 June.
[15] J. Makhoul, “Spectral Linear Prediction: Properties
and Applications,” IEEE Trans ASSP, vol. 23, no. 5,
pp. 283–296, 1975.
[16] J. Laroche, “A New Analysis/Synthesis System of Musical Signals Using Prony’s Method – Application to
Heavily Damped Percussive Sounds,” in Proc. IEEE
Int. Conf. Acoustics, Speech, and Signal Processing,pp. 2053–2056, vol. 3, 1989.

[23] T. Paatero and M. Karjalainen, “Kautz Filters and Generalized Frequency Resolution — Theory and Audio Applications,” in Preprint 5378 AES 110th Convention,(Amsterdam, The Netherlands), 2001 May.
[24] L. W. P. Biscainho, P. S. Diniz, and P. A. A. Esquef,“ARMA Processes in Sub-Bands With Application to
Audio Restoration,” in Proc. IEEE Int. Symp. Circ.
Syst. (ISCAS 2001), (Sydney, Australia), pp. 157–160,
vol. 2, 2001 May.
[25] J. A. Moorer, “Signal Processing Aspects of Computer
Music: A Survey,” Proc. IEEE, vol. 65, pp. 1108–1137,
1977 Aug.
[26] V. Välimäki and T. Tolonen, “Development and Calibration of a Guitar Synthesizer,” J. Audio Eng. Soc.,
vol. 46, no. 9, pp. 766–778, Sept. 1998.
[27] M. W. Macon, A. McCree, W.-M. Lai, and
V. Viswanathan, “Eﬃcient Analysis/Synthesis of
Percussion Musical Instrument Sounds Using an AllPole Model,” in Proc. IEEE Int. Conf. Acoust. Speech,
and Sig. Proc., pp. 3589–3592, vol. 6, 1998.
[28] T. Tolonen, V. Välimäki, and M. Karjalainen, “Modelingof Tension Modulation Nonlinearity in Plucked Strings,”
IEEE Trans. Speech and Audio Processing, vol. 8, no. 3,
pp. 300–310, 2000, May.
[29] M. R. Schroeder, “Statistical Parameters of the Frequency Response Curves of Large Rooms,” J. Audio
Eng. Soc., vol. 35, no. 5, pp. 299–306, 1987.
[30] J. Mourjopoulos, “Digital Equalization Methods for Audio Systems,” in Preprint 2598, AES 84th Convention,(Paris), 1988 March.
[31] J. Kuriyama and Y. Furukawa, “Adaptive Loudspeaker
System,” J. Audio Eng. Soc., vol. 37, no. 11, 1989 Nov.
[32] S. J. Elliott and P. Nelson, “Multiple-Point Equalization
in a Room Using Adaptive Digital Filters,” J. Audio
Eng. Soc., vol. 37, no. 11, Nov. 1989.
[33] J. Mourjopoulos and M. A. Paraskevas, “Pole and Zero
Modelling of Room Transfer Functions,” J. Sound andVibration, vol. 146, pp. 281–302, 1994 Nov.
[34] R. G. Greenﬁeld, Application of Digital Techniques to
Loudspeaker Equalization.
PhD thesis, Department
of Electronic Systems Engineering University of Essex,
Sept. 1991.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

15

KARJALAINEN ET AL.
[35] R. Genereux, “Adaptive Filters for Loudspeakers andRooms,” in Preprint 3375, AES 98th Convention, (SanFrancisco), 1992 Oct.
[36] F. X. Y. Gao and W. M. Snelgrove, “Adaptive Linearization of a Loudspeaker,” in Preprint 3377, AES 98th Convention, (San Francisco), 1992 Oct.
[37] J. Mourjopoulos, “Digital Equalization of Room Acoustics,” J. Audio Eng. Soc., vol. 42, no. 11, 1994 Nov.
[38] G. Karachalios, D. Tsoukalas, and J. Mourjopoulos,“Multiband Analysis and Equalization of Loudspeaker
Response,” in Preprint 3975, AES 98th Convention,(Paris), 1995 Feb.
[39] S. Salamouris, K. Politopoulos, V. Tsakiris, and J. Mourjopoulos, “Digital System for Loudspeaker and Room
Equalization,” in Preprint 3976, AES 98th Convention,(Paris), 1995 Feb.
[40] Y. Haneda, S. Makino, and Y. Kaneda, “Multiple-Point
Equalization of Room Transfer Functions by Using Common Acoustical Poles,” IEEE Trans. Speech and Audio
Processing (SAP), vol. 5, no. 4, pp. 325–333, 1997 July.
[41] M. Karjalainen, E. Piirilä, A. Järvinen, andJ. Huopaniemi, “Comparison of Loudspeaker Response Equalization Using Warped Digital Filters,” J.
Audio Eng. Soc., vol. 47, no. 1/2, pp. 15–31, 1999.
[42] F. Fontana, “Common Pole Equalization of Small Rooms
Using a Two-Step Real-Time Digital Equalizer,” in Proc.
IEEE Workshop on Applications of Signal Proccessing
to Audio and Acoustics (WASPAA’99), pp. 195–198,
1999.
[43] M. Tyril, J. A. Pedersen, and P. Rubak, “Digital Filters
for Low-Frequency Equalization,” J. Audio Eng. Soc.,
vol. 49, no. 1/2, pp. 36–43, 2001 Jan./Feb.
[44] A. Mäkivirta, P. Antsalo, M. Karjalainen, andV. Välimäki, “Low-Frequency Modal Equalization ofLoudspeaker-Room Responses,” in Preprint 5480, AES
111th Convention, (New York, USA), 2001 Nov/Dec.
[45] A. Mäkivirta and C. Anet, “A Survey Study of InSitu Stereo and Multi-Channel Monitoring Conditions,”in Preprint 5496, AES 111th Convention, (New York),
2001 Nov./Dec.

AR/ARMA MODELING[46] AES Technical Committee on Multichannel and Binaural Audio Technology (TC-MBAT), “Multichannel Surround Sound Systems and Operations,” Technical Document, version 1.5, AES, 2001.
[47] EBU Document Tech. 3276-1998 (Second ed.), “Listening Condition for the Assessment of Sound ProgrammeMaterial: Monophonic and Two-Channel Stereophonic,”
EBU, 1998.
[48] J. O. Smith, “Physical Modeling Using Digital Waveguides,” Computer Music J., vol. 16, no. 4, pp. 74–91,
1992.
[49] V. Välimäki, J. Huopaniemi, M. Karjalainen, andZ. Janosy, “Physical Modeling of Plucked String Instruments with Application to Real-Time Sound Synthesis,”
J. Audio Eng. Soc., vol. 44, no. 5, pp. 331–353, May
1996.
[50] J. O. Smith, Techniques for Digital Filter Design and
System Identiﬁcation with Application to the Violin.
PhD thesis, CCRMA, Dept. of Music, Stanford University, Stanford, CA, 1983 June. Report STAN-M-14.
[51] M. Sandler, “Algorithm for High Precision Root Finding from High Order LPC Models,” IEE Proceedings I
on Communications, Speech and Vision, vol. 138, no. 6,
pp. 596–602, 1991 Dec.
[52] M. Sandler, “New Synthesis Techniques for High-Order
All-Pole Filters,” IEE Proceedings I on Communications, Speech and Vision, vol. 139, no. 1, pp. 9–16, 1992
Feb.
[53] J. Laroche and J.-L. Meillier, “Multichannel Excitation/Filter Modeling of Percussive Sounds with Application to the Piano,” IEEE Trans. Speech and Audio
Proc., vol. 2, no. 2, pp. 329–344, 1994 April.
[54] J. Mackenzie, I. Kale, and G. D. Cain, “Applying Balanced Model Truncation to Sound Analysis/Synthesis
Models,” in Proc. Int. Computer Music Conf., pp. 400–
403, 1995 Sept.
[55] M. Karjalainen, V. Välimäki, and T. Tolonen, “PluckedString Models, from the Karplus-Strong Algorithm toDigital Waveguides and Beyond,” Computer Music
Journal, vol. 22, no. 3, pp. 17–32, 1998.

AES 112TH CONVENTION, MUNICH, GERMANY, 2002 MAY 10–13

16