Design of a directive source using the radiation mode method and a soundzone algorithm
Manuel MELON(1) , Maryna SANALATII(2) , Philippe HERZOG(3) , Régine GUILLERMIN(4) , Nicolas POULAIN(2) & Jean
Christophe LE ROUX(2)(1) LAUM UMR CNRS 6613, Av. Olivier Messiaen, 72085 Le Mans Cedex 9, France, manuel.melon@univ-lemans.fr(2) Centre de Transfert de Technologie du Mans (CTTM), 20 rue Thalès de Milet, 72000 Le Mans, France(3) ARTEAC-LAB SAS, Place de l’Innovation, 8 allées Léon Gambetta, 13001 Marseille, France(4) Laboratoire de Mécanique et d’Acoustique (LMA), UMR CNRS 7031, 4 impasse Nikola Tesla, 13000 Marseille, France

Abstract
A common mean to achieve narrow directive patterns consists in using loudspeaker arrays. Generally, theposition of the sources in the array is arbitrarily ﬁxed, for instance by using a regular spacing. The ﬁlters foreach speaker are then calculated to get the desired acoustic ﬁeld. In some cases, the predeﬁned speaker positions
might not be the optimal ones. Moreover, the diffraction on the source body may not be negligible for medium
range applications. In this paper, we propose a two-step approach to overcome these problems. First, theradiation modes of the enclosure are computed. They are used in an inverse scheme, to ﬁnd the most efﬁcientvibration patterns over the enclosure, able to generate a target ﬁeld at some chosen frequencies. These vibration
patterns are then used to deﬁne the positions and sizes of actual loudspeakers on the enclosure. In a second step,a sound zone algorithm is used to calculate the signals driving the speakers. This yields a very narrow directivity
pattern with reduced side lobe amplitudes. The proposed method is tested both numerically and experimentallyand results are reported. A discussion about the pros and cons of the proposed approach is then provided.
Keywords: Sound, Insulation, Transmission

1

INTRODUCTION

The problem of acoustical sound zoning (SZ) is widely studied in the literature. Generally, the SZ techniques
use arrays of loudspeakers which are intuitively placed, what is not always an optimal choice [9, 25, 7, 26, 5].
The problem of sources position optimisation has also been addressed in the literature. Two main strategies
were proposed: (1) selection of the most suitable positions among proposed ones [3, 17, 15] or (2) optimisationof the positions [13, 4, 1, 2, 8, 24] via minimisation of a cost function or by using the gradient method. In thiswork we propose to design the source through the radiation mode (RM) method [23], which allows to obtain
volume velocity distribution over a surface, adapted for a given radiation problem. The capacity of the RM
method to recreate far-ﬁeld source radiation has been shown in previous works [11, 20, 21, 22, 19]. Here we
propose to use the provided velocity distribution over the surface as an indication of the optimal number andpositions of the loudspeakers, as well as their sizes. A brief presentation of the method is given in Sec 2, thedescription of the problem and numerical simulations are presented in Sec 3 and validated by measurements in
Sec. 4. Results are summarized in Sec. 5.

2

THEORY

2.1 Radiation mode method
The “radiation modes” method is inspired by an inverse acoustic imaging technique: the IBEM-NAH regularized by singular value decomposition (SVD) [10] that has been proposed by Veronesi and Maynard [27]. As
the derivation of the method has already been described in [21], only a brief summary will be given here.

1077

The pressure radiated by a vibrating object of surface S can be calculated from the following integral formulation [6, 16]: ∂G∂ ps 
(x −xs ) ps − G(x −xs )p(x) =
dS
(1)
∂n
∂n
S

where x is the location of a point anywhere outside S, xs is the location of an element dS of S, ps = p(xs ) is theparietal pressure at xs , ∂ n denotes the normal derivative on dS and G is the usual Green function G(x − xs ) =1 e jk|x−xs |− jωt time convention.
4π |x−xs | , assuming a e
Formulation 1 may be extended over the surface S, thus relating the surface pressure ps to the normal component of the vibration pattern ∂∂pns over S. This integral equation can then be solved numerically by discretizing
the source surface S with an adequate mesh. Assuming that the pressures and velocities are constant over eachmesh element, the formulation 1 may be reduced to the following matrix form [14, 28]:1
ps = Mqs − Dps2(2)

where qs is the discretization of the normal component of the vibration pattern over S. A pressure ps over S
can then be expressed as ps = Zqs , where the impedance matrix Z is computed as:

Z=

−1
1
I+D
M
2

(3)

where I is the identity matrix.
The impedance matrix Z deﬁned by eq. 3, which is independent of any actual vibration pattern, may then be
expanded trough a SVD. The resulting singular vectors are independent solutions of the associated radiation
problem, called the radiation modes (RMs). The present approach assumes that the target pressure zone is in
the far-ﬁeld, where only the propagative terms are signiﬁcant. The imaginary part of the impedance matrix,
corresponding to the evanescent ﬁeld, can thus be neglected and the RMs can therefore be calculated by SVDof the real part only of the impedance matrix:ℜ(Z) = UΣV∗

(4)

where ∗ means the conjugate transpose of a matrix. U and V are sets of vectors corresponding to the pressure
and the volume velocity of each "mode" respectively, and Σ is a diagonal matrix composed of positive realvalues related to the radiation efﬁciencies of the RMs.
2.2 Expansion coefﬁcients
The far-ﬁeld pressure anywhere around the source pt can therefore be expressed as the superposition of N
RMs: pt = Ht w where Ht is the transfer matrix relating RMs over S to the corresponding pressures pt at
the target locations. The weight coefﬁcients w determine the excitation of each term of the series. Using theMoore-Penrose pseudoinverse of transfer matrix Ht between modes and pressures pt , the coefﬁcients w can be
determined as:
H∗t
pt
(5)
w=
 Ht 2
A simple regularization method consists in truncating the expansion to the N most efﬁcient RMs, where N ischosen according to the desired accuracy (eg by keeping, say, 98% of the cumulated efﬁciency in matrix Σ).
Equation 5 is thus restricted to rank N when computing w.

1078

This allows to estimate a volume velocity q distribution over the surface S, which radiates a desired pressure
pt , with maximum efﬁciency. Here q is simply obtained as the superposition of the RM volume velocities V
weighted by the coefﬁcients w (eq. 6), and it will be used to identify optimal positions of discrete sources forpressure synthesis.
q = Vw(6)
2.3 Sound zone as a pressure target
We now consider two zones: one with high acoustic energy ("bright zone" with pressure pB ) and the other withlower pressure level ("dark zone" with pressure pD ). In the present study, the target pressure pt is deﬁned overthe bright zone, and only a mean pressure level D0 is deﬁned in the dark zone. Therefore, we need to relax
the constraints expressed at all points by equation 5 when computing w. It is here proposed to use a classicalacoustic contrast control (ACC) technique. We chose the pressure matching method [5] which objective is to
impose a target pressure with maximum level in the bright zone, while creating an acoustic contrast between
the bright and the dark zones, by minimizing of the following cost function J:
J = [pB − pt ]∗ [pB − pt ] + [p∗D pD −D0 ](7)

The optimal source strength qopt , that minimizes the cost function J is:
qopt = [G∗B GB + μG∗D GD ]−1 [G∗B pt ](8)here μ is a positive real coefﬁcient which tunes the relative importance of the contrast condition, GB and GD
are the transfer matrices relating source strengths with pressures at MB or MD points, respectively located inbright and dark zones.
This technique can be incorporated into the RM method, according to the similarities between pressure representation through the source strength in ACC techniques and the one through the RMs:
p = Gq = Hw

(9)Here, G is the transfer matrix relating N loudspeakers membrane volume velocities q with pressures p in M
points, and H is the transfer matrix between N surface RMs and pressures p in M points. Hence, instead ofdriving array sources by the adapted signals qopt from eq. 8, one can "drive" the series of modes using adapted
RM coefﬁcients wopt determined as:
wopt = [H∗B HB + μH∗D HD ]−1 [H∗B pt ](10)

q = Vwopt(11)

Eq. 6 is thus rewritten as:

3

SIMULATIONS

3.1 Test caseAs a ﬁrst case study, we use a 53 × 53 × 18 cm enclosure placed at 70 cm in front of a rigid wall. Acousticallybright zone is a 40 × 40 cm square spot on the wall in front of enclosure. The target pressure level is 100 dB.
A 1 m diameter zone around bright zone and 1 m diameter cylinder around enclosure ended by an half-spherecompose the dark zone, in which the pressure level should be attenuated by at least 15 dB. The problem is
depicted in Figure 1.
For the numerical calculation, the enclosure surface is discretized by a mesh with 3 mm elements, bright zone
by MB = 81 points regularly spaced by 5 cm, and dark zones on the wall (D1) and around the surface (D2)

1079

1

Z, [m]

0.5
0
−0.5
−1
1

1

0
Y, [m]

Figure 1. Problem schema

−1

0
−1

X, [m]

Figure 2. Numerical representation of the problem

are presented by MD1 = 864 and MD2 = 2324 points respectively (Figure 2). The rigid wall is modelled as an
inﬁnite reﬂecting plane.
The RMs and their corresponding HB + HD matrices were computed using a custom BEM code developed at
LMA, dedicated to ﬁeld expansion ("FELIN": Field ELement INdependence).
3.2 Modal volume velocity interpretation
For now, analysis of the modal velocity distribution is a tricky task for two main reasons. First, a trade-offmust be found to determine the number N of RMs, as it must allow an accurate enough description without
leading to complex vibrating shapes which could hardly be reproduced in practice. A second challenge is toconvert the frequency-dependant velocity distributions over surface S into a small set of frequency-independent
discrete sources.
First simulations showed that up to 1 kHz, 20 to 35 modes can be sufﬁcient for sound reconstruction. Thus, wetook 35 modes to estimate the distributions of the modal volume velocity over S. They are presented in Figure
3 at several frequencies (200 Hz, 500 Hz, 750 Hz and 900 Hz). Visual analysis helped to establish the generalfeatures of the velocity distributions. Almost for all frequencies there are 5 accentuated zones at the rear sideof the structure: 4 zones in-phase at the corners and the 1 in the center. The locations for the front side is lessobvious: At most frequencies it does not seem to play any important role in the generation of sound. Though,at medium frequencies its contribution becomes very signiﬁcant, with a vibrating central zone of size varying
with frequency.
After few tests that are detailed in [18], a 4 channels system has been deﬁned: (1) 1 loudspeaker in the front
side center (Audax HM 170 GO with membrane diameter of 13 cm), (2) 16 small diameter loudspeakers building
an annular zone around the center (Aura NSW2-326-8A with membrane diameter of 4, 4 cm), (3) 1 loudspeakeron the rear side center Beyma 5P200/Fe, membrane diameter of 12 cm, (4) 4 loudspeakers at the rear side
corners Dayton RS100-8, membrane diameter of 7.7 cm. On the front side, the combination of the center loudspeaker with the ring of smaller ones allows to vary the size of the vibrating zone. This source arrangement is
inspired by a previous work where a loudspeaker ring was used to obtain a wide reproduction zone [12].
After selecting the transducer locations over the enclosure, the driving signals may be computed by one of thesound zone techniques as for a conventional array. In the present work we use again the pressure matching
technique. The optimal source strength vector qopt is thus obtained from equation 8.

1080

−0.2

−0.2

0
0.2
x
Rear : Level

1

−0.2

0
0.2
x
Rear : Phase

−0.2

−0.1

−0.1

0
−0.5

−108 −0.2
−110
−112

−0.2

0
0.2
x
Rear: Level

−0.2
−0.1

0

0

−98

0.1

0.1

−2

−100 0.2

0.2−118 0.1
0.2
−120
−0.2

0.2

−0.2

0
x

0.2

y

−116

y

−1.5

y

0

0
x

0
x

0.2

−114 −0.2
−0.1
−116
0
−118 0.1

−0.2

0
0.2
x
Rear : Phase

−0.2

0

−1

y

−0.1

y

−120 0.2
−0.2

0

−2

0.1
0.2
0
x

0.2

(c) 750 Hz

0.1

−0.2

0
x

0.2

−3

y

−102 0.1
0
−104
−0.1

0.2

−108
−110
−112
−114

3

2

0
−0.1

−106 −0.2

1

−0.2

−0.2

0
0.2
x
Rear : Level

−0.2

0
0.2
x
Rear : Phase

−0.2

−0.2

−0.1

−0.1

0

−116 0.1
0.2
−118
−0.2

y

1

−0.2
0
0.2
x
Rear : Level

0x

Front : Phase
0.2

y

2

−0.1

−0.2

−0.2

−100 0.2

Pressure level [dB]Pressure level [dB]

−112

3

y

0

−110

−2

0.2

Front : Level

Phase [rad]

0.1

0

y

−104 0.1
y

0.2

−0.1
−108 −0.2−1(b) 500 Hz
Front : Phase

0.2

−106

0

0.1

(a) 200 Hz
Front : Level
−102

0
0.2
x
Rear : Phase

−0.1

0

−0.2

−0.2

−0.2

−114

−2.5

1

−0.2

−1

−96

2

0
−0.1

y

−94 −0.2

0.5

0.1

Phase [rad]

−92

−0.2

−104 0.1
0
−106
−0.1

3

−1

0

−2

0.1
0.2
0
x

0.2

−0.2

0

Phase [rad]

−90

−0.1

Pressure level [dB]

−88

−0.1

1.5

0.2
y

0

−102 0.2
y

0.1

0

2

[ d]

0.1

−86

Front : Phase

Front : Level

2.5

Ph

−84

y

0.2

y

0.2

Pressure level [dB]

Front : Phase

Front : Level
−82

0
x

0.2

−3

(d) 900 Hz

Figure 3. Modal volume velocity distribution over the surface for different frequencies(level and phase, front and rear view)
3.3 Design validation
To control the performance of the system, we propose to use 3 metrics. Global acoustical contrast CG , describes
the difference of the acoustical power in the bright and dark zones. Higher value of this metric indicates better
performance. Expressed in decibels, the global contrast is [9]:MD pB ∗ pB
CG = 10 log10(12)×
MB pD ∗ pD
Similarly, the local acoustical contrast CL is deﬁned between the pressure over the bright zone zone, and theone over the dark zone D1 , both in front of the system.
pB ∗ pB
MD1
CL = 10 log10
×
(13)
MB
pD1 ∗ pD1
The normalised reproduction error compares the desired pt and reproduced pB ﬁelds in the bright zone [7]:(pt − pB )∗ (pt − pB )ER = 10 log10(14)pB ∗ pB
For this last metric, a smaller value indicates a better performance, with 0 dB corresponding to a high error
level (equal to the target ﬁeld level).

1081

25

20

20

15

15

G

10
5

5

0

0
−5

−20
−30

−10 2
10

3

10

Frequency, Hz

−10

10

−5
−10 2
10

0

ER, dB25

CL, dBC , dB

Figure 4 presents the simulation of the 4-channel source metrics while Figure 5 shows the pressure ﬁeld radiatedat MB and MD points for several frequencies. We can observe the ability of proposed source to maintain analmost constant contrast level between the zones at all frequencies as well as good synthesis accuracy.

(a) Global contrast

−40 2
10

3

10

Frequency, Hz(b) Local contrast

3

10

Frequency, Hz(c) Reproduction error

Figure 4. Metrics for the numerical model

0
−0.5

85

−1

80
1

1

0
−1

0
−1

(a) 200 Hz

75

0.5

90

0
−0.5

85

−1

80
1

1

0
−1

Pressure level, dB

90

Pressure level, dB

0.5

100

100
1

95

75

0
−1

(b) 500 Hz

1

95

0.5

90

0
−0.5

85

−1

80
1

1

0
−1

Pressure level, dB

1

95

75

0
−1

(c) 750 Hz

95

0.5

90

0

85

−0.5
−1

80
1

1

0
−1

Pressure level, dB

100

100
1

75

0
−1

(d) 900 Hz

Figure 5. Radiated pressure maps of the numerical model

4

REALISATION AND MEASUREMENT

An experimental validation was held in the CTTM semi-anechoic chamber. The system was ﬁxed on a turning table ﬂush mounted in the room ﬂoor. The pressure behind the system was measured thanks to a wooden
support equipped with 18 microphones (Figures 6 and 7). Assuming symmetrical radiation of sources, measure1of space were taken, by rotating the system 90◦ around its axis. Pressure measurements (zones
ments of the4
B and D1) were also recorded on the ﬂoor by using a metallic stem holding 11 microphones spaced by 10 cm(Figure 8).
The driving signals were generated by a RME Fireface US audio interface and ampliﬁed for different groups ofloudspeakers by a multichannel ampliﬁer HPA QA4150. Pressure ﬁeld was measured with 29 14 ” PCB 130D21
microphones and acquired through National Instruments boards NI PXI-1045 using the INTAC software.
The driving commands were computed from eq. 8, using measured transfer functions between 4 channels andmicrophones. Satisfactory results are presented in Figures 9 and 10, showing a constant contrast level and accurate target pressure synthesis. A discrepancy between desired and real pressure generation in the bright zone canhowever be noticed at 200 Hz. It is probably due to the differences between the 4 actual corner loudspeakers,
assumed identical and identiﬁed as a single source during the computation of the commands.

1082

Figure 7. Photo of themeasurements setup
25

20

20

15

15

5

0
−5

10

R

G

10

Figure 8. Photo of thesource on the turning table

E , dB25

CL, dBC , dB

Figure 6. Schema of themeasurement setup

5

0

0
−5
−10 2
10

3

10

−15
−20

−5
−10 2
10

Frequency, Hz−10(a) Global contrast

−25 2
10

3

10

Frequency, Hz(b) Local contrast

3

10

Frequency, Hz(c) reproduction error

Figure 9. Metrics for the experimental mock-up
100

100

100

95

0

85

−0.5

80

−1

75

1

1

0
−1

0
−1

0.5

90

0
85

−0.5

80

−1
1

70

(a) 200 Hz

1

0
−1

0
−1

(b) 500 Hz

75

1

95

0.5

90

0
85

−0.5

80

−1
1

1

0
−1

95

0.5

90

0
85

−0.5
−1

75

0
−1

(c) 750 Hz

80
1

1

0
−1

Pressure level, dB

90

1

95

Pressure level, dB

0.5

Pressure level, dB

1

Pressure level, dB

1

75

0
−1

(d) 900 Hz

Figure 10. Measured maps of the radiated ﬁeld for several frequencies

5

CONCLUSION

The aim of the present work consisted in designing a directive source able to generate a predeﬁned ﬁeld witha high spatial contrast. In order to reduce the number of loudspeakers to the minimum possible, we used theradiation mode method as an indication of the most suitable transducer locations and sizes. The visual analysis
has permitted to design a source, that responds to all imposed conditions. First numerical and experimental
results are very promising, and motivate further investigations. For instance, genetic algorithms could be used
to select the optimal number, diameters and locations of speakers on the enclosure using the RMS.

ACKNOWLEDGEMENTS
Authors would like to thank Marie-Hélène Moulet for conducting the measurements and the CTTM design
engineering department for the source production.

1083

REFERENCES[1] F. Asano, Y. Suzuki, T. Sone, and D. Swanson. Optimization of control source location in active control systems. In Proceedings of Active 95, pages 489–498,
1995.
[2] F. Avanzini, A. Belussi, A. Dal Palù, A. Dovier, and D. Rocchesso. Optimal placement of acoustic sources in a built-up area using CLP (FD). In APPIA-GULPPRODE 2002, Joint Conf. on Declarative Programming, pages 139–154, 2002.
[3] K. Baek. Natural algorithms for choosing source locations in active control systems. Journal of Sound and Vibration, 186(2):245–267, 1995.
[4] E. Benzaria and V. Martin. Constrained optimization of secondary source locations: multipolar arrangements. In Proceedings of Active 95, pages 499–510, 1995.
[5] T. Betlehem, W. Zhang, M. A. Poletti, and T. D. Abhayapala. Personal sound zones: Delivering interface-free audio to multiple listeners. IEEE Signal Processing
Magazine, 32(2):81–91, 2015.
[6] M. Bruneau. Fundamentals of acoustics. CA: ISTE Ltd, 2006.
[7] J. Chang and F. Jacobsen. Sound ﬁeld control with a circular double-layer array of loudspeakers. The Journal of the Acoustical Society of America, 131(6), 2012.
[8] G. N. Charalampopoulos, S. T. Mouzakitis, and C. G. Provatidis. Gradient optimisation methods for the positioning of active noise control actuators in enclosures.
International Journal of Acoustics and Vibrations, 9(4):163–174, 2004.
[9] J.-W. Choi and Y.-H. Kim. Generation of an acoustically bright zone with an illuminated region using multiple sources. The Journal of the Acoustical Society ofAmerica, 111(4):1695–1700, 2002.
[10] G. H. Golub and C. F. V. Loan. Matrix Computations, volume 10. Physics Today, 1996.
[11] P. Herzog, R. Guillermin, P. Lorin, and V. Chritin. Identiﬁcation of a vibration pattern from pressure measurements and radiation modes. In EuroNoise, Maastricht,
Netherlands, 2015.
[12] P. Herzog and M. Melon. An active cell to improve transmission loss of partition walls at lower frequencies. In Euronoise, Prague, Czech Republic, 10-13 June
2012.
[13] J. K. Kim and J.-G. Ih. On the positioning of control sources for active noise control in 3-dimensional enclosed space. In Proceedings of Active 95, pages
510–518, 1995.
[14] S. Kirkup. The Boundary Element Method in Acoustics. Integrated Sound Software, 1998.
[15] D. A. Manolas, T. Gialamas, and D. T. Tsahalis. A genetic algorithm for the simultaneous optimization of the sensor and actuator positions for an active noise
and/or vibration control system. In Proceedings of Internoise 96, pages 1187–1192, Liverpool, UK, 30 July - 2 August 1996.
[16] P. M. Morse and K. U. Ingard. Theoretical acoustics. Princeton University Press, 1968.
[17] S. Pottie and D. Botteldooren. Optimal placement of secondary sources for active noise control using a genetic algorithm. In Proceedings of Internoise 96, pages
1101–1104, Liverpool, UK, 30 July - 2 August 1996.
[18] M. Sanalatii. Synthèse d’un champ acoustique avec contraste spatial élevé. PhD thesis, Le Mans University, 2018.
[19] M. Sanalatii, P. Herzog, R. Guillermin, M. Melon, N. Poulain, and J.-C. Le Roux. Estimation of loudspeaker frequency Response and directivity using theradiation-mode method. J. Audio Eng. Soc., 67(3):101–115, 2019.
[20] M. Sanalatii, P. Herzog, M. Melon, R. Guillermin, J.-C. Le Roux, and N. Poulain. Identiﬁcation du proﬁl vibratoire d’un haut-parleur monté sur une enceinte
close a partir de mesures de pression acoustique en champ proche. In Congrès Français d’Acoustique, number 215, Le Mans, France, 2016.
[21] M. Sanalatii, P. Herzog, M. Melon, R. Guillermin, J.-C. Le Roux, and N. Poulain. Measurement of the frequency and angular responses of loudspeaker systems
using radiation modes. In the 141th Audion Engineering Society Convention, number preprint 9615, Los Angeles, USA, 2016 September 29th - October 1st.
[22] M. Sanalatii, L. Vindarola, C. Vasseur, P. Herzog, M. Melon, R. Guillermin, N. Poulain, and J.-C. Le Roux. Assessment of the radiation mode method for in situ
measurements of loudspeaker systems. In the 142nd Audion Engineering Society Convention, number preprint 9615, Berlin, Germany, 2017 May 20-23.
[23] A. Sarkissian. Acoustic radiation from ﬁnite structures. The Journal of the Acoustical Society of America, 90(1):574–578, 1991.
[24] P. Sergent and D. Duhamel. Optimal placement of sources and sensors with minimax criterion for active contriol of a one-dimentional sound ﬁeld. Journal ofSound and Vibration, 207(4):537–566, 1997.
[25] M. Shin, S. Q. Lee, F. M. Fazi, P. A. Nelson, D. Kim, S. Wang, K. H. Park, and J. Seo. Maximization of acoustic energy difference between two spaces. The
Journal of the Acoustical Society of America, 128(1):121–131, 2010.
[26] M. F. Simón Gálvez, S. J. Elliott, and J. Cheer. A superdirective array of phase shift sources. The Journal of the Acoustical Society of America, 132(2):746–756,
2012.
[27] W. Veronesi and J. Maynard. Digital holographic reconstruction of sources with arbitrary shaped surfaces. Journal of the Acoustical Society of America,
85(2):588–598, 1989.
[28] S. F. Wu. The Helmholtz Equation Least Squares Method. Springer, 2015.

1084