Apparatus and method for reproducing a sound field with a loudspeaker array controlled via a control volume

Abstract

Method and, apparatus for implementing the method, the method comprising determining control signal data for an array of loudspeakers, the control signal data being such as to control the loudspeakers to produce a desired sound field associated with an audio signal, the method comprises determining control signal data for different frequency components of the desired sound field in respect of respective different positions in a listening volume of the loudspeaker array, wherein determination of the control signal data comprises sampling the desired sound field at the surface of a control volume (V).

Description

CROSS REFERENCE TO RELATED APPLICATIONS

This application is the U.S. national stage application under 35 U.S.C. §371 of co-pending International Application No. PCT/GB2009/051292, filed Oct. 1, 2009 and designating the U.S., which published as WO 2010/038075 A2 on Apr. 8, 2010, and which claims the benefit of Great Britain Application Serial No. 0817950.9, filed Oct. 1, 2008. Each of the foregoing patent applications and patent application publications are expressly incorporated herein by reference in their entireties.

TECHNICAL FIELD

The present invention relates to an apparatus and method for sound reproduction.

BACKGROUND

In the prior art, systems are known for the recording and the reproduction of sound. Such systems provide, with an array of loudspeakers, the physical reconstruction of a desired sound field over a region of space. The sound field generated by the loudspeakers should give to the listeners located in the listening area the realistic perception of a desired virtual sound source or of a virtual sound scene.

One well known technology of this kind is Wave Field Synthesis disclosed in patent applications US 2006/0098830 A1, WO2007/101498 A1, US 2005/0175197 A1, US 2006/0109992 A1, for example. The technology uses the Kirchhoff Helmholtz equation which implies, in theory, the use of both dipole-like and monopole-like secondary sources (the loudspeakers), the strength of which (that are proportional to the loudspeaker signals) is explicitly given by the value of the sound field on the integration contour and its normal derivative, respectively.

Another known technology is Ambisonics, for example as disclosed in U.S. Pat. No. 5,757,927—1998, High Order Ambisonics as disclosed in US 2006/0045275 A1 and another technology disclosed in US 2005/0141728 A1. The theoretical formulation of these methods involve a large use of cylindrical or spherical harmonics and Legendre polynomials, in the same way that the use of sines and cosines or complex exponentials arise in the theoretical formulation of any technology based on the traditional Fourier transform. These prior art technologies involve the use of a matrix-based processing to encode the recorded signals and generate an intermediate format, and a matrix-based system to decode this intermediate format and generate the driving signals for the loudspeakers.

We seek to provide an improved apparatus and method for reproduction of sound.

SUMMARY

According to one aspect of the invention there is provided a method of determining control signal data for an array of loudspeakers, the control signal data being such as to control the loudspeakers to produce a desired sound field associated with an audio signal, the method comprises determining control signal data for different frequency components of the desired sound field in respect of respective different positions in a listening volume of the loudspeaker array, wherein determination of the control signal data comprises sampling the desired sound field at the surface of a control volume.

According to another aspect of the invention there is provided sound reproduction apparatus for processing an audio signal, the apparatus configured to output control signal data for an array of loudspeakers to produce a desired sound field associated with the audio signal, wherein the apparatus configured to determine the control signal data for different frequency components of the desired sound field in respect of respective different positions in a listening volume of the loudspeaker array, wherein determination of the control signal data comprises sampling the desired sound field at the surface of a control volume.

A further aspect of the invention relates to a signal processor configured to process the audio signal of the above aspects of the invention and output the control signal data. The signal processor may be configured by suitable machine-readable instructions, and the instructions may be realised in the form of a signal or a data carrier device.

BRIEF DESCRIPTION OF THE DRAWINGS

Various embodiments of the invention will now be described by way of example only, with reference to the following drawings in which:

FIG. 1 shows a co-ordinate system,

FIG. 2 shows of different three-dimensional regions,

FIG. 3 is a schematic diagram showing various functional components of a sound reproduction system.

FIG. 4 is a block diagram of a single input processor arrangement,

FIG. 5 is a block diagram of a multiple input processor arrangement,

FIG. 6 is a block diagram of auralisation processor arrangement,

FIGS. 7 a and 7b are perspective views of a microphone array,

FIG. 7 c is a detailed view of part of the microphone array of FIGS. 7a and 7b, and

FIG. 7 d is a plan view of microphone array of FIGS. 7a and 7b

DETAILED DESCRIPTION

The theoretical background on which the inventive apparatus is based is as follows. The signals driving the array of loudspeakers are designed to be such that the difference (more specifically the L² norm of the difference) between the desired or target sound field and the sound field generated by the array of loudspeakers is minimized on the surface of a three dimensional region, herein called the control volume. The problem is formulated mathematically as an integral equation of the first kind, and its solution is used in order to suitably configure the signal processing apparatus embedded in the system. The solution of the integral equation is often an approximated solution because of the mathematical ill-posedness of the problem. As described in detail below, the solution of the integral equation can be computed either with an analytical approach or with a numerical method. The two cases define two different approaches for the design of the signal processing apparatus. Both methods are described in this application. An aspect of the below described embodiments is that different choices of some parameters and/or solution methods are chosen depending on the frequency of the sound to be reproduced.

When the integral equation is solved with a numerical method, the latter requires that the target sound field is defined at a finite number of points on the boundary of the control volume, which is therefore sampled. The numerical solution of the integral equation is influenced by two sources of error: one due to ill-conditioning and one due to spatial aliasing. The former is more likely to affect the performance of the system at low frequencies, while the latter degrades the performance at high frequencies. These undesired effects can be avoided or contained by a wise choice of the control volume and of the sampling scheme used on the boundary of the control volume. It can be shown that at low frequencies the effects of ill-conditioning can be limited by choosing a large control volume. On the other hand the effects of the spatial aliasing can be reduced by choosing a regular sampling scheme of the surface of the control volume, for which the average distance between any two neighboring sampling points is less than half of the wavelength considered. For this reason, we define a control volume which is a function of the frequency of the sound to be reproduced: a large control volume at low frequencies and a gradually smaller control volume for higher frequencies. Even though the size of the control volume varies with the frequency its shape does not vary, and the number and geometrical arrangement of the samples of the boundary of that region generally does not vary, but in some cases it may be advantageous to allow such a variation. This approach allows the definition, for each frequency, of a control volume which is large enough to avoid the problems arising from ill-conditioning, keeping at the same time the distance between neighboring sampling points below the spatial aliasing limit.

A further aspect of the embodiments below relate to the solution of problems which can arise at certain frequencies that correspond to the so-called Dirichlet eigenvalues of the control volume. These problems arise when an attempt is made to reconstruct a sound field which is defined, at those frequencies, only on the boundary of the control volume. These critical frequencies are determined by the shape and size of the control volume. Our choice of a frequency dependent control volume has been applied to overcome these difficulties. In addition to that, we have chosen to define or measure the sound field at one or more locations in the interior of the control volume. In the design of the microphone array discussed below, when the microphones are arranged on spherical layers, we have chosen to include a microphone also in the centre of the microphone array. This choice has proven to overcome the problems due to the first critical frequency.

Another aspect of the embodiments below is related to the reproduction of a high frequency sound field due to a point source at a given location. In this case the approach used at low frequencies would lead all or almost all the loudspeakers of the array to contribute to the reproduced sound field. They would generate a large amount of acoustic energy, and the sound fields due to the different loudspeakers would generate a complex pattern of destructive and constructive interferences which would reconstruct the target field just over a limited region of the space. We have observed that the digital filters, which are part of the signal processing apparatus and correspond to those loudspeakers which are closer to the location of the virtual source, exhibit an asymptotic behavior at high frequencies. The amplitude of these filters tends to be constant and their phase tends to become linear. This fact has been proven mathematically in the special case when the loudspeakers are regularly arranged on the surface of a sphere, but similar conclusions can be inferred also for other geometric arrangements. For this reason, for the Single Input Mode described below, the processing applied to a high frequency component of the audio signal is different than that applied to the low frequency component of the signal. For the high frequency processing, only the loudspeakers which are the closest to the location of the virtual source are activated, and the corresponding digital filters have been substituted with simple gains, which resemble the asymptotic behavior of those filters. The Single Input Mode of the system may be viewed as a hybrid of a sound field reconstruction system and a three dimensional sound panning system. This approach has multiple advantages: in the first place, the amount of acoustic energy generated by the system is reduced. As a second instance, the real time calculation of the gains is computationally considerably less expansive than the real time calculation of an equivalent number of digital filters. Finally, the fact that only the loudspeaker, or those loudspeakers, close to the virtual source location are activated provides better cues at high frequencies for the human localization of a virtual sound source.

An especially innovative aspect of the embodiments described below is that the control volume can be chosen to be dependent on the frequency of the sound to be reproduced, and different processing steps are in general applied to the signals at different frequency bands. It is also worthwhile mentioning that the recording devices and reproduction devices of this system are designed to operate as a unique system, and there is no intermediate audio format on which the information between the recording part and the reproduction part of the system is transmitted. Finally, a further innovative feature of this invention is constituted by the Auralisation Mode (and the Auralisation Processing Unit), which allows to apply the theory of the reconstruction of a sound field to the design of a multi-channel reverberation simulator.

The principles described above have been applied to the design of the different components which constitute the signal processor apparatus of the sound reproduction system. The signal processor apparatus can be used to generate sound fields of three different characteristics:

• 1. A sound field generated in the free field by a virtual point source, whose location in the space and whose radiation pattern are selected.
• 2. A sound field generated by an omnidirectional point source in a reverberant environment, whose reverberant field is described by a set of measured or simulated impulse responses.
• 3. A generic sound field, described by a set of recorded sound signals acquired using a microphone array.

These three different characterizations define three different operating modes for the system, which are here defined as the Single Input Mode, the Auralisation Mode and the Multiple Input Mode, respectively. The general layout of the sound reproduction system comprises sound recording devices, modular signal processors and an array of loudspeakers (and their associated amplifiers). While the loudspeaker array is the same for the three different operational modes, the input devices and the signal processing are different for each of the modes.

In the Single Input Mode, the input signal is a single audio signal. The latter is obtained by capturing with a microphone the sound field generated by a sound source in an anechoic environment, or alternatively in a moderately reverberant environment, with the microphone located at a short distance from the source of sound. Additional data are provided, which describe the location of the virtual source (azimuth, elevation and distance) and preferably its radiation pattern as a function of two angles and of the frequency. The signal is digitally processed by the Single Channel Processing Unit, which is described in detail below. The location of the virtual source can be within or outside the loudspeaker array, but always outside of the control volume.

In the Multiple Input Mode, the input is a set of audio signals acquired with an array of microphones, designed for this purpose. The signals are processed by the Multiple Input Processing Unit, which is described in detail later and which is constituted by a matrix of digital filters. These are computed from a set of impulse responses describing the transfer function between each loudspeaker of the loudspeaker array and each microphone of the microphone array. These impulse responses can be either measured or computed from a theoretical model.

In the Auralisation mode the input signal is the same as for the Single Input Mode (a single anechoic audio signal), while the signal is digitally processed by the Auralisation Processing Unit described later. The latter is constituted by a set of digital filters, which are computed from two sets of impulse responses. The first set is constituted by the impulse responses of the considered reverberant environment measured with a specially designed microphone array. The second set of impulse responses is the same as that described for the Multiple Input Mode (describing the transfer functions between loudspeakers and microphones of the two arrays).

There is no a priori constraint on the loudspeaker arrangement, apart from the mild constraint that they should be arranged on a three dimensional surface, which surrounds the listening area and their axes are pointed towards the geometrical centre of the array. However, the system gives the best performance when a relatively large number of loudspeakers (eight or more) is used, they exhibit a preferably omnidirectional radiation pattern and they are regularly arranged on the surface of a sphere or of a hemisphere. Regular arrangement means here that the average distance between any two neighboring loudspeakers is constant or almost constant.

An example is now given to highlight the differences between the multiple input mode and the auralisation mode. One desires to reproduce the effect of the sound field generated by a violin playing in a concert hall (which is a reverberant environment). One can either

1. Use multiple input mode by recording the violin playing in the hall with the microphone array and then use the output of the microphone array as the input to the Multiple Input Processing Unit. The information on the sound of the violin and on the reverberant field of the hall are captured at the same time by the array and cannot, in principle, be divided, or

2 Use the auralisation mode by recording with the microphone array a set of impulse responses describing the reverberant field of the hall considered. These measurements do not depend on the sound of the violin and are used to calculate the digital filters in the Auralisation Processing Unit. Then one records the violin playing in an anechoic environment (an artificial non-reverberant environment, without reflections) and then use that single audio signal as the input to the Auralisation Processing Unit. The information on the sound of the violin and on the reverberant filed of the hall are captured separately.

In summary, the input to the multiple input unit are the signals from the microphone array. No assumption regarding the form of the field is made and the reverberant characteristic of the hall is not involved on the calculation of the digital filters. They are computed only from the knowledge of the characteristics of the microphone array and of the loudspeaker array geometry.

On the other hand, the input to the Auralisation Processing unit is a single audio signal (the sound of the violin playing in a non-reverberant environment). The reverberant characteristics of the hall considered, represented by a set of impulse responses, are this time involved in the calculation of the digital filters of the Auralisation Processing Unit.

While the reverberant field and the sound of the violin are captured together for the Multiple Input Mode, they are in the Auralisation Mode captured with separate procedures and are artificially merged by the auralisation processing unit.

The Multiple Input Mode has the advantage of capturing and reproducing a natural and real sound of an acoustic source in a given reverberant or non-reverberant environment.

The Auralisation mode has the advantage that the reverberant characteristic of the room and the sound from the given acoustic source are acquired separately and merged together later. With the Auralisation mode, it is possible to use the same reverberant hall and to change the sound source (for example a violin, a piano etc playing in a given concert hall) or conversely to use the same source of direct sound (the violin) and change the reverberant environment (with the artificial effect of having the same musician playing his/her violin in different concert halls).

FIGS. 7 a to 7d show microphone array 100 comprising a framework for supporting a set of thirty microphones. The array is designed to be used for the auralisation mode, for measuring the impulse responses of the reverberant environment to be emulated. The framework comprises two substructures 110 and 120 of substantially arcuate form. The substructure 110 comprises two spaced apart arcuate members 112 and 113 connected by a bridging member 114. The substructure 120 comprises two spaced apart members 122 and 123, the spaced apart members connected by a bridging member 124. Each of the substructures is disposed at different radii 121 and 122 from a centre point. The substructures 110 and 120 are pivotably mounted about a pole 150. Apertures 140 are provided in the arcuate members and the bridging members and are dimensioned to locate the microphone devices (not illustrated). A removable rigid sphere 160 is located inside the array. It is to be noticed that the two layers of microphones do not extend over a entire sphere but extend only over sectors. For this reason, the measurement of the impulse responses must be repeated eight times. At each measurement, thirty signals are acquired by the thirty microphones. After each measurement step, the array is rotated of 45° around its vertical axis. The eight measurements allow one to obtain a set of two hundred and forty impulses, captured at locations corresponding to an approximately regular sampling of two concentric spherical surfaces. Because the two hundred and forty signals can not be captured at the same time, this microphone configuration can be used for the auralisation mode only and not for the multiple input mode, for which the microphone signals should be acquired simultaneously.

Theoretical Analysis

In what follows, vectors are represented by lower case bold letters. The convention for the spherical co-ordinates r_x, θ_x, φ_x of a given vector x is illustrated in FIG. 1. Matrices are represented by capital bold letters. The imaginary unit is given by i=√{square root over (−1)}. The symbol [•]* represents the complex conjugate and the symbol [•]^H represents the Hermitian transpose (or conjugate transpose) of a matrix. δ_nm is the Kronecker delta (δ=1 if n=m, δ_nm=0 if n≠m). Let V⊂R³ be a bounded and simply connected volume of the three dimensional space, with boundary a ∂V (of class C²). A cross section of this volume is represented in FIG. 2. This region is referred to as the control volume. The assumption is made that the target sound field p(x,t) is due to sources of sound that are not contained in V, and that no scattering object is contained within this region. p(x,t) satisfies the homogeneous wave equation

$\begin{array}{cc}{\nabla }^{2}p\left(x,t\right)-\frac{1}{c^2}\frac{{\partial }^{2}p\left(x,t\right)}{\partial t^2}=0& \left(1\right)\end{array}$ in V, where c is the speed of sound, considered to be uniform in V. When a single frequency ω is considered and p(x,t)=Re{p(x)e^−iωt}, equation (1) can be reformulated as the Helmholtz equation∇²p(x)+k²p(x)=0xεV  (2)where k=ω/c is the wave number and the time dependence e^−iωt has been omitted. Let now A⊂R³ be a bounded and simply connected region of the space, with boundary ∂Λ of class C², that fully encloses V. Assume now that a continuous distribution of an infinite number of secondary, monopole-like sources is arranged on the boundary ∂Λ. This continuous distribution of secondary sources is the ideal model of the loudspeaker array, which is useful for the mathematical formulation of the problem. Later on the assumption that the number of secondary sources is infinite will be removed.

The assumption is made that the sound field p_y(x) generated by the secondary source located at yε∂Λ can be represented by a Green function G(x|y), thus satisfying the inhomogeneous Helmholtz equation∇²p_y(x)+k²p_y(x)=−a(y)δ(x−y)xε V, yε∂Λ  (3)where the function a(y) represents the complex strength of the secondary sources. In a free field, p_y(x) can be represented by the free field Green function

$\begin{array}{cc}{p}_y\left(x\right)=a\left(y\right)g\left(x|y\right)=a\left(y\right)\frac{e^{iky-x}}{4\pi y-x}x\in R^3,y\in \partial \Lambda & \left(4\right)\end{array}$

In a reverberant environment, the Green function has a more complex expression, which strongly depends on the shape of the reverberant enclosure and on its impedance boundary conditions.

The sound field {circumflex over (p)}(x) generated by the infinite number of secondary sources uniformly arranged on ∂Λ can be expressed as an integral, representing the linear superposition of the sound fields generated by the single sources:{circumflex over (p)}(x)=(Sa)(x)=∫_∂ΛG(x|y)a(y)dS(y)xεV  (5)

In what follows, {circumflex over (p)}(x) is also referred to as the reconstructed or reproduced sound field, while the integral introduced is often referred to as a single layer potential and a(y) is called the density of the potential.

Let now consider the following differential equation∇²p(x)+k²p(x)=0 xεV p(x)=ƒ(x) xε∂V  (6)where the function ƒ(x) describes the value field p(x) on the boundary ∂V. This differential equation is known as the Dirichlet problem (and the related boundary condition is called after the same name). As the Kichhoff-Helmholtz integral equation suggests, the knowledge of both the sound field and its normal derivative on the boundary (the Cauchy boundary condition) uniquely defines a sound field in the interior region V. However, this condition can be relaxed. In fact, under the above mentioned conditions and with the appropriate regularity assumption on the function ƒ(x), the Dirichlet problem (5) has a unique solution. This implies that the knowledge of the acoustic pressure on the boundary ∂V of the control volume is enough to define completely the sound field in the interior of volume V. This holds as long as the wave number k is not one of the Dirichlet eigenvalues k_n. The latter are defined as the infinite and countable wave numbers k_n such that the differential equation (5) with homogeneous boundary conditions ƒ(x)=0 is satisfied.

It is worth mentioning that the sound field in V can be ideally extended by analytical continuation to the exterior of the control volume, provided that the considered region does not contain any source of sound. This implies that if the acoustic field is reconstructed perfectly on the boundary of the control volume V, then it is reconstructed perfectly also in its interior and partially in the neighboring exterior region.

The determination of the density a(y) can be formulated asp(x)=(Sa)(x)=∫_∂ΛG(x|y)a(y)dS(y)xε∂V  (7)where p(x) is given and the density a(y) is the unknown of the problem. Note that the integral operator (Sa)(x) has been defined here as the restriction of the single layer potential (5) to the boundary ∂V. Equation (7) is an integral equation of the first kind and the determination of a(y) from the knowledge of p(x) on the boundary ∂V represents an inverse problem. It is important to highlight that equation (7) represents a problem that is, in general, ill-posed. This implies that a solution a(y) might not exist, and even if it exists it might be non-unique or not continuously dependent on the data p(x). The latter concept implies that small variations or errors on p(x) can result in very large errors in the solution a(y), which is therefore said to be unstable. It is however always possible to compute an approximate but robust solution by applying a regularization scheme.

It is important to highlight that the integral equation (7) is different from the Kirchhoff-Helmholtz equation (often also called Green formula) on which the technology called Wave Field Synthesis is grounded. The first main difference is that the integrand in the Kirchhoff-Helmholtz equation involves both monopole-like and dipole-like secondary sources, and the expression of their strength is expressed explicitly. On the other hand, the proposed approach relies on the use of monopole-like secondary sources only, and the determination of their strength is determined by the solution of the integral equation. The second main difference is that the field is known on the boundary ∂V, which is not the same as the boundary ∂Λ on which the secondary sources are arranged. This point describes also the main difference between the proposed approach and the Simple Source Formulation. It is possible to choose the control volume as a function of the frequency, while the secondary sources are arranged on the surface ∂Λ which does not depend on the frequency.

It is now possible to seek a solution to equation (7) using two different methods, one based on an analytical approach and the other based on a numerical solution of discretised version of the integral. As described above, the solution of the integral is applied to the design the signal processing unit embedded in the proposed system. In order to do that, the continuous solution provided by the analytical approach needs to be slightly modified in order to be adapted to the finite number of loudspeakers. The number of loudspeakers is L and the position of the acoustic centre of the l-th loudspeaker is identified by the vector y_l.

Analytical Solution of the Integral Equation

The first method to calculate a solution to equation (7) is an analytical method based on the singular value decomposition of the integral operator and is described in what follows.

To begin with, the scalar product of two square integrable functions ƒ(y) and g(y), with the same domain D, is here defined asƒ|g=∫_Dƒ(x)*g(x)dS(x)  (8)

The adjoint operator S⁺ of S is defined to be such thatg|Sa=S⁺g|a  (9)

It can be easily verified that, for the operator S introduced in equation (7), the adjoint operator S⁺ is given by(S⁺g)(y)=∫_∂VG(y|x)*g(x)dS(x)yε∂Λ  (10)

The physical meaning of the single layer potential (Sa)(x) is represented by the sound field generated by the continuous distribution of secondary sources on ∂Λ, evaluated at xε∂V. Similarly, its adjoint operator (S⁺g)(y) can be regarded as the time reversed version of the sound field generated by a continuous distribution of monopole-like secondary sources on ∂V, evaluated at yε∂Λ. The time reversal is due to the fact that the kernel of the integral (10) is the complex conjugate of the Green function G(y|x). Considering, as an example, the case of the free field Green function g(x|y) introduced in equation (4), it can be easily verified that g(y|x) and g(y|x)* differ because of the sign of the exponential. This implies that while g(•|x) can be considered as the representation of a sound field generated by a source of outgoing (diverging) spherical waves located at x, g(•|x)* could be understood as the case of a source of incoming (converging) spherical waves located at the same position x. Alternatively, the sound field represented by g(ω|x)* could be regarded as the time reversed version of g(•|x), as

$\begin{array}{cc}\mathrm{Re}\left\{e^{-\mathrm{i\omega }t}{g\left(x|y\right)}^{*}\right\}=\mathrm{Re}\left\{\frac{e^{-i\left(kx-y+\omega \left(-t\right)\right)}}{4\pi x-y}\right\}=\mathrm{Re}\left\{e^{-\mathrm{i\omega }\left(-t\right)}g\left(x|y\right)\right\}& \left(11\right)\end{array}$

For both interpretations, the generated spherical wave fronts are converging towards x. As a consequence of what has been said, the adjoint operator S⁺ could be interpreted as a continuous distribution of “sources of incoming waves” on ∂V. It can be shown that the operator S is compact and therefore its adjoint operator is compact too and the composite operator S⁺S is compact and self-adjoint. It is therefore possible to apply the properties of compact self adjoint operators and to perform a spectral decomposition of S⁺S. As a first step, the set of functions a_n (y) is computed, these functions being solutions of the eigenvalue problem(S⁺Sa_n)(y)=λ_na_n(y)yε∂Λ, n=1,2,3 . . . ∞  (12)where the eigenvalues λ_n are real, positive numbers. Considering the interpretation of S and S⁺ introduced above, the effect of the composite operator (S⁺Sa)(y) could be understood as follows: a sound field is generated by the continuous distribution of monopole-like sources on ∂Λ with strength a(y), and this field, on the boundary ∂V, is described by the function {circumflex over (p)}(x). Mathematically this implies that {circumflex over (p)}(x)=(Sa)(x). A different sound field is then generated by a continuous distribution of monopole-like sources on ∂V (the operator S⁺), and the strength of the sources is determined by the function {circumflex over (p)}(x). The sound field is then time reversed (or alternatively the secondary sources on ∂V can be replaced by “sources of incoming waves”), and the function â(y) describes the value of this field on ∂Λ. In mathematical terms, this corresponds to â(y)=(S⁺{circumflex over (p)})(y). Summarizing, the effect of the operator S⁺S can be understood as if the field generated by the secondary sources on ∂Λ was propagated from ∂Λ to ∂V and then propagated back to ∂Λ. Attention is now focused on the relation between â(y) and a(y): if these two functions are such that â(y)=λ_na(y), λ_nεR⁺, then the function a(y) is one of the eigenfunctions a_n(y) of S⁺S, and is a solution of (12). This means that the action of S⁺S on a_n(y) is simply an amplification or attenuation of the function, corresponding to positive real number λ_n.

The set of functions a_n(y) constitutes an orthogonal set of functions for ∂Λ meaning that any square integrable function a(y) defined on ∂Λ can be expressed as

$\begin{array}{cc}a\left(y\right)=\sum _{n=1}^N\langle a_n|a\rangle a_n\left(y\right)+\left(Qa\right)\left(y\right)& \left(13\right)\end{array}$

The integer N depends on the dimension of the range of S and might be infinite. (Qa)(y) is the orthogonal projection of a(y) on the null-space of S. The latter is defined as the set of functions ã(y) such thatN(S)={{tilde over (a)}(y):(Sã)(x)=0}  (14)

If this set is empty, then the set of functions a_n(y) is complete and equation (13) can be regarded as a generalized Fourier series. We can also generate a set of orthogonal functions p_n(x) on ∂V by letting the operator S act on the functions a_n (y) such thatσ_np_n(x)=(Sa_n) n=1,2,3, . . . ∞  (15)where the positive real numbers σ_n=√{square root over (λ_n)} are the singular values of S. It can be also proved thatσ_na_n(x)=(S⁺p_n) n=1,2,3, . . . ∞  (16)

It is possible to express any square integrable function p(x) on ∂V as the infinite series

$\begin{array}{cc}p\left(x\right)=\sum _{n=1}^N\langle p_n|p\rangle p_n\left(x\right)+\left(Rp\right)\left(x\right)& \left(17\right)\end{array}$ (Rp)(x) being the orthogonal projection of p(x) on the null-space of S⁺. Combining equations (13) and (15) and keeping in mind that (S(Qa))(x)=0 because of the definition of the null-space (14), it is possible to express the action of S on a(y) as

$\begin{array}{cc}\left(\mathrm{Sa}\right)\left(x\right)=\sum _{n=1}^N{\sigma }_n\langle a_n|a\rangle p_n\left(x\right)=\hat{p}\left(x\right)x\in \partial V& \left(18\right)\end{array}$

It is now possible to calculate an approximate solution of the integral equation (7) as

$\begin{array}{cc}a\left(y\right)=\sum _{n=1}^N\frac{1}{{\sigma }_n}\langle p_n|p\rangle a_n\left(y\right)& \left(19\right)\end{array}$

This equation provides a very powerful method for solving the problem of sound field reproduction with a continuous layer of monopole-like secondary sources. In fact, considering equations (15), (17) and (19), the single layer potential with the density computed in equation (18) is given by

$\begin{array}{cc}\begin{array}{c}{\int }_{\partial \Lambda }G\left(x|y\right)a\left(y\right)dS\left(y\right)={\int }_{\partial \Lambda }G\left(x|y\right)\left(\sum _{n=1}^N\frac{1}{{\sigma }_n}\langle p_n|p\rangle a_n\left(y\right)\right)dS\left(y\right)=\\ =\sum _{n=1}^N\frac{1}{{\sigma }_n}\langle p_n|p\rangle \left(S{a}_n\right)\left(x\right)\\ =\sum _{n=1}^N\langle p_n|p\rangle p_n\left(x\right)\\ =p\left(x\right)-\left(Rp\right)\left(x\right)\\ =\hat{p}\left(x\right)\end{array}& \left(20\right)\end{array}$

This implies that the reconstructed sound field {circumflex over (p)}(x) is the component of the target field that does not belong to the null-space of the adjoint operator S⁺, or equivalently that {circumflex over (p)}(x) is the projection of p(x) on the sub-space defined by the range of S. In other words, with the usual condition on the Dirichlet eigenvalues, if the target field has a pressure profile p(x) on ∂V that can be expressed as a linear combination of the orthogonal functions p_n(x), then it is ideally possible to determine a density a(y) such that {circumflex over (p)}(x)=p(x) in V. In the case when p(x) is not in the range of S, it can be shown that the reconstructed field {circumflex over (p)}(x) in (20) is the function that belongs to the range of S that minimizes the L² norm of the difference ∥p(x)−{circumflex over (p)}(x)∥_L2 on ∂V.

Considering equation (16), it can be noticed that even if some of the functions p_n(x) do not rigorously belong to the null-space of S⁺, their related singular value σ_n can be so small that it is possible to consider these p_n(x) as if (S⁺p_n)(x)≈0. This sheds some light on the ill-conditioning of the inverse problem represented by the integral equation (7). In fact, if for a given n the corresponding singular value σ_n is very small, then its reciprocal is very large and the norm of the density a(y) computed with the series (19) might become unreasonably large. Furthermore, the factor 1/σ_n is related to the amplification of errors contained in the data p(x) and therefore to the stability of the system. In order to prevent this ill-conditioning problem, it is possible to regularize the solution (19), applying a cut-off to the spectrum of the operator, using the Tikhonov regularization or any other regularization technique. In the design of the digital filters used in the system described in this application, a combination of spectral cut-off and Tikhonov regularization is used and the application of this is described later.

It has been shown that the possibility of calculating a density a(y) for reconstructing the target sound field p(x) with the single layer potential (5) from the knowledge of the boundary values of p(x) on ∂V is strongly related to the null-space of the adjoint operator S⁺. It is important to emphasize that if it is not possible to calculate exactly the density a(y), this does not imply that the latter does not exist and that the target sound field can not be perfectly reconstructed by the continuous distribution of secondary sources. An analogous if not identical problem arises in the field of Acoustic Holography: suppose that an attempt is made to determine the surface normal velocity of a radiating plate from the knowledge of the radiated sound field on a given region. It turns out that the vibration modes that generate evanescent waves cannot in practice be determined from far field measurements. However, this does not mean that these low-efficiency vibro-acoustic modes do not exist. This highlights the fact that the feasibility of the reconstruction of the target sound field not only depends on the arrangement ∂Λ of the secondary sources, but also is very much related to the surface ∂V on which the sound field has been defined or measured.

This method of solution has the big advantage that the density a(y) has an analytical expression and the design of the digital filters implemented in the system does not require any numerical matrix inversion. However, there are two disadvantages of this method. The first is that the eigenfunctions a_n (y) strongly depend on the geometry of V and Λ and their explicit calculation is usually not trivial. Their formulation is known for a limited number of geometries. The second disadvantage is that, when the continuous distribution of secondary sources is substituted by an array of a finite number of loudspeakers, the performance of the system whose signal processing units have been designed with this method are more effective if the distribution of the secondary sources (the loudspeakers) is regular.

As an example, the special case is now considered, where the boundaries of the two volumes V and Λ are two concentric spheres, with radius R_V and R_Λ, respectively. The assumption of a free field is also made here. For the case under consideration the functions a_n(y) and p_n(x) and the singular values σ_n can be expressed analytically as

$\begin{array}{cc}{a}_n\left(y\right)=\frac{1}{R_{\Lambda }}{Y}_n^m\left({\theta }_y,{\varphi }_y\right)p_n\left(x\right)=\frac{{\gamma }_n}{R_V}{Y}_n^m\left({\theta }_x,{\varphi }_x\right){\sigma }_n=R_{\Lambda }{R}_{V}kh_n^{\left(1\right)}\left(k{R}_{\Lambda }\right)j_n\left(k{R}_V\right)& \left(21\right)\end{array}$ j_n(•) is the spherical Bessel function of order n, h_n^(l)(•) is the spherical Hankel function of the first kind and order n. The spherical harmonics Y_n^m(θ,φ) are defined as

$\begin{array}{cc}{Y}_n^m\left(\theta ,\varphi \right)=\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4\pi \left(n+m\right)!}}{P}_n^m\left(\cos \theta \right)e^{im\varphi }& \left(22\right)\end{array}$ where P_n^m(•) are associated Legendre functions. The factor γ_n, having unitary norm, is given by

$\begin{array}{cc}{\gamma }_n=i\frac{h_n^{\left(1\right)}\left(k{R}_{\Lambda }\right)j_n\left(k{R}_V\right)}{h_n^{\left(1\right)}\left(k{R}_{\Lambda }\right)j_n\left(k{R}_V\right)}& \left(23\right)\end{array}$ and represents a phase shift applied to each spherical harmonic of order n due to the action of S. Equation (21) can be verified by substituting it into equation (12), (15) or (16) and applying the spherical harmonic expansion of the free field Green function (A1) together with the completeness and orthogonality relations of the spherical harmonics, equations (A5) and (A6) respectively.

It can be noticed that the spherical harmonics Y_n^m(θ,φ) have two indices, while the singular values σ_n have only one index. This is due to the degeneracy of the singular values, and it implies that one eigenspace of dimension (2n+1) is associated with the singular values σ_n. Hence, for each order n, it is possible to generate a set of (2n+1) orthogonal spherical harmonics which span that subspace. In other words, all the spherical harmonics of order n and degree m are associated with the same singular values σ_n. This degeneracy is typical of symmetrical geometries (such as the sphere), and arises in many other fields of physics (a well known example in quantum physics is the degeneracy of two electronic configurations, which have the same energy level).

As a result of the orthogonality of the spherical harmonics (A5), it is easy to verify the mutual orthogonality of the set of functions p_n(x) and a_n(y) defined by equation (21). Using the expansion of the free field Green function (A1), it possible to show that the functions a_n(y) and p_n(x) and the singular values σ_n satisfy equations (12), (15) and (16).

The spherical wave spectrum S_nm(r) of the target sound field p(x), calculated at r=R_V, is defined asS_nm(R_V)=∫₀^2πdφ∫₀^πp(R_V,θ_x,φ_x)Y_n^m(θ,φ)*sin(θ)dθ  (24)

It is now possible to calculate analytically the density a(y) as

$\begin{array}{cc}a\left(y\right)=\sum _{n=0}^{\infty }\sum _{m=-n}^n\frac{Y_n^m\left({\theta }_y,{\varphi }_y\right)}{R_{\Lambda }^{2}ik{h}_n^{\left(1\right)}\left(k{R}_{\Lambda }\right)j_n\left(k{R}_V\right)}{S}_{\mathrm{nm}}\left(R_V\right)& \left(25\right)\end{array}$

It can be observed that the denominator of (22) equals zero for those wave numbers k_n such that j_n(k_nR_V)=0. These wave numbers correspond to the Dirichlet eigenvalues introduced previously, and identify the frequencies of resonance of a spherical cavity with radius R_V and pressure release boundaries (sometimes called the Dirichlet sphere). Under these circumstances, it is not possible to calculate uniquely the density a(y) with (25) because of the non uniqueness of the solution of the Dirichlet problem (6). Considering equation (16), it is easy to notice that the function Y_n^m(θ_x,φ_x) belongs to the null-space of S⁺ and therefore is not in the range of the single layer potential S. Once again, it is important to emphasize that this problem in the reconstruction is not due to the arrangement of the layer of secondary sources on ∂Λ, but it is due to the boundary ∂V where the target sound field has been defined. For that reason, the density a(y) can be determined by changing the radius of the control volume V.

The determination of the spherical spectrum of the target sound field involves the calculation of the integral (24). For some special cases, this integral can be solved analytically. Considering the expansion of the free field Green function (A1) it is possible to derive the following relation for the spherical spectrum of an omnidirectional point source located at z>R_V (formally a monopole with volume velocity q_VOL=−(iρck)⁻¹):S_nm^ps(R_V)=ikh_n¹)(kr_z)j_n(kR_V)Y_n^m(θ_z,φ_z)*  (26)

It important to emphasize that the source location z should not be in V but can be within Λ. Considering the spherical harmonic summation formula (A3) and the trigonometric relations (A4), a simplified formula can be derived for the density for the reconstruction of the sound field due to a monopole source:

$\begin{array}{cc}{a}^{\mathrm{ps}}\left(y\right)=\sum _{n=0}^{\infty }\frac{h_n^{\left(1\right)}\left(k{r}_z\right)}{R_{\Lambda }^2{h}_n^{\left(1\right)}\left(k{R}_{\Lambda }\right)}\frac{2n+1}{4\pi }{P}_n\left(\cos \left(\zeta \right)\right)& \left(27\right)\end{array}$ where P_n(•) is the Legendre polynomial of degree n and (see relation (A4))

$\begin{array}{cc}\cos \left(\zeta \right)=\frac{y·z}{R_{\Lambda }{r}_z}=\sin \left({\theta }_y\right)\sin \left({\theta }_z\right)\cos \left({\varphi }_y-{\varphi }_z\right)+\cos \left({\theta }_y\right)\cos \left({\theta }_z\right)& \left(28\right)\end{array}$ and represents the cosine of the angle between the vectors y and z. It can be noticed that the summation over the different degrees m has been reduced to the computation of a single Legendre polynomial.

The application of this solution to the design of the signal processing apparatus of the system requires the series (25) and (27) to be truncated to a given order N. Under this assumption, for high operating frequencies such that kR_Λ, kr_z>>N, then the large argument limits or far field approximation (A7) of the spherical Hankel functions and the finite summation formula for Legendre polynomial (A2) can be used in (27), which can be rewritten as

$\begin{array}{cc}{a}^{\mathrm{hf}}\left(y\right)=\frac{e^{ik\left(r_z-R_{\Lambda }\right)}\left(N+1\right)}{R_{\Lambda }{r}_{z}4\pi }\frac{P_N\left(\cos \left(\zeta \right)\right)-P_{N+1}\left(\cos \left(\zeta \right)\right)}{1-\cos \left(\zeta \right)}& \left(29\right)\end{array}$

This is a very powerful formula, which is extremely useful for the real time implementation of a panning function. In fact, the only frequency dependent part of (29) is the complex exponential e^ik(rz−RΛ), and its Fourier transform corresponds in the time domain to a simple delay equal to the distance r_z−R_Λ. The denominator of (29) is a simple attenuation due to the distance of the point source and the term represented by the series of Legendre polynomials is a gain factor depending on the relative angle between z and y.

In the very special case when the distance of the virtual source r_z equals the sphere radius R_Λ, equation (27) can be rewritten in the very simple formulation

$\begin{array}{cc}{a}^{\mathrm{real}}\left(y\right)=\frac{\left(N+1\right)}{R_{\Lambda }^{2}4\pi }\frac{P_N\left(\cos \left(\zeta \right)\right)-P_{N+1}\left(\cos \left(\zeta \right)\right)}{1-\cos \left(\zeta \right)}& \left(30\right)\end{array}$

The continuous density function calculated with equations (19), (25), (27), (29) or (30) has to be transformed into a finite set of (possibly frequency dependent) coefficients corresponding to the each loudspeaker of the system. This can be done by applying a quadrature of the integral (7). The coefficient a_l corresponding to the l-th loudspeaker is therefore computed asa_l=a(y_l)ΔS_l  (31)where ΔS_l has the dimension of an area and depends on the loudspeaker arrangement. If these are arranged regularly, then

$\begin{array}{cc}\Delta S_l=\frac{A_{\partial \Lambda }}{L}& \left(32\right)\end{array}$ where A_∂Λ is the area of the boundary ∂Λ and L is the total number of loudspeakers of the system. If the boundary ∂Λ is a sphere of radius R_Λ, then

$\begin{array}{cc}\Delta S_l=\frac{4\pi R_{\Lambda }^2}{L}& \left(33\right)\end{array}$ Numerical Solution of the Integral Equation

The second method to solve the integral equation (7) is numerical. As a first step, the boundaries ∂Λ and ∂V must be sampled. The sampling scheme adopted for ∂Λ is given by the loudspeaker array: the boundary ∂Λ is divided into L surfaces, each of them corresponding to a loudspeaker. The surface ΔS_l corresponding to the l-th loudspeaker is the same as in equations (31), (32) or (33). The boundary ∂V is divided into Q surfaces. The q-th surface is identified by its geometrical centre x_q, hereafter called a sampling point. It is recommended that the number of sampling points is chosen to be such that Q>L. The sampling points should be chosen in such a way that the average distance δx between two neighboring points is constant or approximately constant. In order to avoid problems arising from spatial aliasing, it is recommended that, for a given angular frequency ω,

$\begin{array}{cc}\delta x<\frac{c\pi }{\omega }& \left(34\right)\end{array}$

In what follows, this relation is referred to as the aliasing condition. As will be discussed later, it might be desirable to choose a control volume V which is a function of the frequency.

In order to avoid the problems arising from the Dirichlet eigenvalues explained above, it can be useful to add some additional sampling points in the interior of V, thus increasing the number Q of sampling points. In the case of ∂V being a sphere, it might be a wise choice to add an additional sampling point in the centre of the sphere, in order to avoid the problems arising from the first Dirichlet eigenvalue, whose characteristic wave number is identified by the first zero of the spherical Bessel function j₀(kR_V). Another way of avoiding these problems is by choosing, for a given frequency {circumflex over (ω)}, a control volume V({circumflex over (ω)}) for which the frequency considered does not correspond to one of its Dirichlet eigenvalues. Some possible strategies for choosing a frequency dependent control volume are described later.

The vector p is defined as the set of values of the target field p(x) evaluated at the positions of the sampling points. In the case that the target field is due to an omnidirectional point source located at z in the free field, as in the case of equation (26), then the q-th element of the vector p is defined as

$\begin{array}{cc}{p}_q=g\left(z|x_q\right)=\frac{e^{i{\mathrm{kd}}_{\mathrm{qz}}}}{4\pi d_{\mathrm{qz}}}& \left(35\right)\end{array}$ where d_qz=∥z−x_q∥ is the distance between the q-th sampling point on ∂V and the virtual source, and it can depend on the frequency if the control volume V also depends on the frequency. The dimension of p is Q. The operator S can be transformed into matrix s, which is defined as

$\begin{array}{cc}{S}_{\mathrm{ql}}=g\left(y_l|x_q\right)=\frac{e^{i{\mathrm{kd}}_{\mathrm{ql}}}}{4\pi d_{\mathrm{ql}}}& \left(36\right)\end{array}$ where d_ql=∥y_l−x_q∥ is the distance between the q-th sampling point on ∂V and the position of the l-th loudspeaker. This distance can depend on the frequency if the control volume V also depends on the frequency. The dimension of S is therefore Q by L. The regularized pseudo inverse matrix S⁺ (not to be confused with the adjoint operator) is defined asS⁺=(S^HS+βI)⁻¹S^H  (37)where β is a regularization parameter and I is the identity matrix of dimension L by L. The dimension of S⁺ is L by Q. It is important to emphasise that this matrix depends only on the loudspeaker arrangement and on the sampling scheme on the boundary of the control volume, and does not depend on the position of the virtual source. It is now possible to compute the coefficient a_l corresponding to the l-th loudspeaker as

$\begin{array}{cc}{a}_l=\sum _{q=1}^Q{S}_{\mathrm{lq}}^{+}{p}_q& \left(38\right)\end{array}$

It is important to notice that both matrices S and S⁺ and the vector p depend on the wave vector k and hence on the operating frequency ω.

A major feature of the two methods which have been presented for solving the integral equation (7) and to derive the loudspeaker coefficients a_l, is that they give identical results when the number of loudspeakers L tends to infinity.

It is now important to discuss how the control volume V can be modified depending on the operating frequency. For what has been said about the aliasing condition, it may seem to be wise to choose the control volume to be as small as possible. However, the study of the stability of the condition number of matrix S as a function of the operating frequency ω shows that if the control volume is too small, then the conditioning of the matrix S is poor.

In order to respect the sampling condition for all the considered frequencies and to have, at the same time, a well-conditioned matrix S, it might be desirable to choose a frequency dependent control volume V(ω). Suppose that, for a given frequency {circumflex over (ω)}, a control volume V({circumflex over (ω)}) has been chosen which is a star-convex set and a set of sampling points x_l({circumflex over (ω)}) have been chosen which respect the sampling condition, which grants good conditioning of the matrix S and {circumflex over (ω)} does not correspond to one of the Dirichlet eigenvalues for V({circumflex over (ω)}). In this case it is possible to define a frequency dependent control volume V(ω), which has the same shape of V({circumflex over (ω)}) and which is identified by the sampling points

$\begin{array}{cc}{x}_l\left(\omega \right)=\frac{\omega }{\hat{\omega }}{x}_l\left(\hat{\omega }\right)& \left(39\right)\end{array}$

Provided that Q>L, a suitable choice for a frequency dependent control volume V(ω) is a sphere centered in the origin with radius

$\begin{array}{cc}{R}_V=\frac{c\left(\sqrt{L}-1\right)}{\omega }& \left(40\right)\end{array}$

The motivation of this choice is non trivial and is related to the behavior of the spherical Bessel functions. It has been shown that for spherical geometries the Bessel functions appear in the expression (21) of the singular values σ_n, and the radius R_V of the control volume appears in the argument of these functions. It can be inferred that, even in the discrete case, the L eigenvalues of the squared matrix S^HS are related to the spherical Bessel functions. The degeneracy of the eigenvalues σ_n in the case of spherical geometries has also been discussed above. This implies that if S^HS has dimension L by L and always in the case of spherical geometry, then it has ideally N=√{square root over (L)}−1 independent eigenvalues. This approximation is exact when L tends to infinity and the boundaries ∂V and ∂Λ are sampled regularly. The condition number of S^HS is related to the ratio between the largest and the smallest singular values, both related to the spherical Bessel functions. Using these arguments, it can be deduced that an optimal choice for the radius R_V is given by R_V=N/k, which leads to equation (40).

In both cases of equation (39) and equation (40), the additional constraint must be applied that all loudspeakers and the virtual source position must be located outside of the control volume.

It is worth highlighting, once again, that the control volume does not correspond to the listening area, as an accurate reproduction of the target sound field can be achieved, because of analytical continuation, also in the exterior of the control volume.

If the target sound field is due to a point source located at z, it is useful to introduce the factorΔΦ_lz=e^{ik(rl−rz)+iωδt}  (41)

This represents a simple phase compensation factor, which compensates for the delay arising in the filter computation caused by the difference of the radial co-ordinate of the l-th loudspeaker r_l and of the virtual source r_z. The term δt is a small constant quantity, corresponding to a modeling delay, which has been introduced in order to guarantee the causality of the digital filters.

DESCRIPTION OF POSSIBLE IMPLEMENTATIONS

FIG. 3 reports a diagram of the signal processing apparatus of the sound reproduction system. The apparatus comprises functional modules, called Single Input Processing Units (SIPU), Multiple Input Processing Units (MIPU) and Auralisation Processing Units (APU). They are the signal processing units corresponding to the three different operational modes of the system described above. It will be appreciated that the limit on the total number of SIPU, MIPU and APU modules of the system is not given a priori and depends on the computational power of the electronic hardware used for the implementation of the system. The different modules are composed by different components, described in detail in what follows. Digital filters are one of these components. They are defined in what follows with mathematical expressions, and from these expressions digital Finite Impulse Response filters (or possibly Infinite Impulse Response Filters) can be easily designed using standard techniques, such as for example the frequency sampling method. The use of Finite Impulse Response filters is suggested, as the filters described in what follows usually show a strong decay in the time domain. It is also relevant to highlight that most of the components of the different modules, such as digital delays and filters, are often described as separate elements for sake of clarity. However, it can be useful implementing a series of one of more of these elements as a single digital filter. As an example, referring to FIG. 4, it is possible to implement the Low Pass Filter (LPF) and each of the filters in the low frequencies bus (F₁, F₂, etc.) as a single filter.

Single Input Processing Unit (SIPU)

FIG. 4 shows a block diagram of a Single Input Processing Unit. This module is designed to generate the L loudspeaker signals which allow the reproduction of the sound field due to a virtual source at a given position in the free field, with a given radiation pattern and with a given orientation. The SIPU receives as an input a single audio signal and some additional data, describing the location (distance and direction) of the virtual source, its radiation pattern and its orientation. With reference to FIG. 1, the distance of the virtual source r_x is described by a scalar and positive number, while the direction of the virtual source is described by the two angles θ_x and φ_x. The orientation of the virtual source is described by two angles and its radiation pattern by a complex function of the frequency and of two angles.

The audio signal is first filtered by the digital filter F_H, which is defined depending on the orientation and distance of the virtual source. The signal is then divided into two busses, called a high frequency bus and a low frequency bus respectively. A high pass filter HPF(ω) and a low pass filter LPF(ω) are applied to the two signals respectively. The two filters are such thatLPF(ω)+HPF(ω)=e^iωΔt

This means that if one signal is filtered separately by the two filters and the two outputs are summed together, the result is a delayed version of the original signal. The cut-on frequency of the high pass filter (−6 dB) and the cut-off frequency of the low pass filter (−6 dB) is the same and is called ω_c. The latter can be chosen to be the smallest frequency ω satisfying the condition

$\frac{h_N^{\left(1\right)}\left(r_{\mathrm{lmin}}\omega /c\right)}{h_0^{\left(1\right)}\left(r_{\mathrm{lmin}}\omega /c\right)}\approx 1$ where r_lmin is the smallest of all loudspeaker radial co-ordinates r_l, that is to say the radial co-ordinate of the loudspeaker whose position is the closest to the origin.

The signal on the high frequency bus and the signal on the low frequency bus are processed in different ways, applying a set of operation called High Frequency Signal Processing and Low Frequency Signal Processing, respectively. They are described in detail below.

It should be noted that the existence of both the High Frequency Bus and the Low Frequency Bus is not essential. For example, a variant of the Single Input Processing Unit may comprise either the High Frequency Signal Processing or the Low Frequency Signal Processing. In these special cases, it is clear that the High Pass Filter HPF(ω) and the Low Pass Filter LPF(ω) would need to be removed.

The signal on the low frequency bus is filtered in parallel by a bank of L digital filters, labeled F₁, F₂, . . . , F_L in FIG. 4. Each filter corresponds to a different loudspeaker. These filters can be defined in two different ways, called numerical filter computation and analytical filter computation, both depending on the position of the virtual source and on the position of the loudspeaker considered.

1—For the numerical filter computation, a control volume V is chosen as explained above. Its geometrical centre coincides with the origin of the co-ordinate system. As described above, a set of Q regularly arranged sampling points is defined on the control surface. The q-th sampling point is identified by the vector x_q. All loudspeakers and the location of the virtual source lie outside of the control volume. As has been discussed, the control volume and the sampling points can be chosen to be frequency dependent. A suitable choice for a frequency dependent control volume is a sphere with radius

$R_V=\min \left[\frac{c\left(\sqrt{L}-1\right)}{\omega },r_z,r_{\mathrm{lmin}}\right]$

As discussed above, it is beneficial to include some extra points in the interior of V in order to avoid problems arising from the Dirichlet eigenvalues. The frequency dependent vector p(ω) is defined as in equation (35). The frequency dependent matrices S(ω) and S⁺(ω) are defined as in equations (36) and (37). It is important to highlight that these matrices depend only on the loudspeaker arrangement and on the sampling scheme of the boundary of the reconstruction area, and do not depend on the position of the virtual source. For this reason, they can be computed off-line and not in real time. The digital filter corresponding to the l-th loudspeaker is computed from

$F_l\left(\omega \right)={\mathrm{\Delta \Phi }}_{\mathrm{lz}}\left(\omega \right)\sum _{q=1}^Q{S}_{\mathrm{lq}}^{+}\left(\omega \right)p_q\left(\omega \right)$

2—The digital filters in FIG. 4, when they are computed following the analytical approach, have a frequency dependence of the form

$F_l\left(\omega \right)=\Delta S_l\left(\omega \right){\mathrm{\Delta \Phi }}_{\mathrm{lz}}\left(\omega \right)\sum _{n=1}^L\frac{a_n\left(y_l,\omega \right)}{{\sigma }_n\left(\omega \right)+\beta \left(\omega \right)}{\int }_{\partial V\left(\omega \right)}g\left(z|x,\omega \right)p_n\left(x,\omega \right)*dS\left(x,\omega \right)$ where the functions a_n(y,ω) and p_n(x,ω) and the singular values σ_n(ω), defined by equations (12) and (15), are computed analytically as described previously. ∂V(ω) is the boundary of the control volume and β(ω) is a regularization parameter; both can be chosen to be frequency dependent. It is recommended to choose the order N of truncation of the series to be equal to (or possibly smaller than) the number of loudspeakers L. In the case when the loudspeakers are regularly arranged on a sphere with radius R_Λ and the control volume is also a sphere with radius R_V, then following equation (27) the frequency response of the digital filter corresponding to the l-th loudspeaker is defined by

$F_l\left(\omega \right)={\mathrm{\Delta \Phi }}_{\mathrm{lz}}\left(\omega \right)\sum _{n=0}^M\frac{h_n^{\left(1\right)}\left({\mathrm{kr}}_z\right)}{h_n^{\left(1\right)}\left({\mathrm{kR}}_{\Lambda }\right)}\frac{2n+1}{L}{P}_n\left(\cos \left(\zeta \right)\right)$ where cos(ζ) is defined by equation (28) and the order M of truncation of the series must be a natural number and is chosen to beM≦√{square root over (L)}−1

The signal on the high frequency bus is first delayed by an amount of time Δt that takes in consideration the length of the digital filters in the low frequencies bus and the quantity δt introduced in equation (41). The signal is then multiplied by a bank of parallel gains, each of them corresponding to a different loudspeaker. They are labeled G₁, G₂, . . . , G_L in FIG. 4. There are two possible ways of defining these gains. The first method is derived from equation (29), corresponding to the high frequency approximation of equation (27). The gain corresponding to the l-th loudspeaker is defined as

$G_l=\{\begin{array}{cc}\frac{r_l\left(N+1\right)}{r_z}\frac{P_M\left(\cos \left({\zeta }_l\right)\right)-P_{M+1}\left(\cos \left({\zeta }_l\right)\right)}{1-\cos \left({\zeta }_l\right)}\frac{\Delta S_l}{r_l^2}& \mathrm{if}{\zeta }_l<{\zeta }_M\\ 0& \mathrm{if}{\zeta }_l\ge {\zeta }_M\end{array}$ cos(ζ) is defined by equation (28). M depends upon the average distance of the neighboring loudspeakers and for a regular arrangement can be chose as M≦√{square root over (L)}−1. In the case when the loudspeakers are regularly arranged on the surface of a sphere, then ΔS_l/r_l²=1/L. The angle ζ_M defines the semi-aperture of the main lobe of the function (P_M(cos(ζ)−P_M+1(cos(ζ)))/(1−cos(ζ) and can be defined by the relationP_M(cos(ζ_M))+P_M+1(cos(ζ_M))=min[P_M(cos(ζ))+P_M+1(cos(ζ))]

This implies that only the loudspeakers which are closest to the location of the virtual source are active, and that they all operate in-phase.

The second method is derived from the assumption that the digital filters on the low frequency bus designed with the numerical approach and corresponding to the loudspeakers which are closer to the location of the virtual source show an asymptotic behavior at high frequencies. It is supposed that after the cut-off frequency ω_c the magnitude of these filters remains constant and the phase is very close to 0. The gain corresponding to the l-th loudspeaker is therefore defined as

$G_l=\{\begin{array}{cc}\sum _{k=1}^K{S}_{\mathrm{lk}}^{+}\left({\omega }_c\right)p_k\left({\omega }_c\right)& \mathrm{if}{\zeta }_l<{\zeta }_M\\ 0& \mathrm{if}{\zeta }_l\ge {\zeta }_M\end{array}$ where the angle ζ_M, the matrix S⁺(ω) and the vector p(ω) are defined as above.Microphone Array and Multiple Input Processing Unit (MIPU)

The Multiple Input Processing Unit is designed to generate the L loudspeaker signals which allow the reproduction of a sound field which has been captured using a specially designed array of microphones.

The array of microphones is designed in connection with the reproduction system, meaning that the microphone array is to be considered as a part of the whole system. The microphone array comprises a plurality of omnidirectional capsules regularly arranged on multiple surfaces. It will be appreciated that a microphone capsule relates to a portion where a microphone membrane is located. These surfaces define the boundaries of multiple, concentric control volumes. The choice of multiple control volumes arise from the fact that for a given number of sampling points the ideal size of the control volume, which respects the aliasing condition and allow a good conditioning of matrix S(ω), depends on the considered frequency. It is not practicable to choose a control volume which changes continuously as a function of the frequency. It is however possible to choose a finite number of control volumes V₁, V₂ . . . V_F, each of them dedicated to a given frequency range. A set of Q_ƒ omnidirectional microphones are regularly arranged on the boundary ∂V_ƒ of the control volume V_ƒ. The set of all the microphones arranged on the same control volume is referred to as a microphone layer. The total number of microphones Q is given by

$Q=\sum _{f=1}^F{Q}_f$

As explained above, problems can arise at those frequencies, which correspond to the Dirichlet eigenvalues of the control volume. The use of a multiple microphone layer can partially overcome this problem, but it can happen that one of the frequencies corresponding to one of the Dirichlet eigenvalues of the volume V_ƒ belongs to the range of frequencies to which that control volume is dedicated. This problem can be partially, when not totally, overcome by adding one or more microphones in the interior of the control volume. A wise choice for an additional microphone location is in the centre of the microphone array, especially when the control volumes are spherical. This is due to the fact that the first critical frequency for a given volume V_ƒ is identified by the first zero of the spherical Bessel function j₀(R_ƒω/c). In physical terms this means that the microphone array can not detect the component of the sound field corresponding to the zero order spherical harmonic. The microphone in the centre of the array can, on the contrary, detect only that missing component, thus overcoming the problem. A higher number of additional microphones might be needed for higher critical frequencies. It is also possible to use, as an additional sampling point for a given layer ∂V_ƒ, one of the microphones arranged on a different layer ƒ′≠ƒ.

A suitable choice for the different control volumes is given by a set of concentric spheres. As a guideline for the choice of the radius of these spheres, if the upper frequency for the operating frequency range of the control volume V_ƒ is ω_ƒ, the radius R_ƒ of that control volume can be chosen to be

$R_f=\frac{c\left(\sqrt{{\min }_{\mathrm{LQf}}}-1\right)}{{\omega }_f}$ where min_LQƒ is the smallest number between the number of loudspeakers L and the number of Q_ƒ of microphones on that layer. This guideline arises from considerations of the spherical Bessel functions and on the conditioning of matrix S^HS, which have been discussed above (refer to equation (40)). This approach suggests that it is also possible to split a given frequency range corresponding to a control volume V_ƒ with Q_ƒ microphones into two sub-ranges by simply reducing the number of microphones used in the processing of the lower frequency range. This can be especially useful at very low frequencies, where the choice of radius discussed above would result in a very large value. As an example, consider a system composed by forty loudspeakers and a layer of thirty six microphones regularly arranged on the sphere ∂V_ƒ dedicated to the audio frequencies below 500 Hz. Following the above guideline, the radius R_ƒ is approximately 0.5 m. It is possible to choose a subset of eight microphones on that layer and define an additional frequency range with higher limit of approximately 200 Hz.

A microphone array with just one layer can be considered as a special case of the microphone array described. Another variant is constituted by a microphone array having a scattering object (as for example a rigid sphere) in the region of the space contained inside the smallest control volume. The filter computation described in what follows remains the same. It is also straightforward to perform the analytical calculation of the digital filters described later for the case corresponding to a set of microphones arranged around or on the surface of a rigid sphere.

The output signals from the microphone array are processed by the Multiple Input Processing Unit, represented by FIG. 5. This figure illustrates the example corresponding to a microphone array with only two layers, but this approach can be identically extended to configurations with more microphone layers. The output signals of each layer are processed separately, as shown in the figure. As a first step, each signal is filtered through a Band Pass Filter (BPF_ƒ(ω)) whose cut-on and cut-off frequencies are defined by the frequency range corresponding to that microphone layer. For the layer ƒ, the cut-on frequency is ω_ƒ−1 and the cut-off frequency is ω_ƒ. These band pass filters must respect the condition

$\sum _{f=1}^{F}BP{F}_f\left(\omega \right)=e^{\mathrm{i\omega \Delta }t}$ analogous to the Low and High Pass Filters in the SIPU. The signals are then process by a matrix of digital filters, labeled F_1,1,1, F_L,1,1, . . . , F_L,Q11, in FIG. 5. These filters can be defined in three different ways: analytically, numerically or with measurements. All these formulation are grounded on the same theoretical background, which has been discussed above.

1—Using a purely numerical approach, the filters corresponds to the elements of the frequency dependent matrix S⁺(ω) defined by equation (37). The filter corresponding to the l-th loudspeaker and the q-th microphone on the layer ƒ is computed fromF_l,q,ƒ(ω)=e^iωδtS⁺_lq(ω)where the small δt has been introduced in order to ensure that the filter is causal (when this is needed).

2.—The filters can be also calculated after having measured, possibly in an anechoic environment, the impulse response between each loudspeaker and each microphone on the given layer. The measurement can be carried out using standard techniques (swept sine, MLS, etc.) and is subject to the well-known sources of error that affect these kinds of measurements. The microphone array must be arranged in such a way that its geometrical centre corresponds to the origin of the co-ordinate system. It is preferable to exclude reflections generated by the surrounding environment in the measurement. This can be done by carrying out the measurements in an anechoic environment or by windowing the measured impulse response in order to take into account only the initial part of the measured signal. The acquired impulse responses need to be transformed in the frequency domain by applying a Fourier transform. The set of acquired measurements constitutes the matrix H(ω). Its element H_qlƒ(ω) represents the transfer function between the l-th loudspeaker and the q-th microphone on the layer ƒ. Following a procedure analogous to equation (37), matrix H⁺(ω) is computed fromH⁺(ω)=(H^H(ω)H(ω)+β(ω)I)⁻¹H^H(ω)where the elements β(ω) and I are defined as for equation (37). It is very important to apply a regularization scheme to the inversion of matrix H^H(ω)H(ω), as the presence of measurement errors can result in the computation of unstable filters. In the proposed approach, Tikhonov regularization with frequency dependent regularization parameter β(ω) is suggested. The filter corresponding to the l-th loudspeaker and the q-th microphone on the layer ƒ is computed fromF_l,q,ƒ(ω)=e^iωδtH⁺_lqƒ(ω)where δt is as defined above.

3—An analytical computation method for the filters can be derived from discretising, in equation (19), the integral p_n|p defined by equation (8). Equation (19) can therefore be reformulated as

$a_f\left(y\right)=\sum _{q=1}^{Q}p\left(x_q\right)\Delta S_q^{\prime }\sum _{n=1}^M\frac{1}{{\sigma }_{n,f}}{p}_{n,f}\left(x_q\right)a_{n,f}\left(y\right)$ where ΔS′_q, analogous to the coefficient where ΔS_l described above, has the dimension of an area and depends on the microphone arrangement on the given layer. The order of truncation M=min_LQƒ is the smallest number between the number of loudspeakers L and the number of Q_ƒ of microphones on the layer ƒ. The subscript index [•]_ƒ is due to the relevant fact that, in general cases, the eigenfunctions and eigenvalues p_n,ƒ(x) a_n,ƒ(y) and σ_n,ƒ can be different for different layers. Considering finally equation (31) the frequency response of the filter corresponding to the l-th loudspeaker and the q-th microphone on the layer considered is therefore defined by

$F_{l,q,f}\left(\omega \right)=e^{i\mathrm{\omega \delta }t}\Delta S_l\Delta S_q^{\prime }\sum _{n=1}^M\frac{p_{n,f}\left(x_q,\omega \right)a_{n,f}\left(y_l,\omega \right)}{{\sigma }_{n,f}\left(\omega \right)+\beta \left(\omega \right)}$

In the special case when both the loudspeakers and the microphones on the considered layer are regularly arranged on two spheres of radius R_Λ and R_ƒ, respectively, the filter is computed from

$F_{l,q,f}\left(\omega \right)=e^{i\mathrm{\omega \delta }t}\frac{4\pi }{{\mathrm{LQ}}_f}\sum _{n=1}^{M^{\prime }}\frac{\left(2n+1\right)P_n\left(\cos \left({\zeta }_{\mathrm{lq}}\right)\right)}{{\mathrm{jkj}}_n\left({\mathrm{kR}}_f\right)h_n^{\left(1\right)}\left(k{R}_{\Lambda }\right)+\beta \left(\omega \right)}$ where equations (21), (33), (A2) and (A3) have been used and cos(ζ_lq) is the cosine of the angle between the vectors identifying the locations of the microphone and of the loudspeaker considered (refer to equation (28)). The order of truncation M′ is chosen to beM′≦√{square root over (min_LQƒ)}−1

The outputs of the digital filters are finally combined as shown in FIG. 5.

Auralisation Processing Unit (APU)

FIG. 6 reports a block diagram of the Auralisation Processing Unit. This module is designed to generate the L loudspeaker signals which allow the reproduction of the sound field due to a point source located at a given position in a given reverberant environment, such as a concert hall or an enclosure. The digital filters labeled G₁, G₂, . . . , G_L in FIG. 6 are computed by combining the filters of the Multiple Input Processing Unit with a set of impulse responses describing the reverberant field of the reverberant environment considered. These impulse responses are the impulse responses measured between a reference source (for example a loudspeaker or an omnidirectional source) and the microphone array 100 described above. Both the reference source and the microphone array are located in the reverberant environment considered and the measurements can be carried out using one of the conventional standard techniques mentioned above. The set of Q measured impulse responses are transformed in the frequency domain by applying a Fourier transform, and are labeled R₁, R₂, . . . , R_Q. The filter corresponding to the l-th loudspeaker is computed from

$G_l\left(\omega \right)=\sum _{q=1}^Q{R}_q\left(\omega \right)F_{l,q,f}\left(\omega \right)BP{F}_{f\left(q\right)}\left(\omega \right)$ where the filters F_l,q,ƒ(ω) are defined in the same way as for the MIPU. The Band Pass Filter BPF_ƒ(q)(ω) depends on the layer on which the q-th microphone is arranged.

When designing a Finite Impulse Response filter from the formulation of G_l(ω) given above, it is important to consider that while the filters F_l,q,ƒ(ω) and BPF_ƒ(q)(ω) are in general short in the time domain, the impulse responses R_q(ω) are in general very long, their length depending on the reverberation time of the measured reverberant environment. This factor is vital when defining the filter length, which must be the same if the filters are defined in the frequency domain. In order to avoid this difficulty it is also possible to define the filters in the time domain as

$G_l\left(t\right)=\sum _{q=1}^Q{𝔍}^{-1}\left[R_q\left(\omega \right)\right]\otimes {𝔍}^{-1}\left[F_{l,q}\left(\omega \right)BP{F}_{f\left(q\right)}\left(\omega \right)\right]$ where the operator ℑ⁻¹[•] represents the inverse Fourier transform and the symbol represents a convolution in the time domain.

The Auralisation Processing Unit shares some strong conceptual similarities with the Multiple Input Processing Unit, but while the input to the latter is a stream of Q audio channels which are processed by a matrix of Q by L by F digital filters, the input to the APU is a single audio signal, processed by a bank of L filters. The latter are computed from set of measurements, but their computation can be made off-line. This implies that the real time implementation of an MIPU is much more computationally expensive than that of an APU.

APPENDIX

Spherical Harmonic Expansion of the Free Field Green Function

$\begin{array}{cc}\frac{e^{ikx-y}}{4\pi x-y}=\sum _{n=0}^{\infty }ik{h}_n^{\left(1\right)}\left(ky\right)j_n\left(kx\right)\sum _{m=-n}^n{Y}_n^m\left({\theta }_x{\varphi }_x\right){Y_n^m\left({\theta }_y,{\varphi }_y\right)}^{*}& \left(A1\right)\end{array}$ Finite Summation of Legendre Polynomials

$\begin{array}{cc}\sum _{n=0}^N\left(2n+1\right)P_n\left(x\right)=\left(N+1\right)\frac{P_N\left(x\right)-P_{N+1}\left(x\right)}{1-x}& \left(A2\right)\end{array}$ Summation Formula for the Spherical Harmonics

$\begin{array}{cc}\sum _{m=-n}^n{Y}_n^m\left(\theta ,\varphi \right){Y_n^m\left({\theta }^{\prime },{\varphi }^{\prime }\right)}^{*}=\frac{2n+1}{4\pi }{P}_n\left(\cos \left(\zeta \right)\right)& \left(A3\right)\end{array}$ where P_n(•) is the Legendre polynomial of degree n and ζ is the angle between the directions identified by θ,φ and θ′,φ′. It holds thatcos(ζ)=cos(φ)sin(θ)cos(φ′)sin(θ′)+sin(φ)sin(θ)sin(φ′)sin(θ′)+cos(θ)cos(θ′)=sin(θ)sin(θ′)cos(φ−φ′)+cos(θ)cos(θ′)  (A4)

Orthogonality of the Spherical Harmonics∫₀^2πdφ∫₀^πY_n^m(θ,φ)Y_n′^m′(θ,φ)*sin(θ)dθ=δ_nn′δ_mm′  (A5)

Completeness Relation of the Spherical Harmonics

$\begin{array}{cc}\sum _{n=0}^{\infty }\sum _{m=-n}^n{Y}_n^m\left(\theta ,\varphi \right){Y_n^m\left({\theta }^{\prime },{\varphi }^{\prime }\right)}^{*}=\delta \left(\varphi -{\varphi }^{\prime }\right)\delta \left(\cos \left(\theta \right)-\cos \left({\theta }^{\prime }\right)\right)& \left(\mathrm{A6}\right)\end{array}$ Large Argument Approximation of Spherical Hankel Functions (x→∞)

$\begin{array}{cc}{h}_n^{\left(1\right)}\left(x\right)\approx {\left(-i\right)}^{n+1}\frac{e^{ix}}{x}& \left(A7\right)\end{array}$

Claims

1. A method of determining control signal data for an array of loudspeakers, the control signal data being such as to control the loudspeakers to produce a desired sound field associated with an audio signal, the method comprising: determining control signal data for different frequency components of a desired sound field in respect of respective different positions in a listening volume of a loudspeaker array;wherein, for each of the different frequency components, determination of the control signal data includes virtual sampling of the desired sound field at the surface of a frequency dependent virtual control volume (V); andwherein the size of the control volume (V) is dependent on the frequency component being sampled such that a larger control volume (V) is used for lower frequency components and a smaller volume is used for higher frequency components.

2. The method of claim 1 wherein different sizes of the control volume (V) are used for different frequency bands.

3. The method of claim 1 wherein different sizes of the control volume (V) are of substantially the same shape.

4. The method of claim 1 wherein the different positions define sampling points and wherein substantially the same number of sampling points is used for different sizes of the control volume (V).

5. The method of claim 1 wherein the average distance between sampling points is less than half the wavelength considered.

6. The method of claim 1, further comprising sampling the desired sound field internally of the control volume (V).

7. The method of claim 6 wherein sampling internally of the control volume (V) is effected at those frequencies which correspond to Dirichlet eigenvalues of the control volume.

8. The method of claim 6 wherein a sample internally of the control volume (V) is taken at the geometric centre of the control volume.

9. The method of claim 1 wherein determining the control signal data includes determining the control signal data such that only a loudspeaker, of the loudspeaker array, closest to a virtual source generated by the loudspeaker array is activated.

10. The method of claim 9 wherein determination of the control signal data is used for high frequency component signal processing.

11. The method of claim 1 wherein, in addition to an input audio signal, data indicative of location, radiation pattern and orientation of a virtual source, generated by the loudspeaker array, are used to determine the control signal data.

12. The method of claim 11 wherein the audio signal comprises a signal representative of the sound of an acoustic source and a signal representative of characteristic data of an acoustic environment.

13. The method of claim 12 wherein the sound of the acoustic source and the data of the acoustic environment are included in a signal obtained by sampling the acoustic source in the acoustic environment.

14. The method of claim 12 wherein the sound of the acoustic source and the data of the acoustic environment are separate signals.

15. A sound reproduction apparatus for processing an audio signal, comprising: a microphone array for obtaining an audio signal; anda signal processor arranged to process the audio signal and to output control signal data for an array of loudspeakers to produce a desired sound field associated with the audio signal;wherein the signal processor is configured to determine the control signal data for different frequency components of the desired sound field in respect of respective different positions in a listening volume of the loudspeaker array;wherein, for each of the different frequency components, determination of the control signal data comprises virtual sampling of the desired sound field at the surface of a frequency dependent virtual control volume (V); andwherein the size of the control volume (V) is dependent on the frequency component being sampled such that a larger control volume (V) is used for lower frequency components and a smaller volume is used for higher frequency components.

16. The apparatus of claim 15 wherein the signal processor includes one or more impulse response filters arranged to process the audio signal and output the control signal data.

17. The apparatus of claim 15 wherein the microphone array includes a plurality of microphones, each at a different distance from the centre of the array.

18. A method for processing an audio signal to produce a desired sound field via an array of loudspeakers, comprising: distinguishing different frequency components of a desired sound field;determining control signal data for each of the different frequent frequency components of the desired sound field in respect of different positions in a listening volume of a loudspeaker array by virtual sampling of the desired sound field at the surface of a frequency dependent virtual control volume (V), the size of the control volume (V) is dependent on the frequency component being sampled such that a larger control volume (V) is used for lower frequency components and a smaller volume is used for higher frequency components; andoutputting the control signal data to the array of loudspeakers.