EP2104374B1

EP2104374B1 - Spatially robust audio precompensation

Info

Publication number: EP2104374B1
Application number: EP08102812A
Authority: EP
Inventors: Lars Johan Brännmark; Mikael Sternad; Anders Ahlén
Original assignee: Dirac Research AB
Current assignee: Dirac Research AB
Priority date: 2008-03-20
Filing date: 2008-03-20
Publication date: 2010-05-05
Anticipated expiration: 2028-03-20
Also published as: EP2104374A1; DE602008001155D1; ATE467316T1

Abstract

A discrete-time audio precompensation filter is designed based on a linear model that describes the dynamic response of a sound generating system at p > 1 listening positions. The filter construction is based on providing information (S2) representative of n non-minimum phase zeros { z i } that are outside of the stability region l z l = 1 in the complex frequency domain. A causal Finite Impulse Response (FIR) filter, of user-specified degree d , having coefficients corresponding to a causal part of a delayed non-causal impulse response is determined (S4) based on the information representative of n non-minimum phase zeros. The resulting precompensation filter is determined (S5) as the product of at least two scalar dynamic systems, represented by an inverse of a characteristic scalar magnitude response (S3) in the frequency domain that represents the power gains at the listening positions, and the causal Finite Impulse Response (FIR) filter designed (S4) to approximately invert only non-minimum phase zeros that can be safely inverted.

Description

TECHNICAL FIELD OF THE INVENTION

The present invention generally concerns digital audio precompensation, and more particularly the design of a digital precompensation filter that generates an input signal to a sound generating system, with the aim of modifying the dynamic response of the compensated system, as measured in several measurement positions.

BACKGROUND OF THE INVENTION

A system for generating or reproducing sound, including amplifiers, cables and loudspeakers, will always affect the spectral properties of the sound, often in unwanted ways. The reverberation of the room where the equipment is placed adds further modifications. Sound reproduction with very high quality can be attained by using matched sets of cables, amplifiers and loudspeakers of the highest quality, but this is cumbersome and very expensive. The increasing computational power of PCs and digital signal processors has introduced new possibilities for modifying the characteristics of a sound generating or sound reproducing system. The dynamic properties of the sound generating system may be measured and modeled by recording its response to known test signals, as well known from the literature. A precompensation filter, R in Fig. 1, is then placed between the original sound source and the audio equipment. The filter is calculated and implemented to compensate for the measured properties of the sound generating system, symbolized by H in Fig. 1. In particular, it is desirable that the phase and amplitude response of the compensated system is close to a prespecified ideal response, symbolized by D in Fig. 1. In other words, it is thus required that the compensated sound reproduction y(t) matches the ideal y_ref(t) to some given degree of accuracy. The pre-distortion generated by the precompensator R cancels the distortion due to the system H, such that the resulting sound reproduction has the sound characteristic of D. Up to the physical limits of the system, it is thus, at least in theory, possible to attain a superior sound quality, without the high cost of using extreme high-end audio equipment. The aim of the design could, for example, be to cancel acoustic resonances caused by imperfectly built loudspeaker cabinets. Another application could be to minimize low-frequency resonances due to the room acoustics, in different places of the listening room. Yet another aim could be to obtain tonal balance and good staging.
The problem of removing undesired distortions introduced by the electro-acoustical signal path of a sound generating system is commonly called equalization, and also sometimes called dereverberation. An aim could be that the reproduced sound y(t) at a particular listening position should exactly equal the original sound w(t), but we allow it to be delayed by d samples to improve the attainable result. It is then desired that y(t) = w(t-d). Equalization by the use of digital filters has been extensively studied for about two decades, with an increasing concern in recent years for the problem of spatial robustness: A behavior close to the desired should be attained not only at one single measuring point in space, but within an extended spatial volume. Of particular importance for the present work are the time-domain properties of the impulse response of the compensated system: Differences of the dynamic responses at different listening positions may result in an adequate result for some listening positions, while the response deviates at other positions. In particular, significant sound energy may arrive before the intended delay d. Such "pre-ringings" or " pre-echoes" are considered very undesirable if their amplitudes are too large. Parts of the impulse response that are later than the target delay d may also be affected differently at different listening positions. Such " post-ringings" may significantly color the perceived spectrum and tonal balance of the sound.
In the literature, the work on robustness of equalization essentially falls into three categories.
In the first category, the goal of the filter design is a complete signal dereverberation at a single position in a room. A subsequent robustness analysis then investigates the equalizer performance at other spatial positions, or under slightly modified acoustic circumstances. It is well known that this kind of filter design is highly non-robust and causes severe signal degradation when the receiver position changes [1], and even for fixed receiver position, due to the " weak nonstationarity" of the acoustical paths in the room [2].
In the second category, the design objective is not a complete dereverberation, but rather a reduction of linear distortion under the constraint that audio performance should not be degraded by changes of listening position. The standard approach in this category is to design a filter based on averaging and/or smoothing of one or several transfer functions and then perform a robustness analysis of the filter [3]. Such methods, and in particular the complex smoothing operation proposed in [3], provide no possibilities to predict and explicitly control the amount of pre-ringing in the compensated system.
The third category imposes robustness directly on the design by employing a multi-point error criterion to optimize the sound reproduction in a number of spatial positions, either by using measured room transfer functions (RTFs) [4] or by direct adaptation of the inverse [5]. The optimization is in general based on minimum mean square error (MSE) criteria, or the sum of the power spectral densities of the compensation errors at different listening positions. MSE and power spectral density criteria do unfortunately not take the time domain properties of the compensated system into account adequately. Errors due to pre-ringings and post-ringings may result in the same MSE, although their perceptual effect can be very different. There also exists a fundamentally different multi-point scenario, where signals are filtered on the receiver side by a unique equalizer at each receiver point. Spatial robustness in this setting has been studied in [6] and [7]. This approach is however not applicable in the present pre-compensation setting, where a single filter operating on the input of a sound generating system, is designed to equalize the audio response in an extended volume in space.
Equalizers can be designed to compensate for distortions of the received energy at different frequencies. This type of filter will below be called a minimum phase inverse, or resulting in a minimum phase equalizer or magnitude equalizer. A minimum phase inverse compensates for magnitude distortions of the received signal, but does not take the phase properties (the delays of individual frequencies) of the signal into account. In the time domain, a minimum phase inverse will never create pre-ringings at any listening positions, but it may create severe post-ringings. It may even make phase and delay distortions more severe, as compared to the uncompensated system.
Both phase and magnitude distortions can be taken into account by using linear-quadratic Gaussian feedforward filter design or Wiener design, as outlined in e.g. [8]. This method has been used in [9] for designing a general class of audio precompensators. See [10]-[11] for some other FIR-filter-based methods. Design methods and resulting filters that are intended to compensate for both magnitude and phase distortions of the sound generating system will be called mixed phase methods, resulting in mixed phase equalizers. When a mixed phase equalizer is designed to compensate for non-minimum phase zeros of a transfer function, and these zeros differ in the design model and the true system, pre-ringings will unfortunately occur. The currently known mixed-phase designs provide inadequate tools for limiting the resulting pre-ringing effects.
Many researchers have concluded that mixed phase equalizers seem less robust than minimum phase equalizers from a perceptual standpoint. Their inevitable side effects in the form of pre-ringing are perceived to be more objectionable than post-ringings. Since minimum-phase equalizers create no pre-ringings, a common strategy at present for robust and perceptually acceptable equalization is therefore to use minimum phase filters only. This solution is, however, unsatisfying, as it is known to generate large phase distortions in the form of post-ringings and it cannot handle the non-minimum phase part of the audio response at all. The reference [12] proposes as a solution to limit the delay d sufficiently to make the pre-ringing inaudible. This is ineffective, since a small delay limit will for many audio systems severely restrict the ability to perform useful phase correction, in particular at the low frequencies where this is most perceptually important.
EP 1355509 A2 generally relates to digital audio precompensation and particularly the design of digital precompensation filters. Briefly, filter parameters are determined based on a weighting between, on one hand, approximating the precompensation filter to a fixed, non-zero filter component and, on the other hand, approximating the precompensated model response to a reference system response. For design purposes, the precompensation filter is preferably regarded as additively comprising a fixed, non-zero component and an adjustable compensator component. The fixed component is normally configured by the filter designer, whereas the adjustable compensator component is determined by optimizing a criterion function involving the above weighting. The weighting can be made frequency- and/or channel-dependent to provide a very powerful tool for effectively controlling the extent and amount of compensation to be performed in different frequency regions and/or in different channels.

SUMMARY OF THE INVENTION

Design techniques and convenient tools for avoiding these drawbacks are thus needed.
The present invention overcomes the difficulties encountered in the prior art.
It is a general objective of the present invention to provide an improved design scheme for audio precompensation filters.
A specific objective of the invention is to obtain a technique that can perform a mixed-phase compensation of non-minimum phase dynamics that can be safely compensated without causing any significant pre-ringings at any listening position, thereby obtaining superior compensation as compared to the known minimum phase compensators.
Another specific objective of the invention is to obtain a technique that uses similarities of different transfer functions, in the form of almost common factors of their transfer functions, or tight clusters of zeros in the complex domain, to obtain mixed phase compensators that attain a given limit on the pre-ringing effects at all or at least a subset of the listening positions.
A basic idea of the present invention is to design a discrete-time audio precompensation filter based on a Single-Input Multiple Output (SIMO) linear model ( H ) that describes the dynamic response of an associated sound generating system at a number p > 1 listening positions, for which the dynamic response differs for at least two of these listening positions. The novel filter design and construction is based on providing information representative of n non-minimum phase zeros {z_i } that are outside of the stability region |z| = 1 in the complex frequency domain, where 1 ≤ n < m, with m being the smallest number of zeros outside |z| = 1 of any of the p individual scalar models from the single input to the p outputs of the linear model H .
A characteristic of these non-minimum phase zeros is that their inversion by the precompensation filter would result in only acceptably small pre-ringings of the compensated impulse response, smaller than a prespecified limit.
The precompensation filter is then calculated as the product of at least two scalar dynamic systems, represented by:

an inverse of a characteristic scalar magnitude response in the frequency domain that represents the power gains at all or a subset of the p listening positions according to the model H ;
a causal Finite Impulse Response (FIR) filter, of user-specified degree d, having coefficients corresponding to a causal part of a delayed non-causal impulse response that is based on the information representative of n non-minimum phase zeros.

Preferably, the causal FIR filter is determined based on the information representative of n non-minimum phase zeros in the form of a design polynomial that has these n non-minimum phase zeros.
The different aspects of the invention include a method, system and computer program for designing an audio precompensation filter, a so designed precompensation filter, an audio system incorporating such a precompensation filter as well as a digital audio signal generated by such a precompensation filter.
The present invention offers the following advantages:

Optimally precompensated audio systems, resulting in superior sound quality and experience.
Provides a time-domain precision that is not attainable with minimum phase filters, while providing means to control the residual pre-ringings normally associated with mixed phase designs so that such effects can be limited to non-perceptible levels at all listening positions.

Other advantages and features offered by the present invention will be appreciated upon reading of the following description of the embodiments of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention, together with further objects and advantages thereof, will be best understood by reference to the following description taken together with the accompanying drawings, in which:

Fig. 1 is a general description of a compensated sound generating system.
Fig. 2 is a schematic flow diagram illustrating an example of the overall flow of a filter design method according to an exemplary embodiment of the invention.
Fig. 3 illustrates clusters of zeros close to the unit circle |z| =1 of the complex plane, obtained from the different transfer functions to 18 different listening positions. Zeros are represented by circles, where different radii are used to distinguish individual microphone positions. The two diagrams represent zoomed segments of the complex plane near the unit circle, at frequencies 100-150 Hz (left) and 150 - 200 Hz (right).
Fig. 4 is a schematic flow diagram illustrating an exemplary procedure for forming a design polynomial used by the present invention.
Fig. 5 shows a zoomed segment of the complex plane near the unit circle, showing the zeros of 9 room transfer functions (RTF)s marked as 'o' and their complex spatial average, marked as 'x'.
Fig. 6 is a schematic block diagram of an example of computer-based system suitable for implementation of the invention.
Fig. 7 illustrates an exemplary audio system incorporating a precompensation filter configured according to the design method of the invention.
Fig. 8 is a schematic block diagram of a single-channel setting including an equalizer filter.
Fig. 9 is a diagram illustrating a segment of the complex plane near the unit circle showing the zeros of a number of RTFs.
Fig. 10 illustrates regions in the complex plane defining the maximum tolerable zero cluster size for different zero locations.
Fig. 11 illustrates the positions of control points for design (white) and validation (black).
Fig. 12 illustrates frequency responses for filters A-F in the design points (top) and validation points (bottom).
Fig. 13 illustrates maximum level envelopes of original (grey) and equalized (black) impulse responses for filters A-F in design points (left) and validation points (right).
Fig. 14 illustrates energy step responses of original and equalized responses for filters A-C; full bandwidth (upper) and below 320 Hz (lower); design points (left) and validation points (right).
Fig. 15 illustrates average Schroeder decay sequences of original and equalized responses for filters A-C; full bandwidth (upper) and below 320 Hz (lower); design points (left) and validation points (right).
Fig. 16 illustrates energy step responses of original and equalized responses for filters D-F; full bandwidth (upper) and below 320 Hz (lower); design points (left) and validation points (right).
Fig. 17 illustrates average Schroeder decay sequences of original and equalized responses for filters D-F; full bandwidth (upper) and below 320 Hz (lower); design points (left) and validation points (right).

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

It may be useful to start with an overview of the overall flow of an exemplary filter design method according to an embodiment of the invention, with reference to the schematic flow diagram of Fig. 2. This flow diagram not only illustrates the actual design steps, but also optional pre-steps (dashed lines) that are preferably used together with the present invention, and hence represents an example of the general steps of designing a precompensation filter of the invention, starting from an uncompensated audio system and ending with an implemented filter.
We will mainly consider a sound generating system with one input signal and p different spatial listening positions, for which the dynamics response is not equal at all positions.
In step S1, a model of the audio system is provided. The model may be determined based on methods well-known for a person skilled in the art, e.g. by determining the model based on physical laws or by conducting measurements on the audio system using known test signals. Preferably the model is a Single Input Multiple Output (SIMO) model that describes the dynamic response of the associated sound generating system (i.e. the audio system) at a number p > 1 listening positions, for which the dynamic response differs for at least some of these listening positions. In step S2, there is provided information representative of n non-minimum phase zeros {z_i } that are outside of the stability region |z| = 1 in the complex frequency domain, where 1 ≤ n < m, with m being the smallest number of zeros outside |z| = 1 of any of the p individual scalar models from the single input to the p outputs of the linear SIMO model. The n non-minimum phase zeros are furthermore selected such that they have the property that their inversion by the precompensation filter results in pre-ringings of the compensated impulse response that are smaller than a prespecified limit. Step S3 involves determining a characteristic scalar magnitude response in the frequency domain that represents the power gains at all the p listening positions according to the SIMO model, and taking the inverse of this characteristic scalar magnitude response. In step S4, a causal Finite Impulse Response (FIR) filter, of user-specified degree d, and having coefficients corresponding to a causal part of a delayed non-causal impulse response is determined based on the information representative of n zeros. In step S5, the precompensation filter is determined as the product of at least the inverse of the characteristic scalar magnitude response and the causal FIR filter. In step S6, the filter parameters of the determined precompensation filter are implemented into filter hardware or software.
If required, the filter parameters may have to be adjusted. The overall design method may then be repeated, or certain steps may be repeated as schematically indicated by the dashed lines.
For example, the information of non-minimum phase zeros may be provided in the form of a design polynomial having n non-minimum phase zeros. The causal FIR filter may then be determined by:

forming a non-causal all-pass filter based on the design polynomial,
multiplying the non-causal impulse response of this non-causal all-pass filter by a time delay of d samples, to obtain a delayed non-causal impulse response, and
selecting a causal part of the delayed non-causal impulse response to obtain said FIR filter.

For a better understanding, the invention will now be described in more detail with reference to various exemplary embodiments. In the following, section 1 introduces the overall problem formulation and the polynomial notation used to describe the assumed discrete-time linear dynamic systems. Section 2 presents exemplary compensator design equations motivated by an ideal case when the systems to the different listening positions have partially common dynamics, Section 3 then generalizes the illustrative method to the case when transfer function zeros are not equal. Section 4 describes some implementation aspects. Further details of the invention, example implementations of algorithm steps and performance examples are provided in Appendix 1.

1. SYSTEM MODEL

It is useful to begin by describing the general approach for designing audio precompensation filters.
The sound generation or reproducing system to be modified is normally represented by a linear time-invariant dynamic model H that describes the relation in discrete time between a set of m input signals u(t) to a set of p output signals y(t): $\begin{array}{l} y (t) = H u (t) \\ y_{m} (t) = y (t) + e (t), \end{array}$

where t represents a discrete time index, y_m(t) (with subscript m denoting " measurement") is a p-dimensional column vector representing the sound time-series at p different locations and e(t) represents noise, unmodeled room reflexes, effects of an incorrect model structure, nonlinear distortion and other unmodeled contributions. The operator H is a p×m-matrix whose elements are stable linear dynamic operators or transforms, e.g. represented as FIR filters or IIR filters. These filters will determine the response y(t) to a m-dimensional arbitrary input time series vector u(t). Linear filters or models will be represented by such matrices, which are called transfer function matrices, or dynamic matrices, in the following. The transfer function matrix H represents the effect of the whole or a part of the sound generating or sound reproducing system, including any pre-existing digital compensators, digital-to-analog converters, analog amplifiers, loudspeakers, cables and the room acoustic response. In other words, the transfer function matrix H represents the dynamic response of relevant parts of a sound generating system. When including the dynamics of the listening room, its elements are called room transfer functions, or RTFs. The input signal u(t) to this system, which is a m-dimensional column vector, may represent input signals to m individual amplifier-loudspeaker chains of the sound generating system.
The measured sound y_m(t) is, by definition, regarded as a superposition of the term y(t) = Hu(t) that is to be modified and controlled, and the unmodeled contribution e(t). A prerequisite for a good result in practice is, of course, that the modeling and system design is such that the magnitude |e(t)| will not be large compared to the magnitude |y(t)|, in the frequency regions of interest.
A general objective is to modify the dynamics of the sound generating system represented by (1.1) in relation to some reference dynamics. For this purpose, a reference matrix D is introduced: $y_{ref} (t) = D w (t),$

where w(t) is an r-dimensional vector representing a set of live or recorded sound sources or even artificially generated digital audio signals, including test signals used for designing the filter. The elements of the vector w(t) may, for example, represent channels of digitally recorded sound, or analog sources that have been sampled and digitized. In (1.2), D is a transfer function matrix of dimension m×r that is assumed to be known. The linear system D is a design variable and generally represents the reference dynamics of the vector y(t) in (1.1).
An example of a conceivable design objective may be complete inversion of the dynamics and decoupling of the channels. In cases where r = p, the matrix D is then set equal to a square diagonal matrix with d-step delay operators as diagonal elements, so that: $y_{ref} (t) = w (t - d) .$
The reference response of y(t) is then defined as being just a delayed version of the original sound vector w(t), with equal delays of d sampling periods for all elements of w(t). The desired bulk delay, d, introduced through the design matrix D is an important parameter that influences the attainable performance. Causal mixed-phase compensation filters will attain better compensation the higher this delay is allowed to be.
The precompensation is generally obtained by a precompensation filter, generally denoted by R, which generates an input signal vector u(t) to the audio reproduction system (1.1) based on the signal w(t): $u (t) = R w (t) .$
A commonly used design objective is then to generate the input signal vector u(t) to the audio reproduction system (1.1) so that its compensated output y(t) approximates the reference vector y_ref(t) well, in some specified sense. This objective can be attained if the signal u(t) in (1.1) is generated by a linear precompensation filter R , which consists of a m×r-matrix whose elements are stable and causal linear dynamic filters that operate on the signal w(t) such that y(t) will approximate y_ref(t): $\begin{matrix} y (t) = H u (t) \end{matrix} = HR w (t) ≅ y_{ref} (t) = D w (t) .$
Linear discrete-time dynamic systems are in the following represented using the discrete-time backward shift operator, here denoted by q^-1 . A signal s(t) is shifted backward by one sample by this operator: q^-1s(t) = s(t-1). Likewise, the forward shift operator is denoted q, so that qs(t) = s(t+1), using the notation of e.g. [15]. A single-input single output (scalar) design model (1.1) is then represented by a linear time-invariant difference equation with fixed scalar coefficients: $y (t) = - a_{1} y (t - 1) - a_{2} y (t - 2) - \dots - a_{n} y (t - n) + b_{0} u (t - k) + b_{1} u (t - k - 1) + \dots + b_{b} u (t - k - h) .$
Assuming b₀ ≠ 0, there will be a delay of k samples before the input u(t) influences the output y(t). This delay, k, may for example represent an acoustic transport delay and it is here called the bulk delay of the model. The coefficients a_j and b_j determine the dynamic response described by the model. The maximal delays n and h may be many hundreds or even thousands of samples in some models of audio systems.
Move all terms related to y to the left-hand side. With the shift operator representation, the model (1.5) is then equivalent to the expression: $(1 + a_{1} q^{- 1} + a_{2} q^{- 2} + \dots + a_{n} q^{- n}) y (t) = (b_{0} + b_{1} q^{- 1} + \dots + b_{h} q^{- h}) u (t - k) .$
By introducing the polynomials A(q^-1) = 1 + a₁q^-1 + a₂q^-2 +... + a_nq^-n and B(q^-1) = b₀ + b₁q^-1 + ... + b_hq^-h , the discrete-time dynamic model (1.5) may be represented by the more compact expression: $A (q^{- 1}) y (t) = B (q^{- 1}) u (t - k) .$
The polynomial A(q^-1) is said to be monic since its leading coefficient is 1. In the special case of Finite Impulse Response (FIR) models, A(q^-1) = 1. In general, the recursion in old outputs y(t-j) represented by the filter A(q^-1) gives the model an infinite impulse response. Infinite Impulse Response (IIR) filters represented in the form (1.6) are also denoted rational filters, since their transfer operator may be represented by a ratio of polynomials in q^-1 : $y (t) = \frac{B (q^{- 1})}{A (q^{- 1})} u (t - k) .$
All involved IIR systems, models and filter are in the following assumed to be stable. The stable criterion means that, when a complex variable z is substituted for the operator q, this is equivalent to the equation A(z^-1) = 0 having solutions with magnitude |z| < 1 only. In other words, the complex function A (z^-1) must have all zeros within the unit circle in the complex plane. The roots of A(z^-1) = 0 are the poles of the scalar system. The roots of B(z^-1) = 0 are its zeros. The roots of B(z^-1) = 0 with magnitude |z| < 1 are the minimum phase zeros. Their steady-state influence can be eliminated by a stable compensator system with poles at the corresponding position. The roots of B(z^-1) = 0 with magnitude |z| ≥ 1 are the non-minimum phase zeros. Their influence cannot be eliminated completely by a stable and causal compensator. Their influence can, however, be approximately eliminated by using a stable but non-causal compensator that operates on future time samples. Such compensators can be realized by delaying all involved signals in an appropriate way. In particular, the delay d in (1.3) can be used for this purpose. In the following, we will also assume that all bulk delays are absorbed by the delay d in D, so that k=0 will be used in all transfer functions that are elements of H.
For any polynomial B(q^-1) = b_o + b_Lq^-1 + ... + b_hq^-h , we further define the conjugate polynomial B*(q) = b_o + b_lq + ... + b_h ^h and the reciprocal polynomial B̅ (q^-1 ) = b_n + b_n-1q^-1 + ... + b_oq^-h . The zeros of B*(z) are reflected in the unit circle with respect to the zeros of B(z^-1).

2. DESIGN FOR COMMON MODEL FACTORS

In the following, we will focus on sound reproduction systems with scalar inputs m=1 (i.e. single loudspeaker and amplifier chains), that has multiple listening positions r=p>1. The compensator (1.4) is thus a scalar discrete-time dynamic system. The listening positions are also called control points. The desired response is assumed to be equal to a delay of d samples for all control points, as described by (1.3).
In this example, we use a discrete-time linear single input multiple output (SIMO) model structure y(t) = H(q^-1)u(t) in common denominator form: $y_{i} (t) = \frac{B_{i} (q^{- 1})}{A (q^{- 1})} u (t), i = 1, \dots, p .$
The stable denominator A(q^-1) is thus the same in the models to all p outputs. This property can always be obtained by writing transfer functions on common denominator form. We will furthermore in this section assume that a set of zeros, of which n>1 are non-minimum phase zeros, are common to all transfer functions. The common zeros define the robust part B^r (q ^-1) of the numerator polynomials, while the other zeros comprise the non-robust parts, which are not assumed equal for all outputs: $B_{i} (q^{- 1}) = B^{r} (q^{- 1}) B_{i}^{n} (q^{- 1}), i = 1, \dots, p .$
The robust numerator factor $B_{u}^{r} (q^{- 1})$
furthermore is the product of a non-minimum phase (" unstable") factor $B_{u}^{r} (q^{- 1})$
with n>1 zeros outside the unit circle and a minimum phase (" stable") factor $B_{s}^{r} (q^{- 1}) : B^{r} (q^{- 1}) = B_{u}^{r} (q^{- 1}) B_{s}^{r} (q^{- 1}) .$
The assumed design model can thus be decomposed into a common scalar IIR model in series with a set of FIR models: $\begin{array}{l} z (t) = \frac{B^{r} (q^{- 1})}{A (q^{- 1})} u (t) = \frac{B_{u}^{r} (q^{- 1}) B_{s}^{r} (q^{- 1})}{A (q^{- 1})} u (t) \\ y_{i} (t) = B_{i}^{n} (q^{- 1}) z (t) . \end{array}$
The scalar intermediate signal z(t) has been introduced in the above model structure for pedagogical purposes. It is not assumed to be a physically existing measurable signal.
We furthermore define two scalar design models:

The complex spatial average model, defined as: $\frac{B_{0} (q^{- 1})}{A (q^{- 1})} = \sum_{i = 1}^{p} \frac{B_{i} (q^{- 1})}{A (q^{- 1})} = \frac{B^{r} (q^{- 1})}{A (q^{- 1})} \sum_{i = 1}^{p} B_{i}^{n} (q^{- 1}) .$

The complex gain of B₀(z^-1)/A(z^-1 ) represents the average of the complex gains of the individual models, pointwise in the frequency domain. Note that the robust numerator polynomial B^r(q^-1) of the individual models will always be a factor of B₀(q^-1).
The Root Mean Square (RMS) average model, also known as the minimum phase equivalent model: $\frac{β (q^{- 1})}{A (q^{- 1})},$

where the stable numerator polynomial β(q ^-1) is obtained from the regularized spectral factorization equation: $β_{*} (q) β (q^{- 1}) = ρ + \sum_{i = 1}^{p} B_{i *} (q) B_{i} (q^{- 1}) = ρ + B_{*}^{r} (q) B^{r} (q^{- 1}) \sum_{i = 1}^{p} B_{i *}^{n} (q) B_{i}^{n} (q^{- 1})$

The RMS average model thus sums the power spectra of the p design models, and then calculates a stable minimum phase system that has the corresponding magnitude spectrum. This model thus includes no phase information of the individual models. The optional regularization term ρ corresponds to the use of an energy penalty on the filter output u(t) in an LQG feedforward compensator filter design, see [8]. It can thus be used for limiting the filter output energy. It also limits the depth of notches in the magnitude response of (2.5).
A regularization parameter with frequency-dependent magnitude may be used. This can be done by generalizing the spectral factorization equation (2.6) to $β_{*} (q) β (q^{- 1}) = ρ W_{*} (q) W (q^{- 1}) + \sum_{i = 1}^{p} B_{i *} (q) B_{i} (q^{- 1}) = ρ W_{*} (q) W (q^{- 1}) + B_{*}^{r} (q) B^{r} (q^{- 1}) \sum_{i = 1}^{p} B_{i *}^{n} (q) B_{i}^{n} (q^{- 1}) .$
Here, W(q^-1) is a user-defined FIR weighting filter that provides frequency weighting of the penalty/regularization term. When ρW*(z)W(z^-1) >0 on the unit circle |z|=1, then a stable spectral factor β(q ^-1), which has no zeros on or outside the unit circle, can always be found. With no penalty/regularization ρ=0, a stable spectral factor exists if the polynomials B_i(q^-1) do not all have a common factor with zeros on the unit circle.
Finally, based on the unstable part of the robust numerator, of order n, $B_{u}^{r} (q^{- 1}) = f_{0} + f_{1} q^{- 1} + \dots + f_{n} q^{- n}$

define a design polynomial F(q) by $F (q) = {\overline{B}}_{u *}^{r} (q) = q^{n} B_{u}^{r} (q^{- 1}) = f_{n} + f_{n - 1} q + \dots + f_{0} q^{n}$

and its reciprocal polynomial $\overline{F} (q) = B_{u *}^{r} (q) = f_{0} + f_{1} q + \dots + f_{n} q^{n} .$
The proposed exemplary design of a robust scalar compensator filter for the single-input p output model (1.9) is then given by $u (t) = Rw (t) = Q (q^{- 1}) \frac{A (q^{- 1})}{β (q^{- 1})} w (t)$

where the last factor is the inverse of the RMS average model (2.5) while the polynomial (FIR filter) Q(q^-1) will have degree d equal to the specified design delay. It performs phase compensation, while the inverse of the RMS average performs magnitude compensation only.
In general, a filter Q(q^-1) that minimizes a prescribed quadratic design criterion can be obtained by solving a linear polynomial equation, a Diophantine equation [8],[15].
We in the following assume that the desired response is w(t-d) at all p measuring positions, i.e. that the matrix D is a column vector containing terms q^-d . The bulk propagation delays are assumed to be taken care of separately, so that this formulation is sensible also for large listening volumes. A Linear Quadratic Gaussian design for the system (2.3), aiming at minimizing the quadratic criterion $E {‖ y (t) ‖}^{2} + E ρ {‖ W (q^{- 1}) u (t) ‖}^{2}$

where E(.) represents an average over a white input sequence w(t). The compensator is then given by (2.11), where β(q ^-1) is obtained from the spectral factorization equation (2.7) and Q(q^-1), together with a polynomial L*(q) is obtained as the unique solution to the following scalar Diophantine polynomial equation $q^{- d} B_{0 *} (q) = β_{*} (q) Q (q^{- 1}) + q L_{*} (q),$

see equation (8) in Appendix 1 [13]. For solutions to more general multi-input LQG feedforward filter design problem formulations, see the equation (3.3.11) in [8].
When d is large and the models for the p listening positions are significantly different, a prefilter designed according to (2.11), (2.7) and (2.13) will in general result in compensated impulse responses with significant pre-ringings at least at some listening positions. A single compensator simply cannot compensate the different dynamics represented by $B_{i}^{n} (q^{- 1})$
in (2.3) exactly. In this situation, a useful way forward is to base the design of the RMS inverse part in (2.11) and of the polynomial Q(q^-1) on different system models and criteria. We here propose the following hybrid design:

1. The factor A(q^-1)/β(q^-1) is kept unchanged in (2.11), with β(q ^-1) being obtained from the spectral factorization (2.7). This can be interpreted as a minimum phase equalizer, obtained by minimizing the criterion (2.12) with no allowed delay d=0.
2. The FIR filter Q(q^-1) is obtained by using a modified criterion and concentrating on compensating the robust part of the model (2.3) only.

If the second step is aimed at minimizing the criterion $E {‖ z (t) ‖}^{2} + E ρ {‖ W_{1} (q^{- 1}) u (t) ‖}^{2},$

where z(t) is the (non-measurable) intermediate signal in (2.3) and the frequency weighting factor W₁(q^-1) may differ from the factor W(q^-1) used in (2.12) for the first step, then design of Q(q^-1) involves solving the polynomial spectral factorization $B_{*}^{r} (q) β^{r} (q^{- 1}) = ρ W_{1 *} (q) W_{1} (q^{- 1}) + B_{*}^{r} (q) B^{r} (q^{- 1})$

with respect to the stable " robust spectral factor" polynomial β ^r (q ^-1). The filter polynomial Q(q^-1) is then obtained by solving the modified Diophantine equation $q^{- d} \sqrt{p} B_{*}^{r} (q) = B_{*}^{r} (q) Q (q^{- 1}) + q L_{*} (q) .$
Here, Q(q^-1) has degree d, and L*(q) has degree one less than $B_{*}^{r} (q) .$
By dividing (2.16) with $β_{*}^{r} (q),$
it is evident that the polynomial Q(q^-1) represents the causal part (the part for negative time shifts) of the impulse response of the filter $q^{- d} \sqrt{p} \frac{B_{*}^{r} (q)}{B_{*}^{r} (q)} = Q (q^{- 1}) + q \frac{L_{*} (q)}{B_{*}^{r} (q)} = m_{- d} q^{- d} + m_{- d + 1} q^{- d + 1} + \dots + m_{0} + m_{1} q + m_{2} q^{2} + \dots$

so $Q (q^{- 1}) = m_{0} + m_{1} q + \dots + m_{- d} q^{- d} .$
When using no input penalty in (2.14), the spectral factorization (2.15) simplifies and the relevant expression reduces to an expression of lower degree: $q^{- d} \sqrt{p} \frac{B_{*}^{r} (q)}{β_{*}^{r} (q)} = q^{- d} \sqrt{p} \frac{B_{s *}^{r} (q) B_{u *}^{r} (q)}{B_{s *}^{r} (q) {\overline{B}}_{u *}^{r} (q)} = q^{- d} \sqrt{p} \frac{\overline{F} (q)}{F (q)} .$
The first equality follows from the properties of the spectral factorization and it shows that the factor $B_{s *}^{r} (q)$
can be cancelled when ρ=0. Use of the definitions (2.9) and (2.10) then give the second equality. This impulse response can in this case be expressed by $\sqrt{p} q^{- d} (f) \frac{(_{0} + f_{1} q + \dots + f_{n} q^{n})}{(f_{n} + f_{n - 1} q + \dots + f_{0} q^{n})} = m_{- d} q^{- d} + m_{- d + 1} q^{- d + 1} + \dots + m_{0} + m_{1} q + m_{2} q^{2} + ....$
The purpose of the FIR filter (2.18) is to approximately invert the non-minimum phase dynamics that is represented by the unstable part (2.8) of the robust numerator. The fidelity of this approximation is improved by increasing the delay parameter d, which allows a larger part of the total impulse response (2.17) to be used by the compensator filter (2.18). As d →∞, perfect inversion is approached.
Let us consider three special cases of the proposed compensation filter:

1. When d=0, then Q(q^-1) just equals the scale factor $\sqrt{p} .$
The design reduces to the inverse of the RMS average model.¹ This is the minimum mean square optimal minimum phase equalizer. It equalizes the RMS average model and thus achieves a magnitude response compensation of the average response at the p listening positions. It does not compensate the phase properties of the audio system.
¹ The scale factor $\sqrt{p}$
is needed since the RMS average model as defined here has not been normalized to keep the gain constant as the number of models grows. An alternative formulation where such normalization is used is a trivial modification of the scheme.
2. When p=1, or when all p models are exactly equal, then Bⁿ (q ^-1) = 1 and by (2.4), B₀ (q ^-1) = pB^r (q ^-1). The relevant model then exactly equals the (scalar) robust model. The resulting filter (2.11) is then the MMSE mixed phase filter or LQ optimal feedforward compensator [8] for this model. The solution to this special case is previously known.
3. The case of p > 1 listening position where all p models (2.1) are not equal is the situation of primary interest. For this case, the compensator filter defined by (2.7), (2.11),(2.15) and (2.16), represents a novel solution to the robust digital audio precompensation problem. It performs power compensation of the whole set of models since it contains the inverse of the RMS average model as factor. The factor Q(q^-1) in addition performs robust phase compensation. It represents a causal realizable Wiener design that approximately inverts only the non-minimum phase zeros that are present in the models to all p measurement positions. Pre-ringings of the compensated impulse response coefficients (before coefficient d), due to erroneous compensation will therefore not occur. This design thus combines total power equalization with robust phase compensation for the assumed model structure.

The Appendix 1 [13] provides a further discussion of the properties of these three variants of the filters and the corresponding compensated systems.
While it is useful to design the compensator (2.11) in the form of two separate filters that are connected in series, the actual implementation of the compensator may differ, and be performed in the form most convenient for the application. For example, the inverse of the RMS average factor may be approximated by a FIR filter so that the whole compensator becomes a FIR filter. Another alternative is to implement parts of the compensator as a series connection of second order links (a so-called biquad filter).
The design has here been outlined for a sound generating system with scalar input signal. The design can be obviously generalized to sound generating systems with an arbitrary number m of inputs, by performing the scalar design separately for each input, using an appropriate number of listening positions.
Another useful modification is to use the proposed technique in the modified LQG design technique that is used in conjunction with a specific filter structure in [9]. There, the precompensation filter is regarded as additively decomposed into a fixed filter component in parallel with an adjustable component. The use of this filter structure improves the design utility of frequency weighted input penalty terms in the quadratic criterion.

3. DESIGN FOR APPROXIMATELY COMMON MODEL FACTORS

The design of the previous section was based on the assumption that the robust part of the model (2.3) was exactly the same in all the p models. This situation will unfortunately rarely occur in practice. We now generalize the design to the more relevant case where the models contain almost common factors (or near-common factors).
Basing a robust design on average dynamics, or on properties of the p-output system that are relatively stable in the spectral properties at all the measurement positions, has been discussed as a general principle in the prior art [18]. We have found such an approach to be too simplistic: Just because a property of a spectrum looks rather similar when comparing the different models, this does not imply that it will always be safe to use it in a mixed-phase equalizer: Inversion of a dynamic system is a nonlinear operation, with sometimes surprising results. Instead of focusing on the similarity of properties of different models, the present invention bases the design on a different principle: It is based on performing an explicit pre-analysis of the consequences that would occur, if a specific property of the p models, or of an aggregate model, were used for the compensator design. The causal FIR filter forming part of the overall filter design is preferably formed to approximately invert only non-minimum phase zeros that, by the preceding analysis, can be safely inverted.
For design purposes, we here assume a scalar input signal u(t), p output signals y_i(t) and a model (1.7) in common denominator form. Exact common numerator factors are no longer assumed. The zeros of any numerator B_i (q ^-1) may thus differ from all the zeros in all the p-1 other model numerator polynomials B_j (q ^-1), j ≠ i.²
2 Small deviations of poles can be transformed to approximately equivalent uncertainties/deviations of numerator polynomials, or zeros, of a transfer function, by using the robust design methodology of [17].
The complex spatial average model is, as in section 2, given by: $\frac{B_{0} (q^{- 1})}{A (q^{- 1})} = \sum_{i = 1}^{p} \frac{B_{i} (q^{- 1})}{A (q^{- 1})} .$
As candidate properties for use in the robust compensation of non-minimum phase zeros, we will focus on the non-minimum phase zeros of the complex spatial average model (3.1). In section 2, we saw that the robust part of the numerators, if it exists, would always be a factor of the numerator of the complex spatial average model. We will now denote the polynomial formed by these sufficiently common non-minimum phase zeros by F. Note however that the polynomial F will no longer be exactly related to an unstable part of a robust numerator (2.8), as it was designed in (2.9).
In the RMS average model (2.5), the numerator polynomial β(q ^-1) is, as before, obtained from the regularized spectral factorization equation: $β_{*} (q) β (q^{- 1}) = ρ W_{*} (q) W (q^{- 1}) + \sum_{i = 1}^{p} B_{i *} (q) B_{i} (q^{- 1}) .$
As an example of a situation of interest consider Figure 3. It illustrates clusters of zeros close to the unit circle |z| = 1 of the complex plane, obtained from the different transfer functions to 18 different listening positions [16]. Zeros are represented by circles, where different radii are used to distinguish individual microphone positions. The two diagrams represent zoomed segments of the complex plane near the unit circle, at frequencies 100-150 Hz (left) and 150 - 200 Hz (right). As is evident from Figure 3, zeros of the room transfer functions move around as the microphone position changes. In particular, slightly above 200 Hz, a zero moves from the inside of the unit circle to the outside - a typical example of a zero that cannot be inverted without causing severe pre- or post-ringing errors in most listening positions. However, some zeros further out from the unit circle exhibit a more static behavior. Based on these observations, it seems reasonable to assume that some non-minimum phase zeros could be safely inverted under a constraint of maximum tolerable residual pre-ringings.
A proposed exemplary design starts from the zeros of the complex spatial average model (3.1). Reference can be made to the schematic flow diagram of Fig. 4.

1. For the p design models, clusters of non-minimum phase zeros are formed (S11), with at most one zero of a cluster belonging to any model polynomial B_i(z^-1).
2. A zero of the complex average model (3.1) is associated (S12) with each of these clusters.
3. For each non-minimum phase zero of B₀(q^-1 ) in the complex average model, the design algorithm evaluates (EVALUATION in Fig. 4) the pre-ringing that would be obtained if a compensator that uses a given design delay d would include this particular dynamics into the design polynomial F(q), and then approximately invert it according to (2.7)-(2.18). The following procedure is performed for this purpose:
1. a. The resulting pre-ringings are estimated (S13). This estimate may be based on the positions and size of the cluster of model zeros that has been associated with this particular non-minimum phase zero of B₀(q^-1 ).
2. b. The so estimated pre-ringings are compared (S14) against a limit on the impulse response coefficients obtained from a limit on acceptable pre-ringings.
3. c. If these impulse response coefficients, as evaluated based on the p design models, all have magnitudes below this limit, a second order polynomial factor defined by the evaluated zero and its complex conjugate is included (S15) into the polynomial $B_{u}^{r} (q^{- 1}) = f_{0} + f_{1} q^{- 1} + \dots + f_{n} q^{- n}$
  introduced in (2.8).
4. d. If the modeled pre-ringing is too large, the zero is not included in (2.8).
4. After having evaluated all non-minimum phase zeros of B₀(q^-1 ) in the complex average model (3.1), and thus having formed the complete design polynomial F(q), a robust precompensator (i.e. the audio filter) may be designed as specified by (2.7), (2.11), (2.15) and (2.16).

A non-causal exponential decay function, i.e. a function that decays exponentially forward in time, is a specific example of a time domain criterion for the impulse response coefficients that can be used as a limit for the pre-ringings: $20 \log_{10} ({Cr}_{0}^{- κ}) < L_{\min} .$
Here, C and r₀ are scaling constants, κ > 0 is a time constant and L _min is the maximum tolerable pre-ringing level, measured in dB, in the equalized response h(k) at time index k = d - κ. This constraint ensures that the pre-ringing level generated by compensating the evaluated zero is at most L _min dB at all time instants prior to k = d - κ. It should be noted that many zeros are likely to contribute to the total pre-ringing. The individual constraint (3.3) should be adjusted with this in mind.
It has above been assumed that a cluster of zeros can be associated with the zero of B₀(q^-1 ) under evaluation. The limit (3.3) is then evaluated with respect to the positions of zeros belonging to this cluster, and their relations to the zero under evaluation. Figure 5 illustrates such a situation. The right-hand cross is the non-minimum phase zero (outside the unit circle) of the complex average model under evaluation. The rings around it are the cluster of zeros of nine individual room transfer functions. Sect. V.C of Appendix 1 [13], presents equations (38)-(40) that approximately quantify the amount of pre-ringings that would result: For a conjugate pair of nominal zeros zo and conj(z₀ ) parameterized as $z_{0} = r_{0} \exp (jω) (_{0}), conj (z) (_{0}) = r_{0} \exp (- {jω}_{0}), where r_{0} = |z_{0}|,$

they, together with the pre-ringing envelope constraint (3.3) and its parameters C, L _min and κ, implicitly define a region around z₀ within which a zero cluster must be contained in order for z₀ to be considered a common, robustly invertible, zero of all the involved room transfer functions. It can from Fig. 3 in Appendix 1 be observed that zero clusters can be allowed to be larger if the nominal zero z₀ is located further away from the unit circle, than when z₀ is close to the unit circle.
The above exemplary algorithm preferably works in conjunction with a scheme for clustering of near-common excess phase zeros that belong to different room transfer function models. One particular such scheme is outlined in section V.F of Appendix 1.
We thus obtain a method for trading phase compensation fidelity against bounds on the amount of pre-ringing at any listening position. For a given set of p models and given design delay d, the technique of this invention provides a tool for optimizing frequency domain design criteria while satisfying a time-domain pre-ringing constraint.
Finally, let us point out that the RMS spatial average model defined by (2.5) and (3.2) may optionally be processed by smoothing in the frequency domain before its inverse is used in the compensator filter (2.11). The reason for this is that the number p of transfer function measurements is limited. The number is in most practical filter designs much smaller than that needed to interpolate the whole listening space at the frequencies of interest. Therefore, the RMS spatial average model will not represent the true RMS average over all possible listening positions. In particular, the sample RMS average tends to have a more irregular frequency response than the true RMS average. A method can therefore be used to better estimate the true RMS average based on a limited number of measured room transfer functions. One possible and commonly used approach, which produces a less irregular frequency response, is to use a smoothing in the frequency domain of the finite sample RMS average (2.4). It may be advantageous to use a frequency-dependent degree of smoothing, as provided by the well-known octave smoothing filter, and as is also discussed in [14].³
³ The reference [14] proposes the principle of using smoothing that is variable with frequency across the signal spectrum as a patented claim. This is however a design technique that is previously well-known in the form of e.g. octave smoothing.

4. IMPLEMENTATIONAL ASPECTS

Typically, the design equations are solved on a separate computer system to produce the filter parameters of the precompensation filter. The calculated filter parameters are then normally downloaded to a digital filter, for example realized by a digital signal processing system or similar computer system, which executes the actual filtering.
Although the invention can be implemented in software, hardware, firmware or any combination thereof, the filter design scheme proposed by the invention is preferably implemented as software in the form of program modules, functions or equivalent. The software may be written in any type of computer language, such as C, C + + or even specialized languages for digital signal processors (DSPs). In practice, the relevant steps, functions and actions of the invention are mapped into a computer program, which when being executed by the computer system effectuates the calculations associated with the design of the precompensation filter. In the case of a PC-based system, the computer program used for the design of the audio precompensation filter is normally encoded on a computer-readable medium such as a DVD, CD or similar structure for distribution to the user/filter designer, who then may load the program into his/her computer system for subsequent execution. The software may even be downloaded from a remote server via the Internet.
Fig. 6 is a schematic block diagram illustrating an example of a computer system suitable for implementation of a filter design algorithm according to the invention. The system 100 may be realized in the form of any conventional computer system, including personal computers (PCs), mainframe computers, multiprocessor systems, network PCs, digital signal processors (DSPs), and the like. Anyway, the system 100 basically comprises a central processing unit (CPU) or digital signal processor (DSP) core 10, a system memory 20 and a system bus 30 that interconnects the various system components. The system memory 20 typically includes a read only memory (ROM) 22 and a random access memory (RAM) 24. Furthermore, the system 100 normally comprises one or more driver-controlled peripheral memory devices 40, such as hard disks, magnetic disks, optical disks, floppy disks, digital video disks or memory cards, providing non-volatile storage of data and program information. Each peripheral memory device 40 is normally associated with a memory drive for controlling the memory device as well as a drive interface (not illustrated) for connecting the memory device 40 to the system bus 30. A filter design program implementing a design algorithm according to the invention, possibly together with other relevant program modules, may be stored in the peripheral memory 40 and loaded into the RAM 22 of the system memory 20 for execution by the CPU 10. Given the relevant input data, such as a model representation and other optional configurations, the filter design program calculates the filter parameters of the precompensation filter.
The determined filter parameters are then normally transferred from the RAM 24 in the system memory 20 via an I/O interface 70 of the system 100 to a precompensation filter system 200. Preferably, the precompensation filter system 200 is based on a digital signal processor (DSP) or similar central processing unit (CPU) 202, and one or more memory modules 204 for holding the filter parameters and the required delayed signal samples. The memory 204 normally also includes a filtering program, which when executed by the processor 202, performs the actual filtering based on the filter parameters.
Instead of transferring the calculated filter parameters directly to a precompensation filter system 200 via the I/O system 70, the filter parameters may be stored on a peripheral memory card or memory disk 40 for later distribution to a precompensation filter system, which may or may not be remotely located from the filter design system 100. The calculated filter parameters may also be downloaded from a remote location, e.g. via the Internet, and then preferably in encrypted form.
In order to enable measurements of sound produced by the audio equipment under consideration, any conventional microphone unit(s) or similar recording equipment 80 may be connected to the computer system 100, typically via an analog-to-digital (A/D) converter 80. Based on measurements of (conventional) audio test signals made by the microphone 80 unit, the system 100 can develop a model of the audio system, using an application program loaded into the system memory 20. The measurements may also be used to evaluate the performance of the combined system of precompensation filter and audio equipment. If the designer is not satisfied with the resulting design, he may initiate a new optimization of the precompensation filter based on a modified set of design parameters.
Furthermore, the system 100 typically has a user interface 50 for allowing user-interaction with the filter designer. Several different user-interaction scenarios are possible.
For example, the filter designer may decide that he/she wants to use a specific, customized set of design parameters in the calculation of the filter parameters of the filter system 200. The filter designer then defines the relevant design parameters via the user interface 50.
It is also possible for the filter designer to select between a set of different preconfigured parameters, which may have been designed for different audio systems, listening environments and/or for the purpose of introducing special characteristics into the resulting sound. In such a case, the preconfigured options are normally stored in the peripheral memory 40 and loaded into the system memory during execution of the filter design program.
The filter designer may also define the reference system by using the user interface 50. In particular, the delay d of the reference system may be selected by the user, or provided as a default delay.
Instead of determining a system model based on microphone measurements, it is also possible for the filter designer to select a model of the audio system from a set of different preconfigured system models. Preferably, such a selection is based on the particular audio equipment with which the resulting precompensation filter is to be used.
Preferably, the audio filter is embodied together with the sound generating system so as to enable generation of sound influenced by the filter.
In an alternative implementation, the filter design is performed more or less autonomously with no or only marginal user participation. An example of such a construction will now be described. The exemplary system comprises a supervisory program, system identification software and filter design software. Preferably, the supervisory program first generates test signals and measures the resulting acoustic response of the audio system. Based on the test signals and the obtained measurements, the system identification software determines a model of the audio system. The supervisory program then gathers and/or generates the required design parameters and forwards these design parameters to the filter design program, which calculates the precompensation filter parameters. The supervisory program may then, as an option, evaluate the performance of the resulting design on the measured signal and, if necessary, order the filter design program to determine a new set of filter parameters based on a modified set of design parameters. This procedure may be repeated until a satisfactory result is obtained. Then, the final set of filter parameters are downloaded/implemented into the precompensation filter system.
It is also possible to adjust the filter parameters of the precompensation filter adaptively, instead of using a fixed set of filter parameters. During the use of the filter in an audio system, the audio conditions may change. For example, the position of the loudspeakers and/or objects such as furniture in the listening environment may change, which in turn may affect the room acoustics, and/or some equipment in the audio system may be exchanged by some other equipment leading to different characteristics of the overall audio system. In such a case, continuous or intermittent measurements of the sound from the audio system in one or several positions in the listening environment may be performed by one or more microphone units or similar sound recording equipment. The recorded sound data may then be fed into a filter design system, such as system 100 of Fig. 6, which calculates a new audio system model and adjusts the filter parameters so that they are better adapted for the new audio conditions.
Naturally, the invention is not limited to the arrangement of Fig. 6. As an alternative, the design of the precompensation filter and the actual implementation of the filter may both be performed in one and the same computer system 100 or 200. This generally means that the filter design program and the filtering program are implemented and executed on the same DSP or microprocessor system.
A sound generating or reproducing system 300 incorporating a precompensation filter system 200 according to the present invention is schematically illustrated in Fig. 7. An audio signal w(t) from a sound source is forwarded to a precompensation filter system 200, possibly via a conventional I/O interface 210. If the audio signal w(t) is analog, such as for LPs, analog audio cassette tapes and other analog sound sources, the signal is first digitized in an A/D converter 210 before entering the filter 200. Digital audio signals from e.g. CDs, DAT tapes, DVDs, mini discs, and so forth may be forwarded directly to the filter 200 without any conversion.
The digital or digitized input signal w(t) is then precompensated by the precompensation filter 200, basically to take the effects of the subsequent audio system equipment into account.
The resulting compensated signal u(t) is then forwarded, possibly through a further I/O unit 230, for example via a wireless link, to a D/A-converter 240, in which the digital compensated signal u(t) is converted to a corresponding analog signal. This analog signal then enters an amplifier 250 and a loudspeaker 260. The sound signal y_m(t) emanating from the loudspeaker 260 then has the desired audio characteristics, giving a close to ideal sound experience. This means that any unwanted effects of the audio system equipment have been eliminated through the inverting action of the precompensation filter.
The precompensation filter system may be realized as a stand-alone equipment in a digital signal processor or computer that has an analog or digital interface to the subsequent amplifiers, as mentioned above. Alternatively, it may be integrated into the construction of a digital preamplifier, a computer sound card, a compact stereo system, a home cinema system, a computer game console, a TV, an MP3 player docking station or any other device or system aimed at producing sound. It is also possible to realize the precompensation filter in a more hardware-oriented manner, with customized computational hardware structures, such as FPGAs or ASICs.
It should be understood that the precompensation may be performed separate from the distribution of the sound signal to the actual place of reproduction. The precompensation signal generated by the precompensation filter does not necessarily have to be distributed immediately to and in direct connection with the sound generating system, but may be recorded on a separate medium for later distribution to the sound generating system. The compensation signal u(t) in Fig, 1 could then represent for example recorded music on a CD or DVD disk that has been adjusted to a particular audio equipment and listening environment. It can also be a precompensated audio file stored on an Internet server for allowing subsequent downloading of the file to a remote location over the Internet.
The embodiments described above are merely given as examples, and it should be understood that the present invention is not limited thereto. Further modifications, changes and improvements that retain the basic underlying principles disclosed and claimed herein are within the scope of the invention.

REFERENCES

[1] B. D. Radlovic, R. C. Williamson and R. A. Kennedy, " Equalization in an acoustic reveberant environment: Robustness results," IEEE Transactions on Speech and Audio Processing, vol. 8, no. 3, pp. 311-319, May 2000.
[2] P. Hatzinatoniou and J. Mourjopoluos, " Errors in real-time room acoustics deveberation," J, Audio Engineering Society, vol. 52, no. 9, pp. 883-899, September 2004.
[3] P. Hatzinatoniou and J. Mourjopoluos, " Real-time room equalization based on complex smoothing: Robustness results," Proc. of the 116th AES Convention, 2004.
[4] R. Wilson, " Equalization of loudspeaker drive units considering both on- and off-axis responses," J. Audio Engineering Society, vol. 39, no. 3, pp. 127-139, 1991.
[5] S. J. Elliott and P. A. Nelson, " Multiple-point equalization in a room using adaptive digital filters," J. Audio Engineering Society, vol. 37, no. 11, pp. 899-907, 1989.
[6] F. Talantzis and D. B. Ward, " Robustness of multichannel equalization in an acoustic reverberation environment," Journal of the Acoustic Society of America, vol. 114, no. 2, pp. 833-841, 2003.
[7] F. Talantzis and L. Polymenakos, " Robustness of non-exact multichannel equalization in reveberant environments," in Artificial Intelligence and Innovations 2007: From Theory to Applications, C. Boukis, A. Pnevmatikakis and L. Polymerankos, eds, vol. 247 of IFIP International Federation for Information Processing, pp. 315-321, Springer, Boston, 2007.
[8] M. Sternad and A. Ahlén, " LQ controller design and self-tuning control," in Polynomial Methods in Optimal Control and Filtering, K. Hunt ed., pp.56-92, Peter Peregrinus, London, 1993.
[9] US Patent 7,215,787 .
[10] US Patent 4,683,590 .
[11] US Patent 5,727,066 .
[12] US Patent 5,627,899 .
[13] L-J Brännmark and A. Ahlén, " Spatially Robust Audio Compensation Based on SIMO Feedforward Control", Unpublished manuscript, enclosed as .
[14] US Patent 6,760,451 .
[15] Åström K.J. and B. Wittenmark (1997) Computer-Controlled Systems, 3rd ed. Prentice-Hall, Englewood Cliffs, NJ.
[16] L-J. Brännmark and A. Ahlén, " Robust loudspeaker equalization based on position-independent excess phase modeling," 2007 IEEE ICASSP, Las Vegas, April 2008, To Appear.
[17] M. Sternad and A. Ahlén, " Robust filtering and feedforward control based on probabilistic descriptions of model errors", Automatica, vol. 29, pp. 661-679, 1993.
[18] US Patent 5,511,129 .

APPENDIX 1

Spatially Robust Audio Compensation Based on SIMO Feedforward Control

Lars-Johan Brännmark*, Student Member, IEEE, and Anders Ahlén, Senior Member, IEEE

Abstract

This paper introduces a SIMO feedforward approach to the single-channel loudspeaker equalization problem. Using a polynomial multivariable control framework, a spatially robust equalizer is derived based on a set of room transfer functions (RTFs) and a multiple-point mean square error (MSE) criterion. In contrast to earlier multiple-point methods, the polynomial approach provides analytical expressions for the optimum filter, involving the RTF polynomials and certain spatial averages thereof. A direct use of the optimum solution is however questionable from a perceptual point of view. Despite its multiple-point MSE optimality, the filter exhibits similar, albeit less severe, problems as those encountered in nonrobust single-point designs. First, in the case of mixed phase design it is shown to cause residual "pre-ringings" and undesirable magnitude distortion in the equalized system. Second, due to insufficient spatial averaging when using a limited number of RTFs in the design, the filter is over-fitted to the chosen set of measurement points, thus providing insufficient robustness. A remedy to these two problems is proposed, based on zero clustering and pertinent modifications of the MSE optimal filter. The outcome is a mixed phase compensator with a time domain performance preferable to that of the MSE optimal design.
The authors are with the Signals and Systems Group, Dept. of Engineering Sciences, Uppsala University, PO Box 534. SE-751 21 Uppsala, Sweden. Phone: +46-18-4713071 Fax: +46-18-555096 e-mail: ljb@signal.uu.se; aa@signal.uu.se. L. Brännmark is also with Dirac Research AB, Hansellisgatan 6, SE-754 50 Uppsala, Sweden.

I. INTRODUCTION

The problem of single-channel loudspeaker equalization by the use of digital filters has been extensively studied for about two decades, with an increasing concern in recent years about spatial robustness. In a broad sense, the aim of all audio channel equalization schemes is to remove undesired convolutional distortions introduced by the electro-acoustical signal path of a sound system. In the literature, the work on robustness of equalization essentially falls into three categories. In the first category, the goal of filter design is a complete signal dereverberation at a single position in a room. The subsequent robustness analysis then investigates equalizer performance at other spatial positions, or under slightly modified acoustical circumstances. It is well known that this kind of filter design is highly non-robust and causes severe signal degradation when the receiver position changes [1], and even for fixed receiver position, due to the "weak nonstationarity" of the acoustical paths in the room [2]. In the second category, the design objective is not a complete dereverberation, but rather a reduction of linear distortions, under the constraint that audio performance should not be degraded by changes of listener position. The standard approach in this category is to design a filter based on averaging and/or smoothing of one or several transfer functions and then perform a robustness analysis of the filter [3]. The third category imposes robustness directly on the design by employing a multiple-point error criterion to optimize sound reproduction in a number of spatial positions, either by using measured RTFs [4] or by direct adaptation of the inverse [5]. We mention here parenthetically a fundamentally different multiple-point scenario, where signals are filtered on the receiver side by a unique equalizer at each receiver point. Spatial robustness in this setting has been studied in [6] and [7]. This approach is however not applicable in the pre-compensation setting, where a single filter operates on the input to the system.
In the present paper, the problem formulation relates closest to the third of the above categories. We shall start by defining a multiple-point Mean Square Error (MSE) criterion for spatially robust filter design in a Single-Input Multiple-Output (SIMO) setting. Using a polynomial approach to the multivariable feedforward control problem [8], a linear filter is designed to minimize the multiple-point MSE criterion. The arising equations allow for mixed phase as well as minimum phase inverse design. In contrast to the FIR Wiener/Levinson and LMS approaches used in e.g. [4], [5], the polynomial approach imposes no restrictions on filter order or structure, and the analytical form of the solution is amenable to interpretation in terms of certain spatial averages of the RTFs. MSE optimality does, however, not necessarily imply a good perceptual behavior, which calls for a solution based on refined perceptual considerations. By lack of degrees of freedom in the SIMO setting, ideal equalization in all measurement positions is not possible. Consequently, there will be an equalization error in every position, contributing to the difference between the reconstructed signal and desired signal. Correlations between this error and the desired signal for negative time lags should be limited, as they will be identified by a listener as "pre-ringings" in the equalized system. By inspection of the design equations we develop a method for avoiding the pre-ringing problem, without necessarily resorting to a pure minimum phase inversion. An early version of this approach was introduced by the authors in [9].
The filter design and analysis presupposes an arbitrarily large number of available RTF measurements. For a practical filter design, a spectral smoothing operation has shown to be a valuable complement to the insufficient spatial averaging that arises from using a limited number of RTF measurements. Furthermore, if the sound system subject to equalization has limited bandwidth, some limitation on the filter gain may be necessary in order not to boost frequencies outside the working range of the loudspeaker. Perceptual issues of more intricate nature such as desired tonal coloration etc. can be straightforwardly included in the design. To keep the discussion focused, such issues will however not be considered here.
The paper is organized as follows. Section II formulates the robust audio compensation problem in our SIMO feedforward setting. In Section III the problem is stated and solved mathematically, and the special cases of minimum and mixed phase inversion with ideal target dynamics are studied. In Section IV, qualitative aspects of the filters are investigated for different design scenarios, and some perceptual problems are pointed out. In Sections V and VI, these problems are analyzed and remedies are proposed. In Section VII, the methods of previous sections are evaluated using RTFs acquired in a real room. Finally, Section VIII concludes the paper and points out some directions for further research.

Notation and Terminology

Throughout this paper, we shall use the following notation and terminology: Scalar and vector valued discrete-time signals are denoted by normal and boldface italic letters, like s(k) and s(k), respectively. In the style of [10], transfer functions are represented by polynomial and rational matrices in the backward shift operator q ^-1, defined by q^-1 s(k) = s(k - 1), corresponding to z ^-1 in the frequency domain. All signals and transfer function coefficients are assumed to be real-valued.
Constant matrices are denoted by boldface capital letters as, for example, P. Scalar polynomials are denoted by capital letters in italic as P(q ^-1) = p₀ + p₁ q ^-1 + ... + p_n _Pq ^{-n _P}. Polynomial matrices are denoted by boldface capital letters in italic as P (q ^-1) = P ₀+P ₁ q ^-1+ ···+P_n _Pq ^{-n _P}. Rational matrices are denoted by boldface calligraphic letters as G (q ^-1), and are represented on right Matrix Fraction Description (MFD) form [11]: G = QP ^-1 which for SIMO systems is equivalent to the common denominator form G (q^-1 ) = Q (q ^-1)/P(q ^-1), where Q (q ^-1) is a polynomial matrix and the scalar monic polynomial P(q ^-1) is the least common denominator of all rational elements in G (q ^-1). For scalar rational functions, normal calligraphic letters are used, like G(q ^-1). The arguments q ^-1 and z ^-1 will often be omitted, unless there is a risk for confusion. All polynomials are assumed to have real coefficients. For any polynomial matrix P (q ^-1), or scalar polynomial P(q ^-1), we define their conjugates as $P_{*} (q) = P^{T} (q) = P_{0}^{T} + P_{1}^{T} q + \dots + P_{n_{P}}^{T} q^{n_{P}},$
or $P_{*} (q) = P (q) = p_{0} + p_{1} q + \dots + p_{n_{P}} q^{n_{P}} .$
$ℜ (z)$
and $ℑ (z)$
denote the real and imaginary parts respectively of a complex number z.
A Room Transfer Function (RTF) is a linear time-invariant model of the signal path between the source (sound system input) and receiver (microphone output) in a room. In the general case the RTF between the system input and receiver position i is represented in discrete time by a scalar rational transfer function $H_{i} (z^{- 1}) = B_{i} (z^{- 1}) / A_{i} (z^{- 1});$
i ∈ {1,...,p}, where p is the number of receiver positions. In the sequel, the receiver positions are referred to as control points. We will frequently use transfer operators, e.g. H_i (q ^-1) as a representation of RTFs. For simplicity we will however refer to both as RTFs, or simply transfer functions, and when H_i (q ^-1) is used in this context we mean that q ^-1 is substituted for the complex variable z ^-1. For FIR models (i.e. A_i (q ^-1) = 1), the polynomial notation B_i (q ^-1) is used instead of H_i (q ^-1). The time-domain impulse response related to a transfer function H(z ^-1) is denoted h(k). The complex spatial average model B ₀(q ^-1) refers to the polynomial obtained by taking the coefficient-wise sum of the FIR transfer functions B ₁,..., B_p : $B_{0} (q^{- 1}) = \sum_{i = 1}^{p} B_{i} (q^{- 1}) .$
The RMS spatial average model β(q ^-1) refers to the minimum phase polynomial obtained by spectral factorization of the coefficient-wise sum of the power responses $B$ $_{1 *} B_{1}, \dots,$
$B_{p *} B_{p}$
associated with the FIR models B ₁,..., B_p : $β_{*} (q) β (q^{- 1}) = \sum_{i = 1}^{p} B_{i *} (q) B_{i} (q^{- 1}) .$
The minimum phase equivalent β _i (q ^-1) of an FIR transfer function B_i (q ^-1) is the minimum phase polynomial obtained by spectral factorization of the power response $B_{i *} B_{i} .$
The excess phase part of the same transfer function is the all pass response obtained as B_i (q ^-1)/β _i (q ^-1). A zero cluster is a set of polynomial zeros {z ₁,..., z_p }, located within in a small neighborhood N_∈(z ₀) in the complex plane, where each zero z_i belongs to exactly one RTF B_i (z ^-1). If the region N_∈(z ₀) is sufficiently small, then the RTFs are said to have a near-common zero at z ₀. Zeros outside the unit circle are referred to as excess phase zeros.

II. THE ROBUST AUDIO COMPENSATION PROBLEM

We consider a single-channel setting, where the equalizer filter R(q ^-1) is assumed to operate on a scalar input signal w(k), see Fig. 1. The filtered signal is emitted by a loudspeaker and is received by a listener in one out of (infinitely) many locations in a room. Each receiver location is associated with an individual RTF, and the filter should be designed so as to improve sound reproduction over a whole set of control points. The control points are selected so as to cover a spatial region of hypothetical listener positions, henceforth referred to as the listening region. In deriving the equations, the number of control points, p, is assumed large but finite. Theoretically, a finite p imposes no essential restriction, since by the limited range of wavelengths a discrete grid of points is sufficient to represent the complete sound field within the region of interest. The dense spatial sampling required for such a complete representation is however infeasible in a practical situation, and the optimization will in general be based on a rather low number of RTFs. As we shall see, this restriction can be quite problematic and calls for a solution, if true robustness within the whole listening region is to be obtained.
Now, with each RTF described as a rational function H_i (q ^-1), the signals at the control points can be viewed as the p-dimensional output of a SIMO linear system of dimension p|1, having transfer function matrix $H (q^{- 1}) .$
Similarly, the desired responses $D_{i} (q^{- 1})$
can be stacked in a p|1 matrix $D (q^{- 1}) .$
If the criterion to be minimized is chosen as the sum of the mean squared errors E{|y_i (k)|²}, with y_i (k) being the difference between the received filtered signal and the desired signal, y_i (k) = D_i (q ^-1)w(k) - H _i (q^-1)R(q ^-1)ωw(k), then the problem is equivalent to a SIMO LQ feedforward control problem as depicted in the block diagram of Fig. 1. The sound propagation to the control points is affected by propagation delays of Δ _i samples. While the "true" RTF in position i is $q^{- Δ_{i}} H_{i} (q^{- 1}),$
we shall assume that the individual acoustic delay q ^{-Δ _i} associated with each $H_{i} (q^{- 1})$
is removed prior to the filter design, so that all impulse responses h_i (k) are aligned and start at k = 0. An equivalent but notationally more cumbersome approach would be to include the delays q ^{-Δ _i} in the desired responses D_i (q ^-1).

III. SIMO LQ FEEDFORWARD CONTROL

A. The SIMO Optimum Controller Equations

It is assumed that w(k) is a scalar stationary white noise sequence with zero mean and covariance E{w²(k)} ≙ ψ. The stable rational matrices H (q ^-1) and D (q ^-1), representing respectively the original RTFs and the desired system responses at p spatial control points, are described by right MFD models as $H = B / A = [\begin{matrix} B_{1} \\ ⋮ \\ B_{p} \end{matrix}] \frac{1}{A}; D = D / E = [\begin{matrix} D_{1} \\ ⋮ \\ D_{p} \end{matrix}] \frac{1}{E}$

so that $y (k) = \frac{D}{E} w (k) - \frac{B}{A} Rw (k)$

with A and E being stable monic polynomials. The robust SIMO compensator is defined as the filter which minimizes the sum of the powers of the signals in y (k). That is, the scalar rational filter R(q ^-1) is to be designed so that the criterion $J = E \{{‖ y (k) ‖}_{2}^{2}\} = E \{tr (y (k) y^{T} (k))\}$

is minimized, under the constraints of stability and causality of R(q ^-1). Formulated as above, this problem is readily seen to be a special case of the general MIMO feedforward problem treated in section 3.3 of [8]. Following that derivation and using our specialization of the problem, the optimum compensator filter is given by $R = \frac{QA}{βE}$

where β(q ^-1) is the minimum phase polynomial defined by $β_{*} β = B_{*} B = \sum_{i = 1}^{p} B_{i *} B_{i}$

and Q(q ^-1) = q₀ + q₁ q ^-1 + ··· + q_n _Qq ^-n_Q along with the polynomial L _*(q)=l₀+l₁ q+···+l_n _Lq ^{n _L} constitute the unique solution to the scalar polynomial Diophantine equation $B_{*} (q) D (q^{- 1}) = β_{*} (q) Q (q^{- 1}) + q L_{*} (q) E (q^{- 1})$

with polynomial degrees $n_{Q} = \max \{n_{D}, n_{E} - 1\}; n_{L} = n_{B} - 1.$

B. Optimum Mixed and Minimum Phase Designs

In this subsection we study the filter equations (6)-(8) more closely for the two important special cases of minimum and mixed phase inversion using ideal target dynamics. For clarity of presentation and ease of interpretation we assume that the system H (q ^-1) and the target dynamics D (q ^-1) are polynomial matrices of dimension p|1, containing the RTFs B_i (q ^-1) and target responses D_i (q^-1 ), respectively. Hence A(q ^-1) = E(q ^-1) = 1 in (3) and subsequent equations. This restriction is of no practical importance, since the FIR models B_i (q ^-1) and D_i (q ^-1) are allowed to be of arbitrarily high degree¹.
¹Note that in some situations, it may be more efficient to include A(q ^-1) and E(q ^-1), for example in the modeling of very large or undamped rooms. However, to keep the discussion as clear as possible we use A = E = 1.
We begin the analysis by concluding from (7) that the polynomial β(q ^-1) in the denominator of (6) is identical to the RMS average model (2). Further, we note that if D_i (q ^-1) = q^-d , i.e. the desired response at position i is a pure delay of d samples (in addition to the acoustic delay q ^{-Δ _i} discussed in Section II), then (8) can be rewritten as $\sum_{i = 1}^{p} B_{i *} q^{- d} = β_{*} Q + q L_{*} \Rightarrow q^{- d} B_{0 *} = β_{*} Q + q L_{*}$

where B ₀(q ^-1) is the complex spatial average (1). The delay q^-d represents the number of "future" input signal samples used by the filter. Exchanging q ^-1 for q in (10) and dividing by β(q ^-1) gives the equivalent equation $q^{d} \frac{B_{0} (q^{- 1})}{β (q^{- 1})} = Q_{*} (q) + q^{- 1} \frac{L (q^{- 1})}{β (q^{- 1})} .$
Since β(q ^-1) is minimum phase, we can define the power series Γ(q ^-1) and Λ(q ^-1): $Γ (q^{- 1}) ≜ \frac{B_{0} (q^{- 1})}{β (q^{- 1})} = \sum_{k = 0}^{\infty} γ_{k} q^{- k}; |γ_{k}| < c_{γ} r_{\max}^{k}$
$Λ (q^{- 1}) ≜ q^{- 1} \frac{L (q^{- 1})}{β (q^{- 1})} = \sum_{k = 1}^{\infty} λ_{k} q^{- k}; |λ_{k}| < c_{λ} r_{\max}^{k}$

where c _γ and c _λ are positive constants and r_max < 1 is the maximum radius for any zero of β(z ^-1). Equation. (11) can then be written $\begin{array}{l} q^{d} Γ (q^{- 1}) = Q_{*} (q) + Λ (q^{- 1}) \\ \Leftrightarrow q^{d} \sum_{k = 0}^{\infty} γ_{k} q^{- k} = Q_{*} (q) + \sum_{k = 1}^{\infty} λ_{k} q^{- k} \\ \Leftrightarrow \sum_{k = 0}^{d} γ_{d - k} q^{k} + \sum_{k = 1}^{\infty} γ_{d + k} q^{- k} = Q_{*} (q) + \sum_{k = 1}^{\infty} λ_{k} q^{- k} . \end{array}$
Since Q _*(q) is a polynomial in nonnegative powers of q only, identifying the coefficients for positive and negative powers of q in (14) yields $\begin{matrix} \sum_{k = 1}^{\infty} γ_{d + k} q^{- k} = \sum_{k = 1}^{\infty} λ_{k} q^{- k} \end{matrix} = Λ = q^{- 1} \frac{L (q^{- 1})}{β (q^{- 1})}$
$Q (q^{- 1}) = \sum_{k = 0}^{d} γ_{d - k} q^{- k} .$
We know however from (12) that γ _k is an exponentially decaying sequence, so by increasing the delay d, the coefficients of Λ can be made arbitrarily small. Let the left hand side of (14) be denoted Q̃ _*(q,q ^-1). Then, for the special case when d is very large, $\begin{array}{l} Q (q^{- 1}) & \approx \tilde{Q} (q^{- 1}, q) = \sum_{k = 0}^{d} γ_{d - k} q^{- k} + \sum_{k = 1}^{\infty} γ_{d + k} q^{k} \\ = q^{- d} \sum_{k = 0}^{\infty} γ_{k} q^{k} = q^{- d} B \frac{_{0 *} (q)}{β_{*} (q)} . \end{array}$
With Q̃(q ^-1) ≈ Q̃(q ^-1, q) we mean that Q(q ^-1), which is a polynomial in q ^-1 only, has almost the same impulse response as does Q̃(q ^-1, q), which is a rational function in q and q ^-1. In fact, Q(q ^-1) is the causal part of Q̃(q ^-1, q). The impulse response of Q, when approximated as above, is seen to be the time-reversed and delayed impulse response of the ratio between the complex and RMS spatial averages. Although technically the correct expression for Q is (16), we shall be using the approximation (17) in the following, since it allows the pre-ringing part of the inverse filter to be interpreted as a non-causal filter containing excess phase poles. Using the approximation (17) for Q in (6), and assuming A = E = 1, the optimal compensator filter can be written $R \approx q^{- d} B \frac{_{0 *}}{β_{*}} \frac{1}{β}$

and the equalized system response $H_{i}^{eq} (q^{- 1})$
at position i becomes $H_{i}^{eq} = R B_{i} \approx q^{- d} B \frac{_{0 *}}{β_{*}} \frac{1}{β} B_{i} .$
Note that B _0*/β_* can be expressed as a decaying series in positive powers of q and therefore its impulse response has a noncausal decay.
A second special case of particular interest occurs when d = 0. Equation (9) with degree n_D = 0 then gives that Q must have zero degree and Q = γ₀, with γ₀ obtained from (12), so that $R = \frac{γ_{0}}{β} .$
The equalized system at position i can then be expressed as $H_{i}^{eq} = R B_{i} = γ_{0} \frac{B_{i}}{β}$

whose impulse response decays forward in time only. We shall follow the common terminology of the field and refer to the filters in (18) and (20) as the mixed phase and minimum phase inverse filters, respectively.

IV. QUALITATIVE ASPECTS

Based on the analysis in the previous section, we now state some important qualitative properties of the optimum filter for different scenarios, some of which are of considerable perceptual importance. The system and target dynamics are modeled as in Section III-B. That is, H_i (q ^-1) = B_i (q ^-1), and D_i (q ^-1) = D_i (q ^-1) = q^-d , with d being either zero or very large.

A. Single-Point Mixed Phase Design

In the case of a single-point design, i.e. p = 1 (or if all transfer functions B_i are identical, as may be the case in the far-field of a loudspeaker in an anechoic chamber), then for any i ∈ {1,..., p}, B ₀ = pB_i and $β = \sqrt{p} β_{i} .$
Thus we obtain $Q \approx q^{- d} \sqrt{p} \frac{B_{i *}}{β_{i *}}$
$R \approx q^{- d} \frac{B_{i *}}{β_{i *}} \frac{1}{β_{i}}$
$H_{i}^{eq} \approx q^{- d} \frac{B_{i *} B_{i}}{β_{i *} β_{i}} = q^{- d} = D_{i} .$
We note from (22) that Q approximates an all pass filter (since the magnitude responses of B_i and β_i are identical) scaled by a constant $\sqrt{p},$
and the equalization in (24) is perfect. We recognize R as the time-reversed and delayed excess phase part of B_i in series with the minimum phase inverse 1/β _i . This case is in general of little practical interest and will not be further considered.

B. Multiple-Point Mixed Phase Design

In a multiple-point design (p >> 1) in a normal room, perfect equalization cannot be expected in any point due to the phase and magnitude variability among the RTFs. This variability causes the optimal filter to behave differently from the single-point/anechoic case, and its behavior can be quite problematic from a perceptual perspective. First, Q no longer has all pass character, because the magnitude responses |B ₀(e ^-jw)| and |β(e ^-jw)| differ by more than a constant factor. To see this, suppose that at two separate frequencies w ₀ and w ₁, the magnitudes of all RTFs are equal to one, |B_i (e ^{-jw 0})| = |B_i (e ^{-jw 1})| = 1, while the phases are equal at w ₀, ∠B_i (e ^{-jw 0}) = φ, but random and uniformly distributed at w ₁, ∠B_i (e ^{-jw 1}) = $φ_{i} \sim u [0, 2 π] .$
Then $|β (e^{- j ω_{0}})| = |β (e^{- j ω_{1}})| = \sqrt{p},$
and |B ₀(e ^{-jw 0})| = p. However, due to phase cancellations at w ₁ we have |B ₀(e ^{-jw 1})| << p. Therefore $|Q (e^{- j ω_{0}})| = \sqrt{p},$
but $|Q (e^{- j ω_{1}})| ≪ \sqrt{p},$
and $|R (e^{- j ω_{0}})| = 1 ≫ |R (e^{- j ω_{1}})| .$
Hence, at frequencies where phase variability among the RTFs is large, the MSE optimal filter strategy leads to attenuation of the signal, resulting in a magnitude distortion not suitable for e.g. music listening. Second, the equalized responses $H_{i}^{eq}$
of (19) will contain residual pre-ringings, since the impulse response of $B_{0 *} / β_{*}$
decays noncausally. In Section V we show that the two problems above are interconnected, and a remedy is proposed.

C. Minimum Phase Design

For the case d = 0, the filter
= γ₀/β is minimum phase and has the same character regardless of any possible similarities or dissimilarities among the RTFs. Perfect equalization is obtained only if all B_i are minimum phase and identical. By the strict causality of $H_{i}^{eq}$
in (21), the minimum phase filter is guaranteed to generate no pre-ringing artifacts. It has therefore become common practice in loudspeaker equalizer design to use variants of this filter, with more or less sophisticated processing of the RMS average prior to inversion. It should be noted that there is a significant risk of introducing artificial post-ringings with this type of filter, since all notches in the average frequency response are inverted by minimum phase poles, regardless of whether they were caused by minimum phase zeros or not. We shall be using this filter type for comparison purposes in the experiments in Section VII.

V. TREATMENT OF THE PRE-RINGING PROBLEM

As stated in Section IV-B, the optimum multiple-point mixed phase inverse causes residual pre-ringings in the equalized system, due to the noncausal component B _0*/β_* in (19). In this subsection we analyze this further and propose a remedy to alleviate the pre-ringings. A key issue turns out to be the possible existence of common excess phase zeros, as is shown next.

A. The Origin of Pre-Ringing

Suppose that all RTFs in the listening region share a common factor which is independent of spatial position. Each RTF B_i (q ^-1) can then be decomposed into a robust factor B^r (q ^-1) and a nonrobust factor $B_{i}^{n} (q^{- 1})$
as $B_{i} (q^{- 1}) = B^{r} (q^{- 1}) B_{i}^{n} (q^{- 1}) .$
The corresponding decompositions of the complex and RMS spatial averages then become $B_{0} = \sum_{i = 1}^{p} B^{r} B_{i}^{n} = B^{r} \sum_{i = 1}^{p} B_{i}^{n} = B^{r} B_{0}^{n}$
$β_{*} β = \sum_{i = 1}^{p} B_{*}^{r} B_{i *}^{n} B^{r} B_{i}^{n} = B_{*}^{r} B^{r} \sum_{i = 1}^{p} B_{i *}^{n} B_{i}^{n} = β_{*}^{r} β^{r} β_{*}^{n} β^{n}$

where $B_{0}^{n}$
is the complex spatial average of $B_{1}^{n}, \dots, B_{p}^{n},$
β^r is the minimum phase equivalent of B^r , and βⁿ is the RMS spatial average of $B_{1}^{n}, \dots, B_{p}^{n} .$
Using (27) in (17) and (19) yields $Q \approx q^{- d} B \frac{_{0 *}}{β_{*}} = q^{- d} \frac{B_{*}^{r} B_{0 *}^{n}}{β_{*}^{r} β_{*}^{n}}$
$H_{i}^{eq} \approx q^{- d} \frac{B_{*}^{r} B_{0 *}^{n} B^{r} B_{i}^{n}}{β_{*}^{r} β_{*}^{n} β^{r} β^{n}} = q^{- d} B \frac{_{0 *}^{n} B_{i}^{n}}{β_{*}^{n} β^{n}} .$
The part of (29) that causes the pre-ringings is seen to be $B_{0 *}^{n} / β_{*}^{n},$
which is a factor of Q̃(q ^-1,q) and therefore approximately contained in Q(q ^-1). Note that the noncausally decaying $B_{0 *}^{n} / β_{*}^{n}$
will always occur in the equalized system as soon as the RTFs have non-common zeros, whether they be minimum phase, excess phase or both.

B. A Proposed Improvement

We now consider a modification of
, obtained by avoiding the factor $B$ $_{0 *}^{n} / β_{*}^{n}$
from appearing in Q̃(q ^-1,q). In order to be consistent with (22), where all B_i are identical, i.e. B^r = B_i for any i ∈ {1,..., p}, we include the factor $\sqrt{p}$
in Q. We will thus investigate the precompensator given by $Q \approx q^{- d} \sqrt{p} \frac{B_{*}^{r}}{β_{*}^{r}} \Rightarrow R \approx q^{- d} \sqrt{p} \frac{B_{*}^{r}}{β_{*}^{r}} \frac{1}{β} .$
Clearly, with the proposed modification Q consists of the all pass filter B^r /β ^r which has been time-reversed, scaled with $\sqrt{p}$
and then delayed with q ^-d. We interpret B^r /β ^r as the common excess phase part of the RTFs B ₁,...,B_p . Note that in the quotient B^r /β ^r , all minimum phase zeros of B^r (z ^-1) are cancelled by zeros of β ^r (z ^-1), and the remaining zeros of B^r (z ^-1) and β ^r (z ^-1) are located at reciprocal positions outside and inside the unit circle, respectively. Therefore, Q in (30) can be expressed as $Q (q^{- 1}) \approx q^{- d} \sqrt{p} \frac{{\overline{F}}_{*} (q)}{F_{*} (q)}$

where F _*(q) = f₀+f_1q+···f _nFq^nF is a polynomial constructed from the excess phase zeros of B^r (z ^-1), and F _*(z) is the reciprocal of F _*(q), i.e. F _*(q)=f_{n _F}+f_{n _F -1}q+···+f₀ q ^{n _F}. In general, n_F << n_Br , due to said cancellation of minimum phase zeros. The modified
can be interpreted as the minimum phase inverse of the RMS spatial average in series with an inverse of the common excess phase part. With
as in (30), the equalized system response in position i is $H_{i}^{eq} \approx q^{- d} \sqrt{p} \frac{B_{*}^{r}}{β_{*}^{r}} \frac{1}{β^{r} β^{n}} B^{r} B_{i}^{n} = q^{- d} \sqrt{p} \frac{B_{i}^{n}}{β^{n}}$

which contains no factors that may cause pre-ringings. While discarding the factor $B$ $_{0 *}^{n} / β_{*}^{n}$
from (28) may seem ad-hoc, R in (30) turns out to be the MSE optimal filter if the target responses are chosen as $D_{i} = q^{- d} \sqrt{p} B_{i}^{n} / β^{n} .$
(The perceptual consequences of this choice of target response deserve a further investigation.) Note that the magnitude distortion introduced by Q in (17) of Section III-B is also alleviated by this modification since Q in (30) is all pass. Both of the problems caused by the MSE optimal filter that were discussed in Section IV-B-pre-ringing and magnitude distortion-are thus related to the existence of non-common factors among B₁ ,..., B_p .
One can of course not expect it to occur in practice that all B_i share a truly common excess phase part B^r /β ^r . Nevertheless, in [9] it was demonstrated by the authors that an approximately common excess phase part can be found by detection of zero clusters in the set of RTFs. With the clusters represented by nominal zeros located at the cluster centra, a cluster is classified as invertible if the pre-ringing that results from placing a pole at the nominal zero location is kept below a pre-defined envelope constraint. We now relate this concept to the present work by using the excess phase zeros of B ₀ as nominal zeros.
In order to construct the modified compensator polynomial Q of (30), the near-common all pass factor B^r /β ^r has to be found. This is equivalent to finding the excess phase zeros of B^r . In the case of exactly common zeros, this can be accomplished by discarding all zeros of B ₀ which are not common to all B_i . For this argument to be transferable to the case when zeros are only near-common, we need to know whether the zero clusters of B ₁,..., B_p are represented by zeros in B ₀ which in some sense are close to the zero clusters. An empirical verification of this property is provided in Fig. 2 where the zeros of B ₀ are located approximately at the center of each zero cluster. While a rigorous proof of this property may be quite involved, we motivate it here by a heuristic argument as follows. Let B_i (z ^-1), i ∈ {1,...,p} represent the individual RTFs, and let B ₀(z ^-1) be the coefficient-wise sum of all B_i (z ^-1). Suppose that there is a complex number z ₀ and a small neighborhood N _ε(z ₀) around it, such that each B_i (z ^-1) has a zero z_i within N _ε(z ₀). Define the polynomials G_i by factoring out the zero z_i as B_i = (z - z_i )G_i . Then $\begin{array}{l} B_{0} (z^{- 1}) & = \sum_{i = 1}^{p} B_{i} = \sum_{i = 1}^{p} (z - z_{i}) G_{i} = z \sum_{i = 1}^{p} G_{i} - \sum_{i = 1}^{p} z_{i} G_{i} \\ = (z - \frac{\sum_{i = 1}^{p} z_{i} G_{i}}{\sum_{i = 1}^{p} G_{i}}) \sum_{i = 1}^{p} G_{i} . \end{array}$
Suppose further that the zero cluster contained in N_ε(z ₀) is well separated from all other zeros of B ₁, ..., B_p so that the polynomials G_i do not contain zeros in the vicinity of N_ε(z ₀). Then each G_i can be approximated by a constant for all z ∈ N _ε(z ₀), so that G_i (z ^-1) ≈ g_i , ∀z ∈ N_ε(z ₀). We then obtain $B_{0} (z^{- 1}) \approx (z - \frac{\sum_{i = 1}^{p} z_{i} g_{i}}{\sum_{i = 1}^{p} g_{i}}) \sum_{i = 1}^{p} g_{i}, \forall z \in N_{ϵ} (z_{0})$

i.e. the polynomial B ₀ has a zero which is a weighted average of the zero locations z_i of the individual polynomials B_i . The near-common excess phase zeros of B ₁,..., B_p can hence be found by inspecting each excess phase zero of B ₀ and requiring it to be located within a cluster containing one zero of each B_i . It is intuitively clear that if a zero cluster is small enough, the corresponding zero of B ₀. should be regarded as belonging to B^r in (30). Upon inversion of B^r , the remaining mismatch between B^r and the true zeros of B ₁,..., B_p then causes pre-ringings with negligible amplitudes. In the next subsection we establish a relation between zero cluster size and pre-ringing amplitude.

C. Quantification of Pre-Ringing Error

Suppose that a noncausal filter with transfer function R(z ^-1, z) has been designed to be the inverse of a system H(z ^-1), but with a small mismatch, so that the excess phase poles of R(z ^-1,z) do not completely cancel the excess phase zeros of H(z ^-1). The residual pre-ringing that results can be quantified as follows.
Let a zero of $H (z^{- 1})$
be represented by z ₀ = r ₀ e ^jω0 and a perturbation to this zero by ε = ρe^jθ where r ₀ > 1; 0 < ρ << 1; 0 < ω₀ < π; -π ≤ θ ≤ π. Suppose that H(z ^-1) contains a complex conjugate pair of zeros at z ₀ + ε and z ₀ + ε, so that $H (z^{- 1}) = H_{1} H_{2} = (z - (z_{0} + ϵ)) (z - (\overline{z_{0} + ϵ})) H_{2} .$
Furthermore, suppose that the compensator $R (z^{- 1}, z)$
contains the pole pair z ₀ and z ₀ , $R (z^{- 1}, z) = R_{1} R_{2} = \frac{1}{(z - z_{0}) (z - \overline{z_{0}})} R_{2} .$
The total transfer function of the equalized system thus becomes [9] $\begin{array}{l} H^{eq} (z^{- 1}) & = \frac{(z - (z_{0} + ϵ)) (z - (\overline{z_{0} + ϵ}))}{(z - z_{0}) (z - \overline{z_{0}})} R_{2} H_{2} \\ = (1 - \frac{(ϵ + \overline{ϵ}) z - (ϵ z \overline{_{0}} + \overline{ϵ} z_{0} + ϵ \overline{ϵ})}{(z - z_{0}) (z - \overline{z_{0}})}) R_{2} H_{2} \\ = (1 - 2 ℜ (ϵ) \frac{z - \frac{{|z_{0} + ϵ|}^{2} - {|z_{0}|}^{2}}{2 ℜ (ϵ)}}{(z - z_{0}) (z - \overline{z_{0}})}) R_{2} H_{2} . \end{array}$
Applying the inverse z-transform on each factor in the last line of (37) yields the total impulse response, $h^{eq} (k) = [δ (k) + {Cr}_{0}^{k} \cos (- ω_{0} k + Φ) u (- k)] * r_{2} (k) * h_{2} (k)$

where * denotes convolution, δ(k) is the Kronecker delta function, u(k) is the unit step function, and $Φ = \arctan (\frac{\frac{2 ℜ (ϵ) {|z_{0}|}^{2}}{{|z_{0} + ϵ|}^{2} - {|z_{0}|}^{2}} - ℜ (z_{0})}{ℑ (z_{0})})$
$C = \frac{{|z_{0} + ϵ|}^{2} - {|z_{0}|}^{2}}{{|z_{0}|}^{2} cosΦ} .$

In (37) and (39) we have used the assumptions that Re(ε) ≠ 0, and |z ₀ + ε|² ≠ |z ₀|², which are reasonable for measured data. Equation (38) clearly shows how the pole/zero mismatch ε between H(z ^-1) and R(z ^-1, z) has created a noncausal ringing which affects the total system in a convolutive way.
Suppose now that B_i (z ^-1); i ∈ {1,...,p} represent the set of p RTFs in $H (z^{- 1}),$
each containing M zeros z_im ; m E {1, ..., M}. Furthermore suppose that these zeros are expressed as perturbations, z_im = z _0m + ε _im , of the nominal zeros z _0m; m ∈ {1,..., M _o, M _o + 1,..., M}, where the first M _o nominal zeros are located outside the unit circle in the upper half plane. Once the nominal zeros z _0m = r ₀ _mej ^ω0m and their perturbations ε _im = ρ_ime ^jθim have been determined, equations (38)-(40) with obvious modifications can be used to determine the maximum amplitudes C ₁,..., C _{M o} of the residual pre-ringings caused by placing poles at the nominal zero locations z ₀₁,..., z₀ M _o and their conjugated counterparts z ₀₁,..., z ₀ M _o. We saw in the previous subsection that the excess phase zeros of B ₀ may be used as nominal zeros z _0m. Given that all excess phase zeros are available, what remains is to associate each z_0m with a zero cluster of as small a size a possible.

D. Extraction of Excess Phase Zeros

Suppose that a set of p RTFs B_i (q ^-1); i ∈ {1, ..., p} has been acquired within the listening region. In order to apply the method of the previous subsection, the excess phase zeros of all B_i and of the complex average model B ₀ are required. Considering that the polynomial degree is typically on the order of 10 000-20 000 for FIR models representing full-bandwidth RTFs, finding their zeros is a nontrivial task. However, since only the excess phase zeros of B ₀ and B ₁,..., B_p are sought, they can be found indirectly by identifying the poles of the all pass sequences b̃ ₀(k), b̃ ₁(k),..., b̃ _p (k) defined by the excess phase parts of B ₀, B ₁,..., B_p as ${\tilde{b}}_{0} (k) = \frac{B_{0} (q^{- 1})}{β_{0} (q^{- 1})} δ (k); {\tilde{b}}_{i} (k) = \frac{B_{i} (q^{- 1})}{β_{i} (q^{- 1})} δ (k) .$

The excess phase zeros are then found as the conjugate reciprocals of the pole positions. Note that in b̃ ₀(k) and b̃_i (k) the minimum phase factors of B ₀ and B_i are cancelled by corresponding factors in β₀ and β _i respectively, and the number of poles in b̃ ₀(k) and b̃_i (k) is therefore low compared to the polynomial degrees of B ₀ and B_i . The polynomials β₀ and β _i in (41) can be computed with a suitable spectral factorization algorithm [12], and the poles of b̃ ₀(k) and b̃ _i (k) are then found by performing a model reduction on the systems B ₀/β₀ and B_i /β _i , see e.g. [13].

E. A Pre-Ringing Constraint

With all excess phase zeros given, the next step is to see whether the nominal zeros of B ₀ can be associated with zero clusters of sufficiently small size. "Sufficiently small" here means that the pre-ringing caused by inverting the cluster with a pole at the nominal zero location should not exceed a pre-specified envelope at any control point. If q ^- ^d is the desired system delay included in D(q ^-1), pre-ringings are defined as nonzero values in the equalized system impulse response, |h^eq (k)| > 0, for time indices k < d. We define the maximum tolerable pre-ringing by an exponential envelope constraint as $20 \log_{10} ({Cr}_{0}^{- κ}) < L_{\min}$

where C and r ₀ are as in (38), $κ \in Z^{+}$
is a time constant and L _min is the maximum tolerable pre-ringing level, measured in dB, in the equalized response h^eq (k) at time index k = d - κ. This constraint ensures that the pre-ringing level is at most L _min dB at all time instants prior to k = d - κ Given a nominal zero z ₀, equations (39)-(40) together with the pre-ringing envelope constraint (42) implicitly define a region around z ₀ (and z ₀) within which a zero cluster (and its conjugated counterpart) must be contained in order for z ₀ to be considered a common, robustly invertible, zero of all RTFs. Fig. 3 shows the contours of such regions for different values of z ₀. We note from Fig. 3 that the zero clusters are allowed to be larger if the nominal zero z ₀ is located further away from the unit circle, than when z ₀ is close to the unit circle.

F. Clustering of Near-Common Excess Phase Zeros

We will now describe an algorithm for sorting the excess phase zeros of B ₁,..., B_p into separated clusters, centered around the excess phase zeros of B ₀. The requirement that each cluster must contain exactly one zero from each B_i makes this problem somewhat different from the typical clustering problems encountered in image analysis, data mining etc. No standard off-the-shelf method has been found to be applicable, so the algorithm has been constructed with this specific application in mind. We start with some preliminaries. Suppose that B ₀ contains M _o zeros outside the unit circle in the upper half plane, and that each B_i contains $K_{o}^{i}$
such zeros. Further, assume that $M_{o} \leq K_{o}^{i} \forall i .$
Now arrange these zeros into the sets denoted $Z_{0}$
and $Z_{i}$
respectively: $Z_{0} = \{z_{01} \dots z_{0 M_{o}}\}$
$Z_{i} = \{z_{i 1} \dots z_{i K_{o}^{i}}\}, i \in \{1 \dots p\} .$

The aim of the clustering algorithm is to associate each nominal zero $z_{0 m} \in Z_{0}$
from B ₀ with one zero $z_{ik} \in Z_{i}$
from each B_i . Thereby the zeros are sorted into clusters C_m , defined as $C_{m} = \{z_{1 k_{m}^{1}} z_{2 k_{m}^{}} \dots z_{p k_{m}^{p}}\}, m \in \{1 \dots M_{o}\}$

where the indices $k_{m}^{i}$
determine which of the zeros $z_{i 1}, \dots, z_{i K_{o}^{i}}$
in $Z_{i}$
is to be associated with a certain nominal zero z _0m. We will also make use of a set ${\tilde{Z}}_{0},$
along with an index set µ, defined as $μ = \{μ_{1} \dots μ_{\tilde{M}}\} \subset \{1 \dots M_{o}\}$
${\tilde{Z}}_{0} = \{z_{0 μ_{1}} \dots z_{0 μ_{\tilde{M}}}\} \subset \{z_{01} \dots z_{0 M_{o}}\} = Z_{0}$

where µ is always ordered, i.e. µ _j < µ_i+1, j = 1,..., M̃ - 1. Note that M̃ is the number of elements in µ and ${\tilde{Z}}_{0},$
which varies between different passes through the algorithm. The algorithm is greedy in the sense that, by a principle of "mutually nearest neighbors", it prioritizes dense and well separated clusters instead of minimizing a global criterion based on average distances, as is often the case with other clustering algorithms. The algorithm is described in pseudo code as follows.
Now, with the zeros $z_{i k_{m}^{i}}$
of each cluster C_m expressed as perturbed nominal zeros, $z_{i k_{m}^{i}} =$
$z_{0 m} + ϵ_{k_{m}^{i}},$ one can employ the equations (39), (40) and (42) to decide which zeros in $Z_{0}$
should be included in the inverted common all pass factor $B_{*}^{r} / β_{*}^{r}$
of (30).

VI. SMOOTHING OF THE RMS SPATIAL AVERAGE

In a practical filter design, the number p of transfer function measurements will be limited. Therefore, the RMS spatial average β as defined in (7) will not represent the true RMS average for all possible listener positions, and the filter will be optimal only with respect to the actual measurement positions. For the filter to be truly robust, a method is needed which estimates the true RMS average from a limited number of RTFs. In the present work, this problem has been treated by a smoothing of the frequency response of the finite-sample RMS average, using a 1/6th octave resolution. We motivate this operation by the fact that local irregularities in the RMS frequency response are expected to be smoothed out as the number of RTFs tends to infinity. In Section VII, the benefit of such smoothing is confirmed. A further improved performance is however anticipated with a refined design of the smoothing operation. Another practical issue with influence on the filter design is the bandlimited nature of most loudspeakers. This can be treated by an "amplitude regularization" [14] at the extreme ends of the frequency spectrum, in order to prevent the inverse filter from boosting frequencies outside the working range of the loudspeaker.

VII. A DESIGN EXAMPLE

In this section, we compare the performance of six different equalizer filters designed using the methods of the previous sections. The target dynamics was in all cases set to
= [q ^-d ··· q ^-d] ^T , with either d = 0 or d = 4096, depending on whether the filter is to be minimum or mixed phase. The filters will be referred to with letters from A to F, and they were designed as follows.

A) The MSE optimal mixed phase (d = 4096 samples) filter of equation (18), without any smoothing or regularization of the RMS spatial average.
B) The modified mixed phase (d = 4096 samples) filter of equation (30), without smoothing or regularization.
C) The minimum phase (d = 0) filter of equation (20), without smoothing or regularization.
D) Same as filter A, but with smoothing and regularization of the frequency response of the RMS spatial average prior to computing β. Smoothing resolution was 1/6th octave, and regularization was used below 30 Hz and above 20 kHz.
E) Same as filter B but with the same smoothing and regularization as filter D.
F) Same as filter C but with the same smoothing and regularization as filter D. Filter F represents the "standard" minimum phase approach to loudspeaker equalization.

A. Methods for Evaluation

The performance of a filter will be assessed by studying simulated responses of the equalized system at different control points. These responses are obtained by applying the filter R(q ^-1) to the impulse responses of the RTFs in question: $h_{n}^{eq} (k) = R (q^{- 1}) B_{n} (q^{- 1}) δ (k), n \in \{1 \dots N\} .$

Hence we here rely on the assumption of linearity and time-invariance of the true system, i.e. that the simulated equalized response $h_{n}^{eq} (k)$
is equal to that obtained by a real RTF measurement of the system at position n, using a test signal pre-filtered with R(q ^-1). Robustness is assessed by comparing the performance for two different sets of RTFs. The first set is the design set containing p RTFs which represent the control points that were used for filter design. The second set is the validation set, representing control points within the listening region, but spatially separated from the design set, see Fig. 4. Such a comparison indicates to what extent the filters are over-fitted to the design points. Since our proposed modified design is based primarily on a time-domain argument (avoidance of pre-ringings), the assessment will focus on the time-domain behavior of the filters. We will however start by presenting the average RMS frequency responses of the system, before and after equalization. For graphical evaluation of the time domain properties, we shall use the average Schroeder decay sequence D(k), the average energy step response, or energy build-up, S(k) and the impulse response maximum level envelope L(k), $D (k) = 10 \log_{10} (\frac{1}{N} \sum_{n = 1}^{N} \sum_{l = k}^{M - 1} \frac{h_{n}^{2} (l)}{\sum_{m = 0}^{M - 1} h_{n}^{2} (m)})$
$S (k) = \frac{1}{N} \sum_{n = 1}^{N} \sum_{l = 0}^{k} \frac{h_{n}^{2} (l)}{\sum_{m = 0}^{M - 1} h_{n}^{2} (m)}$
$L (k) = 20 \log_{10} (\max_{n} |h_{n} (k)|)$

defined in (49), (50) and (51) respectively. The Schroeder and energy build-up curves were introduced in [15] and [16] respectively. Here h_n (k); k ∈ {0,..., M - 1} is an impulse response of length M in microphone position n; n ∈ {1, ..., N}. Prior to computation of D(k), S(k) and L(k), all responses are time-aligned and normalized so that max |h_n (k)| = |h_n (k ₀)| = 1 for some time instant k = k ₀. While L(k) is useful as a worst case presentation of pre- or post-ringing problems, S(k) and D(k) indicate how good are the transient properties of the system. In order for a comparison of systems with different pre-ringing behavior to be feasible, a further alignment of the curves D(k) and S(k) is needed. We have chosen to define the starting time, k = 0, of S(k) so that k = 1 occurs at the sample where S(k) for the first time reaches above 5% of its steady state value. For D(k), we define k = 0 so that k = 1 occurs at the sample where the decay for the first time reaches below -0.5 dB. It is sometimes instructive to see how the curves D(k) and S(k) behave in narrow frequency bands, and we shall therefore complement the full frequency band presentations with low pass filtered versions, with a cutoff frequency of 320 Hz.

B. Experimental Conditions

In a room of dimensions 4.5×6×2.6 m, with an average distance between loudspeaker and microphones equal to 2.5 m, 9 measurement positions for filter design (p = 9), and 9 positions for validation were selected according to Fig. 4. This microphone configuration was designed to cover the typical head movements of a normal listener. The RTFs were acquired using a pink-colored random phase multisine signal [17, Chapter 13] with a period time of 3 seconds. The FIR models so obtained were truncated to a length of 0.408 seconds, or 18 000 coefficients at a sampling frequency of 44 100 Hz. This model order is motivated by the reverberation time T₆₀ of the room which is slightly less than 0.4 seconds at low frequencies.
Filters A to F were designed as described in the beginning of this section. The parameters in the pre-ringing constraint (42) were set to L _min = -60 dB and κ = 220 samples. The minimum phase polynomials β(q ^-1), β ₀(q ^-1) and β_i (q ^-1) were obtained by spectral factorization [12], and the poles of the sequences b̃ ₀(k) and b̃_i (k) in (41) were identified using a Hankel matrix based model reduction technique as described in [13]. The accuracy of this method for finding excess phase zeros has been found to be reasonably good when compared to a brute-force polynomial rooting approach. A deeper study of the accuracy of this method is unfortunately beyond the scope of the present paper.

C. Results

In this subsection, we present graphically the time and frequency domain performance of the filters A to F. We begin by stating some properties that are evident from the frequency responses of Fig. 5.

Filters A-C perform as desired only in the design points. Although the general trends in the frequency responses are corrected in the validation points also, the filters seem to cause an increased jaggedness of the curves at high frequencies.
The "attenuation property" of the MSE optimal filter, discussed in subsection IV-B, is evident in the frequency responses of filters A and D. The deep notches at 190, 280, 400 and 600 Hz respectively indicate a large phase variability at those frequencies among the RTFs in the listening region.
In the frequency region between about 30 and 200 Hz, the "unsmoothed" filters A-C perform better than filters D-F, even in the validation points. This suggests that 1/6th octave smoothing is too coarse at those frequencies, motivating a more flexible smoothing operation.
The most desirable overall frequency domain performance is exhibited by filters E and F, which flatten out the response without adding any strange properties to the curves. A further discrimination between filters E and F is not possible based on Fig. 5, since they differ only by an all pass factor.

Next, we turn to studying the time domain properties of the filters. The curves L(k) in Fig. 6 obviously reveal some important properties not visible in Fig. 5. We summarize the details provided by Fig. 6:

The pre-ringings caused by filters A and D are unacceptably high (-40 dB at 20 msec before the maximum peak).
The ratio between the maximum peak and the lower levels seems to be improved, both in the design and validation points, by all filters except filter C. There is however slightly less improvement in the validation points for filters A-C than for filters D-F.
Best overall performance is exhibited by filters B and E, which cause only a very low level of pre-ringing, while substantially amplifying the maximum peak in the responses.

So far, our graphical evaluation suggests filters E and F as the best candidates for a perceptually acceptable loudspeaker compensation, since they are the only filters without any immediately objectionable properties. However, provided that its low-level pre-ringings can be tolerated, filter E seems to possess the most preferable time-domain properties. This is confirmed by a study of the Schroeder decays and energy step responses in figures 7-10. We conclude this section by commenting on the behavior observed in these figures. It should be noted that the scales on the axes of the diagrams in figures 7-10 have been selected so as to display the most interesting parts of the responses in a reasonable resolution.

The top left diagram of Fig. 7 suggests an intuitively appealing ranking of the filters A-C: All of the filters A-C seem to improve the original system, with the fastest energy build-up being provided by filter A, closely followed by filter B, while filter C causes only a moderate improvement. This behavior is however not maintained in the other diagrams. The bottom left and right diagrams show that it no longer holds at low frequencies, because the pre-ringing part due to filter A obviously contains a considerable part of the total low-frequency energy. In the validation points, the pre-ringing problem is evident also in the full bandwidth case. Moreover, in the validation points the unequalized response has, at times, better performance than any of the equalized responses. In particular, in the bottom right diagram at about 23 msec, the original response "catches up" on the step responses produced by filters B and C. By increasing the sound energy in the late parts of the impulse responses, the filters have thus caused artificial post-ringings in the validation points.
Fig. 8 provides essentially the same information as Fig. 7, although the post-ringing problems introduced by filters B and C are even more evident here.
Based on figures 9 and 10, filter D can be ruled out due to its severe pre-ringings at low frequencies. Filter F improves upon the original response everywhere except in the first few samples of the full bandwidth case. Filter E is seen to improve the original response everywhere. It is considerably better than filter F in the earliest parts (0.0-0.3 msec) of the full bandwidth responses, and throughout the low-frequency responses.

VIII. CONCLUSION

A new method for robust mixed-phase audio compensation has been presented. By the use of polynomial multivariable control techniques and a SIMO MSE criterion, analytical expressions for a spatially robust filter were obtained. It was shown that the optimum mixed phase MSE solution involves two kinds of spatial averages, here named the complex and RMS averages respectively, of which the latter is commonly used in minimum phase equalizer design. Due to perceptual shortcomings of the optimum mixed phase MSE filter, a refined mixed phase design was proposed and experimentally shown to possess time domain qualities preferable to those of the MSE optimal mixed and minimum phase filters. It is our opinion that this result motivates a revision of the widespread conclusion that excess phase properties of the RTFs must be neglected in a robust equalizer design.
In order to keep the presentation transparent, RTFs were represented with FIR models H_i (q ^-1) = B_i (q ^-1) in most of the analysis in sections III-V and in the evaluations in section VII. The results and interpretations regarding e.g. spatial averages and clustering of near-common zeros can however be shown to be valid for the general IIR model H_i (q ^-1) = B_i (q-1)/A_i (q ^-1). In particular, the results hold for the common acoustical pole and zero model (CAPZ) [18], where H_i (q ^-1) = B_i (q ^-1 )/A, for a common pole polynomial A. The inverse filters derived in the present work however differ significantly- from that proposed in [18], which consists only of the common denominator, R(q ^-1) = A(q ^-1).
In subsection V-B, we mentioned that the modified mixed phase compensator would be MSE optimal for a particular choice of target responses, namely $D_{i} = q^{- d} \sqrt{p} B_{i}^{n} / β^{n},$
which involves the nonrobust part $B_{i}^{n}$
of the RTF B_i and the RMS average βⁿ of all the nonrobust parts. The role of $B_{i}^{n}$
in this context can be explained by the fact that it is the nonrobust dynamics that causes the pre-ringing and consequently, a strategy to avoid pre-ringings is to leave the nonrobust dynamics untreated in all control points. The choice of β ⁿ as denominator of D_i is harder to motivate and it is not obvious how to change it or remove it from D_i , since this target response is obtained indirectly without knowledge of the full robust/nonrobust decomposition (25). A further investigation of this kind of "clever target assignment" requires a complete knowledge of the non-common and near-common factors among B ₁,..., B_p . Recent results from the field of approximate greatest common divisors (AGCDs) of polynomials [19] may prove fruitful here.
We argued in Section VI that an accurate estimate of the true power response average in the region, based only on a few measurements, is required for a practically useful filter design. We also saw in Section VII that our solution-a 1/6th octave smoothing of the RMS spectrum-was helpful, although far from optimal. A more flexible smoothing operation, taking more acoustical information into account, would probably improve the filter performance.
Finally, we emphasize that the applicability of our proposed mixed-phase method, and its superiority to a standard minimum phase design, heavily depends on the existence of a near-common excess phase part among the RTFs in the listening region. In an arbitrary acoustic environment, there is of course nothing that guarantees the existence of such a common part. Our experience so far has however indicated that it may exist under quite general circumstances. It is an interesting topic for further research to reach a better understanding of the conditions for its existence.

ACKNOWLEDGMENT

This research was supported in part by Dirac Research AB, Uppsala, Sweden.

REFERENCES

[1] B. D. Radlovic, R. C. Williamson, and R. A. Kennedy, "Equalization in an acoustic reverberant environment: robustness results," IEEE Transactions on Speech and Audio Processing, vol. 8, no. 3, pp. 311-319, May 2000.
[2] P. Hatziantoniou and J. Mourjopoulos, "Errors in real-time room acoustics dereverberation," J. Audio Eng. Soc., vol. 52, no. 9, pp. 883-899, September 2004.
[3] P. Hatziantoniou and J. Mourjopoulos, "Real-time room equalization based on complex smoothing: Robustness results," in Proceedings of the 116th AES Convention. AES, 2004.
[4] R. Wilson, "Equalization of loudspeaker drive units considering both on- and off-axis responses," J. Audio Eng. Soc., vol. 39, no. 3, pp. 127-139, 1991.
[5] S. J. Elliott and P. A. Nelson, "Multiple-point equalization in a room using adaptive digital filters," J. Audio Eng. Soc., vol. 37, no. 11, pp. 899-907, 1989.
[6] F. Talantzis and D. B. Ward, "Robustness of multichannel equalization in an acoustic reverberant environment," The Journal of the Acoustical Society of America, vol. 114, no. 2, pp. 833-841, 2003.
[7] F. Talantzis and L. Polymenakos, "Robustness of non-exact multi-channel equalization in reverberant environments," in Artificial Intelligence and Innovations 2007: From Theory to Applications, C. Boukis, A. Pnevmatikakis, and L. Polymenakos, Eds., vol. 247 of IFIP International Federation for Information Processing, pp. 315-321. Springer, Boston, 2007.
[8] M. Sternad and A. Ahlén, "LQ controller design and self-tuning control," in Polynomial Methods in Optimal Control and Filtering, K. Hunt, Ed., pp. 56-92. Peter Peregrinus, London, 1993.
[9] L. Brännmark and A. Ahlén, "Robust loudspeaker equalization based on postition-independent excess phase modeling," in 2008 IEEE International Conference on Acoustics, Speech, and Signal Processing, Proceedings., 2008.
[10] K. J. Åström and B. Wittenmark, Computer-Controlled Systems: theory and design, Prentice-Hall, Upper Saddle River, NJ, 1997.
[11] T. Kailath, Linear Systems, Prentice-Hall, Englewood Cliffs, NJ.
[12] A. H. Sayed and T. Kailath, "A survey of spectral factorization methods," Numerical linear algebra with applications, vol. 8, no. 6-7, pp. 467-496, 2001.
[13] S. Kung, "A new identification and model reduction algorithm via singular value decomposition," in Proceedings of the 12th Asilomar Conference on Circuits, Systems and Computing., 1978, pp. 705-714.
[14] P. Craven and M. Gerzon, "Practical adaptive room and loudspeaker equaliser for hi-fi use," in The Proceedings of the AES DSP UK Conference. AES, 1992, pp. 121-153.
[15] M. R. Schroeder, "New method of measuring reverberation time," The Journal of the Acoustical Society of America, vol. 37, no. 3, pp. 409-412, 1965.
[16] E. A. Robinson, Statistical Communication and Detection, Griffin, London, 1967.
[17] L. Ljung, System Identification - Theory for the User, Prentice-Hall, 2 edition, 1999.
[18] Y. Haneda, S. Makino, and Y. Kaneda, "Multiple-point equalization of room transfer functions by using common acoustical poles," IEEE Transactions on Speech and Audio Processing, vol. 5, no. 4, pp. 325-333, July 1997.
[19] R. M. Corless, S. M. Watt, and Lihong Zhi, "QR factoring to compute the GCD of univariate approximate polynomials," IEEE Transactions on Signal Processing, vol. 52, no. 12, pp. 3394-3402, 2004.

	Lars-Johan Brännmark (S'07) was born in Luleå, Sweden in 1974. He received the M.S. degree from Uppsala University in 2004. He holds a Diploma in Sound Engineering from Piteå School of Music, Luleå University of Technology and has studied Electrical Engineering at Purdue University, Indiana, and Musicology at Uppsala University, for one year each. He has a professional background as a Sound Engineer at the Swedish National Radio and as a Research Engineer at Dirac Research AB. He is currently working towards a PhD degree in Signal Processing at the Department of Engineering Sciences, Uppsala University. His research interests are in the field of signal processing for audio applications

	Anders Ahlén (S'80-M'84-SM'90) is full professor and holds the chair in Signal Processing at Uppsala University where he is also the head of the Signals and Systems Group. He was born in Kalmar, Sweden, and received the PhD degree in Automatic Control in 1986 from Uppsala University. He was with the Systems and Control Group, Uppsala University, from 1984-1989 as an Assistant Professor and from 1989-1992 as an Associate Professor in Automatic Control. During 1991 he was a visiting researcher at the Department of Electrical and Computer Engineering, The University of Newcastle, Australia. In 1992 he was appointed Associate Professor of Signal Processing at Uppsala University. From 1998-2004 he was the editor of Signal and Modulation Design for the IEEE Transactions on Communications. During 2001-2004 he was the CEO of Dirac Research AB, a company offering state-of-art audio signal processing solutions. He is currently the chairman of the board of the same company. His research interests, which include Signal Processing, Communications and Control, are currently focused on Signal Processing for Wireless Communications, Wireless Systems Beyond 3G, Wireless Sensor Networks, and Audio Signal Processing.

Claims

A method for designing a discrete-time audio precompensation filter (200) based on a Single-Input Multiple Output linear model ( H ) that describes the dynamic response of an associated sound generating system at p > 1 listening positions, for which said dynamic response differs for at least two of these listening positions, said Single Input Multiple Output linear model having p individual scalar models from the single input to the p outputs of the linear model, characterized by
providing information representative of n non-minimum phase zeros {z_i } that are outside of the stability region lzl = 1 in the complex frequency domain, where 1 ≤ n < m, with m being the smallest number of zeros outside |z| = 1 of any of the p individual scalar models from the single input to the p outputs of the linear model H , said non-minimum phase zeros being selected such that their inversion by a given compensator filter, for evaluation purposes, results in pre-ringings of the compensated impulse response that are smaller than a prespecified limit;
determining said precompensation filter as the product of at least two scalar dynamic systems, said at least two scalar dynamic systems being represented by:
i) an inverse of a characteristic scalar magnitude response in the frequency domain that represents the power gains at all or a subset of the p listening positions according to the model H ;

ii) a causal Finite Impulse Response (FIR) filter of user-specified degree d, wherein said FIR filter has coefficients corresponding to a causal part of a delayed non-causal impulse response that is based on said information representative of n non-minimum phase zeros.
The method according to claim 1, wherein said step of providing information representative of n non-minimum phase zeros comprises the step of forming a design polynomial F(z) having said n zeros, and said causal FIR filter, of user-specified degree d, is represented by a polynomial denoted Q(q^-1 ) where q^-1 is the backward shift operator q^-1v(t) = v(t-1), and wherein said FIR filter is obtained by
- forming a non-causal all-pass filter having the polynomial F(q) as denominator, where the forward shift operator q defined as qv(t) = v(t+1), is substituted for the complex variable z,

- multiplying the non-causal impulse response of this non-causal all-pass filter by a time delay of d samples, to obtain a delayed non-causal impulse response,

- setting Q(q^-1) equal to a causal part of said delayed non-causal impulse response.
The method according to claim 2, wherein said FIR filter Q(q^-1) approximately inverts the non-minimum phase dynamics of a causal linear discrete-time dynamic system that has F(q) as its transfer function numerator.
The method according to claim 2 or 3, wherein a stable scalar design model is formed as the complex spatial average of the p individual scalar dynamic models in H that describe the separate responses at the p listening positions, and the zeros of the design polynomial F(z) are selected as a subset of the non-minimum phase zeros (zeros outside |z| = 1) of said scalar design model.
The method according to claim 4, wherein the set of zeros {z_i } of the scalar design model that are selected to be zeros of the design polynomial F(z) are selected by judging the magnitudes of the coefficients prior to a delay equal to d of the impulse response vector of the compensated response H(q^-1)R(q^-1) for the p measurement positions, where R(q^-1) denotes said precompensation filter.
The method according to claim 5, wherein each zero of the transfer function of the scalar design model is associated with a cluster of zeros, where each zero of a cluster is equal or approximately equal to a zero of an individual transfer function of the model H and where said clusters are used to estimate the magnitudes of the coefficients prior to a delay d of the impulse response vector of the compensated response to the p measurement positions.
The method according to claim 1, wherein the information representative of n non-minimum phase zeros represents an approximate common factor of the p transfer function numerator polynomials in the model H that is represented by a vector of p rational discrete-time transfer functions.
The method according to claim 1, wherein the characteristic scalar magnitude response is obtained by
- forming a Root Mean Square (RMS) average model by summing the power spectral densities of the p rational discrete-time transfer functions and adding an optional frequency-dependent scalar regularization parameter, and

- optionally performing frequency-domain smoothing of the frequency response curve of said RMS average model.
The method according to claim 2, where the polynomial Q(q^-1) is a scaled by a scalar, c, resulting in the design equation: $c q^{- d} \frac{(f_{0} + f_{1} q + \dots + f_{n} q^{n})}{(f_{n} + f_{n - 1} q + \dots + f_{0} q^{n})} = m_{- d} q^{- d} + m_{- d + 1} q^{- d + 1} + \dots + m_{0} + m_{1} q + m_{2} q^{2} + ....$

where fi are the coefficients of the design polynomial F(q) = f_n + f _n-1 q +...+f ₀ qⁿ and m_i are the coefficients of the delayed non-causal impulse response, and setting Q(q^-1) equal to the causal part of this delayed non-causal impulse response gives: $Q (q^{- 1}) = m_{- d} q^{- d} + m_{- d + 1} q^{- d + 1} + \dots + m_{0} .$
The method according to claim 1, wherein the linear model ( H ) that describes the dynamic response of an associated sound generating system is determined based on measurements of sound at said p > 1 listening positions, said sound being produced by said sound generating system, and said step of calculating said precompensation filter as the product of at least two scalar dynamic systems comprises the step of determining corresponding filter parameters, and said audio precompensation filter is created by implementing the determined filter parameters in a filter structure.
The method according to claim 10, wherein said filter structure is embodied together with said associated sound generating system so as to enable generation of sound influenced by said audio precompensation filter.
A system (100; 200) for designing a discrete-time audio precompensation filter (200) based on a Single-Input Multiple Output linear model ( H ) that describes the dynamic response of an associated sound generating system at p > 1 listening positions, for which the said dynamic response differs for at least two of these listening positions, said Single Input Multiple Output linear model having p individual scalar models from the single input to the p outputs of the linear model, characterized by
means for providing information representative of n non-minimum phase zeros {z_¡ } that are outside of the stability region |z| = 1 in the complex frequency domain, where 1 ≤ n < m, with m being the smallest number of zeros outside |z| = 1 of any of the p individual scalar models from the single input to the p outputs of the linear model H , said non-minimum phase zeros being selected such that their inversion by a given compensator filter, for evaluation purposes, results in pre-ringings of the compensated impulse response that are smaller than a prespecified limit;
means for determining a characteristic scalar magnitude response in the frequency domain that represents the power gains at all or a subset of the p listening positions according to the model H ;
means for determining, based on said information representative of n non-minimum phase zeros, a causal Finite Impulse Response (FIR) filter having coefficients corresponding to a causal part of a delayed non-causal impulse response, said causal FIR filter being of user-specified degree d;
means for determining, for design purposes, said precompensation filter as the product of at least two scalar dynamic systems, said at least two scalar dynamic systems being represented by:
i) an inverse of said characteristic scalar magnitude response; and

ii) said causal Finite Impulse Response (FIR) filter.
The system according to claim 12, wherein said means for providing information representative of n non-minimum phase zeros comprises means for forming a design polynomial F(z) having said n zeros, and said causal FIR filter, of user-specified degree d, is represented by a polynomial denoted Q(q^-1) where q^-1 is the backward shift operator q^-1v(t) = v(t-1), and wherein said means for determining a FIR filter comprises:
- means for providing a non-causal all-pass filter having the polynomial F(q) as denominator, where the forward shift operator q defined, as qv(t)= v(t+1), is substituted for the complex variable z,

- means for multiplying the non-causal impulse response of this non-causal all-pass filter by a time delay of d samples, to obtain a delayed non-causal impulse response,

- means for setting Q(q^-1) equal to a causal part of said delayed non-causal impulse response.
The system according to claim 13, wherein said FIR filter Q(q^-1) is configured to approximately invert the non-minimum phase dynamics of a causal linear discrete-time dynamic system that has F(q) as its transfer function numerator.
The system according to claim 13 or 14, comprising means for providing a stable scalar design model as the complex spatial average of the p individual scalar dynamic models in H that describe the separate responses at the p listening positions, and wherein the zeros of the design polynomial F(z) are selected as a subset of the non-minimum phase zeros (zeros outside lzl = 1) of said scalar design model.
The system according to claim 15, wherein said means for providing a design polynomial is configured for selecting the set of zeros {z_i } of the scalar design model that are selected to be zeros of the design polynomial F(z) by judging the magnitudes of the coefficients prior to a delay equal to d of the impulse response vector of the compensated response H(q^-1)R(q^-1) for the p measurement positions, where R(q^-1) denotes said precompensation filter.
The system according to claim 16, wherein each zero of the transfer function of the scalar design model is associated with a cluster of zeros, where each zero of a cluster is equal or approximately equal to a zero of an individual transfer function of the model H and where said clusters are used to estimate the magnitudes of the coefficients prior to a delay d of the impulse response vector of the compensated response to the p measurement positions.
The system according to claim 12, wherein said means for providing information representative of n non-minimum phase zeros is configured to provide information that represents an approximate common factor of the p transfer function numerator polynomials in the model H that is represented by a vector of p rational discrete-time transfer functions.
The system according to claim 12, wherein said means for determining a characteristic scalar magnitude response comprises
- means for forming a Root Mean Square (RMS) average model by summing the power spectral densities of the p rational discrete-time transfer functions and adding an optional frequency-dependent scalar regularization parameter, and

- means for optionally performing frequency-domain smoothing of the frequency response curve of said RMS average model.
The system according to claim 13, wherein means for determining a causal FIR filter comprises means for scaling the polynomial Q(q^-) by a scalar, c, resulting in the design equation: $c q^{- d} \frac{(f_{0} + f_{1} q + \dots + f_{n} q^{n})}{(f_{n} + f_{n - 1} q + \dots + f_{0} q^{n})} = m_{- d} q^{- d} + m_{- d + 1} q^{- d + 1} + \dots + m_{0} + m_{1} q + m_{2} q^{2} + ....$

where f_i are the coefficients of the design polynomial F(q) = f_n + f _n-1 q +...+f ₀ qⁿ and m_i are the coefficients of the delayed non-causal impulse response, and said means for setting Q(q^-1) equal to the causal part of this delayed non-causal impulse response gives: $Q (q^{- 1}) = m_{- d} q^{- d} + m_{- d + 1} q^{- d + 1} + \dots + m_{0} .$
The system according to claim 12, further comprising means for determining the linear model ( H ) that describes the dynamic response of the associated sound generating system based on measurements of sound at said p > 1 listening positions, said sound being produced by said sound generating system, and said means for calculating said precompensation filter as the product of at least two scalar dynamic systems comprises means for determining corresponding filter parameters, wherein said audio precompensation filter is created by implementing the determined filter parameters in a filter structure.
The system according to claim 21, wherein said filter structure is embodied together with said associated sound generating system so as to enable generation of sound influenced by said audio precompensation filter.
A computer program product for designing, when running on a computer system (100; 200), a discrete-time audio precompensation filter (200) based on a Single-Input Multiple Output linear model ( H ) that describes the dynamic response of an associated sound generating system at p > 1 listening positions, for which the said dynamic response differs for some of these listening positions, said Single Input Multiple Output linear model having p individual scalar models from the single input to the p outputs of the linear model, characterized by
program means for providing information representative of n non-minimum phase zeros {z_i } that are outside of the stability region |z| =1 in the complex frequency domain, where 1 ≤ n < m, with m being the smallest number of zeros outside |z| = 1 of any of the p individual scalar models from the single input to the p outputs of the linear model H ;
program means for determining a characteristic scalar magnitude response in the frequency domain that represents the power gains at all or a subset of the p listening positions according to the model H ;
program means for determining, based on said information representative of n non-minimum phase zeros, a causal Finite Impulse Response (FIR) filter having coefficients corresponding to a causal part of a delayed non-causal impulse response, said causal FIR filter being of user-specified degree d; and
program means for determining, for design purposes, said precompensation filter as the product of at least two scalar dynamic systems, said at least two scalar dynamic systems being represented by:
i) an inverse of said characteristic scalar magnitude response; and

ii) said causal Finite Impulse Response (FIR) filter.
An audio precompensation filter designed by using the method according to claim 1.
An audio system comprising a sound generating system and an audio precompensation filter in the input path to said sound generating system, wherein said audio precompensation filter is designed by using the method according to claim 1.
A digital audio signal generated by a precompensation filter designed by using the method according to claim 1.