Three-dimensional spatial active noise control based onkernel-induced sound field interpolation
Hayato ITO(1) , Shoichi KOYAMA(1),(2) , Natsuki UENO(1) , Hiroshi SARUWATARI(1)(1) The University of Tokyo, Graduate School of Information Science and Technology, Japan, hayato_ito@ipc.i.u-tokyo.ac.jp(2) JST, PRESTO, Japan

Abstract
A method for feedback active noise control (ANC) over a three-dimensional (3D) spatial region is proposed.
Conventional multipoint ANC does not guarantee to reduce the noise between multiple discrete control points.
Several attempts have been made to reduce the noise over the continuous target region. Most methods for
spatial ANC found their basis on the spherical/cylindrical harmonic expansion of the sound field. Therefore,
they can be applied only for arrays of simple geometries such as sphere and cylinder. We propose a 3D spatial
ANC method based on sound field interpolation, which can be applied to arbitrary array geometries. The
mathematical principle of the proposed method is based on our recent studies on kernel ridge regression forsound field interpolation. The difference in the resulting algorithm from the general adaptive control appears
only in the weighting matrix for error signals, which can be computed in advance only by using the relative
positions of the microphones and the target region. Numerical simulation results indicated that a large regional
noise reduction is achieved by the proposed method compared with conventional multipoint ANC.
Keywords: Active noise control, Sound field interpolation, Sound field control, Kernel method

1

INTRODUCTION

Active noise control (ANC), or active noise cancellation, aims to cancel incoming primary noise by using secondary actuators (loudspeakers) (1, 2). ANC has been mainly investigated for reducing the noise at predefined
single or multiple control points using adaptive control techniques when the noise sources and secondary loudspeakers are placed apart from a target control region. Many practical ANC systems have also been developed
based on these techniques (3, 4, 5, 6).
Multichannel ANC algorithms, where multiple secondary loudspeakers are used for controlling noise at multiple
positions, have been developed both in time and frequency domains in the literature (7). The power of theerror microphones placed at the control points is generally used as the cost function, and a type of leastmean-square (LMS) algorithms is typically applied for the adaptive control. Two types of adaptive control havebeen investigated: feedforward and feedback control (8). The feedforward ANC system is applicable to reduce
nonstationary noise by using reference microphones. The feedback ANC is effective for periodic noise signals
without using the reference microphones. The main drawback of the multichannel ANC is its limitation of thenoise cancellation only on the multiple discrete positions; therefore, the region between the error microphones
is not taken into consideration.
Sound field recording and reproduction have been investigated for reconstructing a spatial sound field both in
the theoretical analysis (9, 10, 11, 12) and development of practical systems (13, 14). In recent years, thesestudies have led to researches on spatial ANC, which aims to cancel spatial noise in a continuous target region.
Most methods of the spatial ANC is generally based on the spherical/circular harmonic expansions of the soundfield inside the target region (15, 16, 17, 18); however, they can only be applied to simple array geometry suchas sphere and cylinder.
We propose a method for three-dimensional (3D) spatial ANC using distributed microphones and loudspeakers,
which is based on the interpolation using kernel ridge regression with the constraint that the function to be estimated satisfies the Helmholtz equation in the frequency domain (19, 20). We recently proposed a feedforward1101

Figure 1. Arrangement of secondary loudspeakers and error microphones in feedback spatial ANC
spatial ANC in 2D on the basis of the kernel-induced sound field interpolation (21). In this paper, we formulate
a feedback spatial ANC method, which aims to reduce periodic noise inside a 3D target region. The spatial
noise reduction performance of the proposed method is evaluated by numerical simulations compared with themultipoint feedback control.

2

FEEDBACK ACTIVE NOISE CONTROL BASED ON NLMS ALGORITHM

We first revisit the multichannel feedback control based on the normalized LMS (NLMS) algorithm in thefrequency domain. The objective is to reduce the spatial noise from primary noise sources inside a target regionΩ by using multiple microphones and loudspeakers as shown in Fig. 1. L secondary loudspeakers and M error
microphones are arranged around Ω, and the loudspeaker output is determined by an adaptive algorithm tocancel the primary noise field. The driving signals of the loudspeakers and the signals of the error microphones
at the angular frequency ω are denoted as d (ω ) ∈ CL and e (ω ) ∈ CM , respectively. We hereafter omit ωfor notational simplicity. When the sound pressures of the primary noise field at the positions of the error
microphones are denoted as u p ∈ CM , the signals of the error microphones are represented as
e = G d + up ,(1)

where G ∈ CM×L is the transfer function matrix from the secondary loudspeakers to the error microphones. We
assume that G is given by measurements in advance.
The multipoint pressure control (MPC) for the feedback ANC aims to minimize the square ℓ2 -norm of the error
signals (7). Thus, the cost function of MPC JMPC is formulated as
JMPC = ∥ee∥22 = e H e ,(2)where the superscript (·)H denotes the conjugate transpose. In the complex LMS algorithm, the driving signalsof the loudspeakers are adaptively determined by updating in the steepest-descent direction of JMPC at eachstep as∂JMPC∂d∗= d (n) − η G H e (n)d (n + 1) = d (n) − η(3)where the superscript (·)∗ denotes the complex conjugate, n denotes the iteration index, and η is the step size
parameter. The complex NLMS algorithm (22) is frequently used for faster convergence in which d is updated

1102

as
d (n + 1) = d (n) −η0
G H e (n),GH G ∥2 + ε∥G(4)where ε > 0 is the regularization parameter and η0 ∈ (0, 2) is the normalized step size parameter.
We here introduce a weighting matrix to the cost function. The cost function J of the square ℓ2 -norm of ewith the weighting matrix A ∈ CM×M , which is independent of d , is represented as
J = eHAe.

(5)

The NLMS algorithm for J can be obtained in a similar manner to the case of JMPC as

η0
G H A e (n).
GH A G ∥2 + ε∥Gd (n + 1) = d (n) −(6)

Our proposed method, which is discussed in the next section, is based on the NLMS algorithm using theweighted-square-ℓ2 -norm. Note that the computational cost for the filter update in Eq. (6) is exactly the same
GH A G ∥2 in advance. The specific A for 3D spatial ANC is given
as that in Eq. (4) by calculating G H A and ∥G
in the following section.

3

3D SPATIAL ANC BASED ON KERNEL-INDUCED SOUND FIELD INTERPOLATION

Our objective is to reduce the primary noise in the entire target region Ω. However, in the conventional MPC,
the noise only at the positions of the error microphones is reduced, and the region between these microphones
is not taken into consideration. We consider the regional noise power for the cost function as in our previously
proposed method for feedforward control (21) as
L =

∫
Ω

|ue (rr )|2 drr ,(7)

where ue (rr ) is the continuous pressure field at the position r , which is the superposition of the primary noise
field and the sound field synthesized by the secondary loudspeakers. Obviously, it is impractical to directly
obtain ue (rr ) by densely arranging the error microphones inside Ω. We apply the kernel interpolation technique
for the sound field (19, 20) to estimate ue (rr ) from the error signals e .
3.1 Weighting matrix for 3D sound field interpolation
A method for estimating spherical harmonic coefficients at an arbitrary position from the pressure measurementsof distributed microphones was proposed by Ueno et al. (20), which is based on the spherical harmonic analysis
of infinite order of a sound field. This method can be regarded as the kernel ridge regression for sound fieldinterpolation (19). We apply this method to estimate ue (rr ) from e . When Ω does not include any sources andthe error microphones are omnidirectional, the estimate ûe (rr ) can be obtained as)T(
K + λ I )−1 e ) κ (rr ),ûe (rr ) = (K(8)

where the superscript (·)T denotes the transpose, I is the identity matrix, λ is the regularization parameter.
Here, K and κ for 3D sound field interpolation are defined as
j0 (k∥rr 1 − r 1 ∥) · · · j0 (k∥rr 1 − r M ∥)


..
..
..
K =
(9)

.
.
.
j0 (k∥rr M − r 1 ∥)···j0 (k∥rr M − r M ∥)

1103

and

κ (rr ) = [ j0 (k∥rr − r 1 ∥), · · · , j0 (k∥rr − r M ∥)]T ,(10)

respectively, where j0 (·) is the zeroth-order spherical Bessel function of the first kind, k = ω /c is the wave
number with the sound speed c, and r m (m ∈ {1, · · · , M}) is the position of the mth error microphone. Notethat the interpolant is determined only by the relative positions of the error microphones without setting the
expansion center and truncation order of the harmonic coefficients required in the general spherical-harmonicexpansion-based methods.
By substituting Eq. (8) into Eq. (7), the cost function L can be written as
L =

∫
Ω

|ûe (rr )|2 drr = eH A e,

(∫

where
H

A=P

Ω

∗

(11)

)

κ (rr )κ (rr )drr P
T(12)

K + λ I )−1 . Since L in Eq. (11) has the form of the weighted ℓ2 -norm of e as in Eq. (5), and A inand P ≡ (KEq. (12) is independent of d , the NLMS algorithm in Eq. (6) can be applied.
3.2 Efficient computation of weighting matrix for spherical target region
When the target region Ω is spherical, it is possible to avoid the computation of the numerical integration in
A . The radius of the spherical region Ω is denoted as R. First, in spherical coordinates r = (r, θ , ϕ ), addition
theorem for spherical Bessel function is expressed as (23)∞j0 (k∥rr − r m ∥) = 4π ∑

ν

µ∗

µ

∑ (−1)ν jν (krm )Yν (−r̂r m ) jν (kr)Yν (r̂r )ν =0 µ =−ν
∞≡ 4π ∑ ξν ,m (rr ),

(13)

ν =0

where ξν ,m (rr ) is defined as −ν Y (r̂r )( −ν)∗  ν . 
ν
−
ν
+1
ν
(−r̂r m ) · · · Yν (−r̂r m )  ..  jν (kr) ≡ α ν ,m β ν (rr ),ξν ,m (rr ) = (−1) jν (krm ) Yν (−r̂r m ) Yν|{z}Yνν (r̂r )α ν ,m
|
{z
}

(14)

β ν (rr )µand Yν (·) is the spherical harmonic function of ν th order and µ th degree. Here, r̂r represents the unit vector
pointing to r , i.e., r̂r = r /∥rr ∥2 . Thus, j0 (k∥rr − r m ∥) can be written as
β 0 (rr )β 1 (rr )
()  .. (15)
j0 (k∥rr − r m ∥) = 4π α 0,m α 1,m · · · α k,m · · ·  . .
β (rr )
k


..
.

1104

By substituting Eq. (15) into Eq. (10), κ (rr ) is represented as

α 0,1
 α 0,2κ (rr ) = 4π  .
 ..
α 0,M
|

α 1,1
α 1,2
..
.

α 1,M

···
···

{z



 β 0 (rr )· · · β (rr )
 1 · · ·
.. 

T
. 


 ≡ S ϕ (rr ).
· · · β (rr )
 k 
···
..
}
.
| {z }

α k,1
α k,2
..
.

···
· · · α k,M

ST

(16)

ϕ (rr )

Therefore, the integral term in A becomes
(∫
H H

A=P S

Ω

∗

)

ϕ (rr )ϕ (rr ) drr S P ,
T

(17)

where
∗

β 0 (rr )β 1 (rr )

 ..  ()∗
T
T β (rr )T · · · β (rr )T · · ·ϕ (rr )ϕ (rr ) =  . 1k β 0 (rr )β (rr )
 k 
..
.

β 0 (rr )∗ β 0 (rr )T β 0 (rr )∗ β 1 (rr )T · · · β 0 (rr )∗ β j (rr )Tβ 1 (rr )∗ β 0 (rr )T β 1 (rr )∗ β 1 (rr )T · · · β 1 (rr )∗ β j (rr )T


..
..
..
..
.
=
.
.
.

 β (rr )∗ β (rr )T β (rr )∗ β (rr )T · · · β (rr )∗ β (rr )T
0
i
1
i
j
 i
..
..
..
..
.
.
.
.
and

 −i ∗Yi (r̂r )()
 ..  ∗∗T
β i (rr ) β j (rr ) =  .  ji (kr) j j (kr) Y j− j (r̂r ) · · · Y jj (r̂r ) .
Yii (r̂r )
···· · ·
· · ·
· · ·

..
.

(18)

(19)

Owing to the orthogonality of the spherical Bessel functions, the integral term of this equation is derived as∫Ωϕ ∗ (rr )ϕ T (rr ) drr = diag(γ0 , γ1 , γ1 , γ1 , γ2 , · · · ) ≡ Γ ,

(20)

where

)
R3 (
(21)
jν (kR)2 − jν −1 (kR) jν +1 (kR) .
2
Here, a mathematical formula in Ref. (24) is used to derive Eq. (21). 2ν + 1 elements of γν is arranged foreach ν as the diagonal components of Γ . Finally, the weighting matrix A for the spherical target region Ω isderived as
A = PHSHΓSP.
(22)

γν =

In practice, it is necessary to approximate the infinite matrices S and Γ as finite matrices of a sufficiently large
size. Larger sizes of S and Γ are better for accurate numerical computation of the integration (25). Note that
A can be calculated in advance only from the relative positions of the error microphones and the target region.

1105

15

MPC
Proposed10

Regional reduction (dB)Regional reduction (dB)

0

-5

-10

5
0
-5
-10
-15
-20

MPC
Proposed

-25

-15

100

200

300

400

500

Number of iterations

Figure 2. Regional noise power reduction inside Ω at 250 Hz with respect to number of
iterations.

4

-30
100

200

300

400

500

Frequency (Hz)Figure 3. Regional noise power reduction inside Ω with respect to frequency.

NUMERICAL SIMULATIONS

Numerical simulations were conducted in a 3D sound field to compare the proposed method (Proposed) andmultipoint pressure control (MPC). The target region Ω was set to be spherical with a radius of 0.5 m, andthe coordinate origin was set at its center. M = 60 omnidirectional microphones were used as the error microphones. The secondary loudspeakers were assumed to be point sources and their number was L = 30. The error
microphones were placed on three dual rings set on the surface of Ω. The intervals of the dual rings were0.01 m. The three rings were set at z = −0.3, 0.0, and 0.3 m, respectively. On each ring, 10 error microphones
were equiangularly placed. The secondary loudspeakers were placed on three rings set on the spherical surfaceof the radius 1.0 m. The three rings were at z = −0.7, 0.0, and 0.7 m, respectively. 10 secondary loudspeakers
were equiangularly placed on each ring. The speed of sound was set as c = 340 m/s.
Six point sources were placed as the primary noise sources at (−10.0, 0.0, 0.0), (10.0, 0.0, 0.0), (0.0, −10.0, 0.0),(0.0, 10.0, 0.0), (0.0, 0.0, −10.0) and (0.0, 0.0, 10.0) m in Cartesian coordinates. The source signals were singlefrequency complex sine waves. Observation noise was also added to the error microphones so that the signal-tonoise ratio (SNR) was 40 dB. For evaluation, the performance measure with regard to the regional noise power
reduction inside Ω is defined as(n)∑ |ue (rr i )|2
Pred = 10 log10 i (0),(23)∑i |ue (rr i )|2(n)where ue (rr i ) is the total pressure field at the ith evaluation position r i for the nth iteration. The evaluation
positions were obtained by discretizing Ω with intervals of 0.025 m. The pressure field at the zeroth iteration(0)ue (rr i ) was the original pressure field generated by the primary noise sources.
We first show the results when the frequency was set at 250 Hz and the amplitude of six noise sources was
10, 20, 30, 40, 50, and 60, respectively. The parameter ε was set as 10−2 . The maximum order ν used forconstructing the matrices S and Γ was 40. The regularization parameter λ was set as 10−4 . The performancemeasure Pred is plotted with respect to the number of iterations for the two methods in Fig. 2. Proposed
achieved larger noise power reduction than MPC.
Figure 3 shows Pred with respect to the frequency after 500 iterations. Proposed and MPC achieved similar
noise reduction when the frequency was below 200 Hz. Above 200 Hz, Proposed outperformed MPC.

1106

5

CONCLUSIONSWe proposed a feedback spatial ANC method using distributed secondary loudspeakers and error microphones,
which is based on the kernel ridge regression of a 3D sound field with the constraint that the function to be
estimated satisfies the Helmholtz equation. The cost of the 3D regional error power is defined, and the kernelinduced sound field interpolation is applied to calculate it. The resulting algorithm includes a weighting matrix
for the error signals obtained in advance by using the relative positions of the microphones and the target region.
Numerical simulation results indicated that the proposed method outperforms the conventional multipoint control
in regional noise reduction.

ACKNOWLEDGEMENTS
This work was supported by JST PRESTO Grant Number JPMJPR18J4, and SECOM Science and Technology
Foundation.

REFERENCES[1] Nelson, P. A.; Elliott, S. J. Active control of sound, Academic Press, London, 1st edition, 1992.
[2] Fuller, C. R.; Elliott, S. J.; Nelson, P. A. Active Control of Vibration, Academic, San Diego, CA, USA,
1st edition, 1996.
[3] Kuo, S. M.; Mitra, S.; Gan, W.-S. Active noise control system for headphone applications, IEEE Trans.
Control Syst. Technol., Vol. 14 (2), 2006, pp. 331–335.
[4] Kottayi, S.; Althomali, R.; Thasleema, T. M.; Narayanan, N. K. Active noise control for creating a quiet
zone around mobile phone, Proceedings of Int. Conf. Comm., Signal Process. (ICCSP), Melmaruvathur,Apr. 2016, pp. 73–77.
[5] Sano H.; Inoue, T.; Takahashi, A.; Terai, K.; Nakamura, Y. Active control system for low-frequency roadnoise combined with an audio system, IEEE Trans. Speech Audio Process., Vol. 9 (7), 2001, pp. 755–763.
[6] Diaz, J.; Egaña, J. M.; Viñolas, J. A local active noise control system based on a virtual-microphone
technique for railway sleeping vehicle applications, Mech. Syst. Signal Process., Vol. 20, (8), 2006, pp.
2259–2276.
[7] Elliot, S. J.; Stothers, I. M.; Nelson, P. A. A multiple error LMS algorithm and its application to theactive control of sound and vibration, IEEE Trans. Acoust., Speech, Signal Process., Vol. 35 (10), 1987,
pp. 1423–1434.
[8] Kuo S. M.; Morgan, D. R. Active noise control: A tutorial review, Proc. IEEE, Vol. 87 (6), 1999, pp.
943–973.
[9] Poletti, M. A. Three-dimensional surround sound systems based on spherical harmonics, J. Audio Eng.
Soc., Vol. 53 (11), 2005, pp. 1004–1025.
[10] Ahrens J.; Spors, S. An analytical approach to sound field reproduction using circular and spherical loudspeaker distributions, Acta Acustica United with Acustica, Vol. 94, 2008, pp. 988–999.
[11] Koyama, S.; Furuya, K.; Hiwasaki, Y.; Haneda, Y. Analytical approach to wave field reconstruction filtering in spatio-temporal frequency domain, IEEE Trans. Audio, Speech, Lang. Process., Vol. 21 (4), 2013,
pp. 685–696.

1107

[12] Koyama, S.; Furuya, K.; Hiwasaki, Y.; Haneda, Y.; Suzuki, Y. Wave field reconstruction filtering in
cylindrical harmonic domain for with-height recording and reproduction, IEEE/ACM Trans. Audio, Speech,
Lang. Process., Vol. 22 (10), 2014, pp. 1546–1557.
[13] Brix, S.; Sporer, T.; Plogsties, J. Carrouso - an European approach to 3D-audio, Proceedings of 110th
Conv. AES, Amsterdam, May 2001.
[14] Koyama, S.; Furuya, K.; Uematsu, H.; Hiwasaki, Y.; Haneda, Y. Real-time sound field transmission
system by using wave field reconstruction filter and its evaluation, IEICE Transactions on Fundamentals
of Electronics, Communications and Computer Sciences, Vol. E97-A (9), 2014, pp. 1840–1848.
[15] Spors, S.; Buchner, H. An approach to massive multichannel broadband feedforward active noise control
using wave-domain adaptive filtering, Proceedings of IEEE Int. Workshop Appl. Signal Process. Audio
Acoust. (WASPAA), New Paltz, Oct. 2007, pp. 171–174.
[16] Spors S.; Buchner, H.; Efficient massive multichannel active noise control using wave-domain adaptive
filtering, Proceedings of Int. Symp. Commun. Control Signal Process. (ISCCSP), Malta, Mar. 2008, pp.
1480–1485.
[17] Zhang, J.; Zhang, W.; Abhayapala, T. D. Noise cancellation over spatial regions using adaptive wave
domain processing, Proceedings of IEEE Int. Workshop Appl. Signal Process. Audio Acoust. (WASPAA),
New Paltz, Oct. 2015.
[18] Zhang, J.; Abhayapala, T. D.; Zhang, W.; Samarasinghe, N.; Jiang, S. Active noise control over space: A
wave domain approach, IEEE/ACM Trans. Audio, Speech, Lang. Process., Vol. 26 (4), 2018, pp. 774–786.
[19] Ueno, N.; Koyama, S.; Saruwatari, H. Kernel ridge regression with constraint of helmholtz equation forsound field interpolation, Proceedings of "IWAENC 2018; International Workshop on Acoustic Signal
Enhancement", Tokyo, Sept. 2018.
[20] Ueno, N.; Koyama, S.; Saruwatari, H. Sound field recording using distributed microphones based on
harmonic analysis of infinite order, IEEE Signal Process. Lett., Vol. 25 (1), 2018, pp. 135–139.
[21] Ito, H.; Koyama, S.; Ueno, N.; Saruwatari, H. feedforward spatial active noise control based on kernelinterpolation of sound field, Proceedings of IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP),
Brighton, May 2019 (to appear).
[22] Slock, D. T. M. On the convergence behaviour of the LMS and the normalized LMS algorithms, IEEE
Trans. Signal Process., Vol. 41 (9), 1993, pp. 2811–2825.
[23] Martin, P. A. Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles, Cambridge
University Press, New York, 1st edition, 2006.
[24] Anderson, C. A.; Teal, P. D.; Poletti, M. A. Spatial correlation of radial Gaussian and uniform sphericalvolume near-field source distributions, IEEE/ACM Trans. Audio, Speech, Lang. Process., Vol. 24 (1), 2016,
pp. 143–150.
[25] Ueno, N.; Koyama, S.; Saruwatari, H. Listening-area-informed sound field reproduction based on circularharmonic expansion, Proceedings of IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), New
Orleans, Mar. 2017, pp. 111–115.

1108