Optimized microphone deployment for near-field acoustic holography: To be, or not to be random, that is the question

(1)

Optimized microphone deployment for near-ﬁeld acoustic

holography: To be, or not to be random, that is the question

Mingsian R. Bai

, Jia-Hong Lin, Kwan-Liang Liu

Department of Mechanical Engineering, National Chiao-Tung University, 1001 Ta-Hsueh Road, Hsin-Chu 300, Taiwan

a r t i c l e

i n f o

Article history:

Received 3 February 2009 Received in revised form 26 January 2010 Accepted 26 January 2010 Handling Editor: L.G. Tham

a b s t r a c t

Arrays with sparse and random sensor deployment are known to be capable of delivering high quality far-field images without grating lobes. This raises the question of whether or not this idea can be applied to near-field imaging as well. To answer this question that has not yet been widely investigated in previous research, numerical simulations are undertaken in this paper to optimize the microphone deployment for both far-field and near-field arrays with the latter being the main focus. In the simulation, a recently introduced near-field equivalent source imaging (NESI) technique is employed for the near-field imaging. Global optimization techniques including the simulated annealing (SA) algorithm and the intra-block Monte Carlo (IBMC) algorithm are exploited to find the optimal microphone position efficiently. The combined use of the SA and the IBMC algorithms enables efficient search for satisfactory deployment with excellent beam pattern and relatively uniform distribution of microphones. In the near-field optimization, a special kind of beam pattern and cost function definition is used for the multiple-input-multiple-output (MIMO) imaging problem. As indicated by the simulation results, random deployment of microphones is necessary to avoid grating lobes in far-field imaging. In the near-field simulation, all results suggest that the optimal near-field array is the uniform rectangular array (URA) and the random deployment presents no particular benefit in near-field imaging.

1. Introduction

Array technology has been used in many diverse areas including radar [1], sonar [2], radio astronomy [3], telecommunications [4], and so forth. Its application encompasses purposes including signal enhancement, spatial filtering, direction of arrival (DOA) estimation, etc. Early development of arrays or beamformers was primarily based on the far-field assumption that the source is far away and the waves become spherical or planar at the array position. Far-field beamforming algorithms such as the delay-and-sum algorithm[5], the minimum variance distortionless response (MVDR) algorithm[6], and the multiple signal classification (MUSIC) algorithm[7]have been suggested in the past. Recently, array technology has found application in noise source identification (NSI) with the use of microphones. Microphone arrays serve as a powerful tool for acoustic field visualization that enables effective estimation of the positions and strengths of noise sources[5,8–10]. In comparison with the far-field arrays that are particularly useful for long-distance and large scale sources such as trains and aircrafts, near-field arrays represent a more recent and independent development of source imaging technology. With the advent of the near-field acoustical holography (NAH)[11,12], many near-field acoustic

Contents lists available atScienceDirect

journal homepage:www.elsevier.com/locate/jsvi

Journal of Sound and Vibration

_{Corresponding author.}

(2)

imaging techniques including the inverse boundary element method (IBEM)[13], the Helmholtz equation least squares (HELS)[14], the near-ﬁeld equivalent source imaging (NESI)[15], etc., are introduced. These near-ﬁeld techniques are well suited for imaging small-scale sources such as cars and computers by virtue of high resolution focusing schemes.

In array implementation, transducer deployment has been one of the key issues. It is well known that, for uniform linear arrays (ULA) and uniform rectangular arrays (URA)[16], array deployment must comply with the

l

/2-rule to avoid the spatial aliasing and the grating lobe problems[12]. Consequently, a large number of microphones are required to cover the source area, which can render the array conﬁguration impractical for sources at high frequencies. This prompts the development of non-uniform arrays that are capable of achieving high resolution and aliasing-free imaging with sparse sensors[17].

Given the fact that random deployment can be useful to far-ﬁeld arrays, a question arises naturally. Can this idea of random deployment be carried over to near-ﬁeld arrays? This interesting issue is relatively unexplored in the literature of the past. The

l

/2-rule—a well accepted criterion in deploying NAH arrays—can lead to undesirable high channel count. It is then tempting to ‘‘randomize’’ the sensor positions like in the case of far-ﬁeld arrays and achieve comparable performance with sparse deployment. To explore this conjecture, numerical simulation is undertaken in this study with the aid of optimization techniques for microphone deployment.

Several global optimization techniques are employed to deploy microphones for planar arrays used in far-field and near-field imaging. Monte Carlo (MC) simulation [18–22] is based on straightforward random search. Despite its simplicity, the MC method can be very time-onsuming. A more efficient technique, the simulated annealing (SA) algorithm

[23–25]is also used in the simulation. The SA algorithm relies on a search principle resembling the annealing process in the metallurgy. The search process follows an annealing schedule dictated by a temperature-dependent probability. The probability of accepting ‘‘worse’’ solutions in the initial high-temperature stage lends the SA method an effective approach for problems with many local minima[24]. Several researchers have applied the SA algorithm to optimize far-ﬁeld arrays

[26,27]. In this paper, modifications are made to enhance the search for optimal deployment. The intra-block Monte Carlo (IBMC) method conducts the random search only in the pre-partitioned local regions. This approach enhances search efficiency and often results in relatively uniform sensor deployment. A hybrid approach combining the SA and the IBMC methods is also presented to improve the search performance. Optimization of microphone deployment was carried out for both far-field and the near-field arrays. Near-field imaging is based on the results obtained using NESI and rectangular planar arrays. The main objective of the present study is to find the optimal deployment of the near-field array, while the far-field results are used as a means of testing the diverse optimization algorithms and contrasting the near-field results with the conventional wisdom in far-field imaging. Results are discussed and summarized in the conclusion section.

2. Near-ﬁeld imaging algorithm

In this paper, the recently introduced NESI algorithm[15] is employed for near-ﬁeld imaging. The idea of NESI is depicted inFig. 1, where the sound pressures picked up by the microphone array and the source strengths are related by

pM1¼GMJqJ1; (1)

where J and M denote the numbers of virtual sources and microphones, respectively, pM1 denotes the microphone pressure vector, GMJ denotes the propagation matrix and qJ1denotes the source strength vector. The source strengths

^

qJ1are estimated using the inverse ﬁlters CJMthat fulﬁll CJ MGM JEI, with I being a J J identity matrix, as

^ qJ1¼CJMpM1¼CJMGMJqJ1qJ1: (2) Lr Focal point L d Focal surface URA df Reconstruction Microphone qJ×1 pM×1 GM×J surface

(3)

This problem can be recast into a multiple-input-multiple-output (MIMO) model matching problem[15]. The ill-posedness involved in the inverse filter design must be properly dealt with using regularization methods such as singular value decomposition (SVD)[28]and Tikhonov regularization[29]. The resulting frequency response matrix of the inverse filters is converted to the time-domain filters. In what follows, however, we are content with the frequency-domain formulation in Eq. (3) in discussing array optimization.

Using the source strengths estimated by Eq. (2), we can reconstruct the sound pressure on a plane nearby[30], as shown inFig. 1, by

pr¼Grq; (3)

where Grdenotes the propagation matrix between the focal surface and the reconstruction surface. For more details of the NESI array, one may consult Ref.[15].

3. Beam patterns and cost functions

In this section, the beam patterns and the cost functions are defined to facilitate the following array optimization formulation. Since our focus is the near-field array, we only briefly review the far-field array as follows.

3.1. Far-ﬁeld beam pattern and cost function

For a far-ﬁeld array, the beam pattern can be deﬁned in the wave number domain[5]:

b ¼ 1 M

XM m ¼ 1

eikrm_; ₍₄₎

where k ¼

o

=c is the wave number,

o

is angular frequency, c is the speed of sound, and rmis the position vector of the mth microphone and k ¼ k

j

is the wave number vector of a plane wave incident from the direction represented by the unit vector

j

, as shown inFig. 2.

In optimizing far-ﬁeld performance, the aim is to minimize the maximum side-lobe level (MSL) of the beam pattern[5]. First, a circle with radius rmis drawn on the kxkyplane to deﬁne the scope of the main-lobe, which is a judicious choice

based on the beam pattern observations. The exterior of this circle is considered the side-lobe region. The cost function for far-ﬁeld arrays is deﬁned as

Q ¼m~

~s ; (5)

where ~m and ~s denote the maxima of the main-lobe and the side-lobes, respectively. Because ~m ¼ 1, maximizing the cost function Q amounts to minimizing the MSL.

3.2. Near-ﬁeld beam pattern and cost function

Due to the fact that the source is in the near-field and the focal points are multiple in numbers, the definition of near-field beam pattern is not as straightforward as in the far-field case. This paper suggests a procedure to calculate the

k = −kκ

Plane wave Source

r1 r2 rm

(4)

near-ﬁeld pressure ﬁeld beam pattern by scanning on the focal surface using a test point source. The procedure is depicted inFig. 3and described as follows:

(1) Design the inverse matrix C for the given array conﬁguration. The inverse matrix C is calculated using the propagation matrix G ¼ fGðxm;yj;

o

Þg. The entries in the propagation matrix G consist of the free-space Green’s functions:

Gðxm;yj;

o

Þ ¼ eikrjm

rjm ; (6)

where rjm¼ jxmyjjis the distance between the source point yjand the ﬁeld point xm.

(2) Position the test source at a grid point on the focal surface. Calculate the sound pressure vector p received at the microphones.

(3) Calculate the source strength vector at the focal points using the inverse matrix C:

q ¼ Cp (7)

(4) Calculate the pressure ﬁeld pz

rby propagating the spherical waves emitting from the point source at the zth grid point on the focal surface to the reconstruction surface by Eq. (3).

(5) Move the test source to next grid point on the focal surface and repeat steps (2)–(4). Superimpose the magnitude of the calculated pressure ﬁeld for all Z test source positions to get the following near-ﬁeld beam pattern vector:

b ¼ X Z z ¼ 1 jpz rð1; 1Þj XZ z ¼ 1 jpz rð2; 1Þj . . . XZ z ¼ 1 jpz rð

G

;1Þj " #T ; (8)

where

G

is the number of pixels on the reconstruction surface. Often, Zo

G

is selected.

With the beam pattern obtained above, a cost function is defined for the near-field array. Instead of the maximum main-lobe and side-lobe ratio as in the far-field array, the near-field cost function focuses on the number and positions of the main-lobe peaks. The desired focal points are located right underneath the microphones, whereas the synthesized beam pattern is calculated according to the procedure detailed above.

Test source

b Beam pattern vector C Inverse matrix Gr Propagation matrix G Propagation matrix z pr Reconstruction surface q Source strength p Microphone surface Focal surface

Move the test source to the next position

z + 1 z = Z No Yes (1,1) (2,1) ( ,1) Z _z r p p p Γ ∑ ∑ ∑ b = ... z=1 Z _z r z=1 Z _z r z=1

(5)

In the following, a near-field cost function is devised to ensure that the resulting beam pattern has well defined main-lobe at all focal points. First, the maximum of the beam pattern has been normalized to unity such that the regions with the value greater than 0.707 ( 3 dB) are taken as the main-lobes. Second, define the interior of the circle centered at each focal point with the radius rmas the scope of the main-lobe. The circle is used not only to define the main-lobe but also to check if a main-lobe is defocused. For the jth focal point, find the maximum of a main-lobe mjthat is greater than

0.707 ( 3 dB). The cost function of the near-ﬁeld array is deﬁned as

Q ¼ X J j ¼ 1

mj: (9)

The objective of the present array optimization is to ﬁnd the deployment that maximizes the cost function above.

4. Optimization algorithms

In this section, global optimization methods for microphone deployment are presented.

4.1. Intra-block Monte Carlo (IBMC) simulation

The basic MC algorithm is based on straightforward random search. For M microphones to be allocated to ðm þ 1Þ ðn þ 1Þ rectangular grid points, the number of possible combination is C_Mðm þ 1Þðn þ 1Þ, which is known to be an NP-complete problem [31]. Due to the blind search nature, the MC algorithm can be very inefﬁcient and result in non-uniform distribution of microphones that concentrate at certain areas. To address these problems, a modiﬁed method IBMC is proposed. By ‘‘intra-block’’, we mean the localized region designated to each microphone on the surface, as shown in

Fig. 4(a). The MC search is only conducted within each block with random positions generated inside this designated region. The M microphone elements will be designated to M localized regions. Hence, each region necessarily contains one and only one microphone. The ﬂowchart of IBMC is shown inFig. 5. Initially, m n divisions of a rectangular grid are set up on the microphone surface. Next, M localized search regions are designated to microphones, as shown inFig. 4(a). Each localized region inFig. 4(b) has the dimensions dmx¼dmy¼d, whereas the inter-element spacing of the grid points is chosen to be dx¼4d=m and dy¼5d=n, respectively. The localized regions are centered at the microphone positions of the uniform rectangular array (URA) that is selected to be the initial conﬁguration in the optimization. The associated data including the microphone positions xi, the beam pattern bi, and the cost function Qi are calculated. Next, each of the M microphone positions x is randomly assigned to one of the search points on the localized region. The new beam pattern b and the cost function Q are calculated for the assigned microphone positions x. The optimal solutions xopt, boptand Qoptare

Fig. 4. Schematic diagram of localized regions and grid points. (a) The localized regions (dashed lines) on the microphone surface with the inter-element spacing d =0.6 m. The symbol ‘‘&’’ indicates the microphone position. (b) The grid on a localized region. The microphone is constrained to move to one of nine grid points (including the central one) in the far-ﬁeld optimization. The microphone is assigned to any grid point in the localized region in the near-ﬁeld optimization for the SA search.

(6)

then replaced by the new solutions x, b and Q if Q 4 Qopt; otherwise the solutions are discarded. The simulation is continued until the number of iterations I exceed the preset value IIBMC.

The IBMC algorithm is more efﬁcient than the MC algorithm in that the search area for each microphone is far smaller. In addition, the IBMC algorithm generally results in microphone positions that are more uniformly distributed than those of the MC algorithm.

4.2. Simulated annealing (SA) technique

The MC algorithm can be very time-consuming and result in deployment that is far from optimal. Instead of blind search like the MC method, another efﬁcient SA algorithm is used in this study. SA is a generic probabilistic meta-algorithm for the global optimization problem, namely locating a good approximation to the global optimum of a given function in a large search space[23–26]. SA is well suited for solving problems with many local optima. Each point in the search space is analogous to the thermal state of the annealing process in metallurgy. At high temperatures, atoms with high internal energy are free to move to the other positions. As temperature drops, the internal energy is decreased to a lower state to gradually form a crystalline structure. The objective function Q to be maximized is likened to the internal energy in that state. One important feature of the SA approach is that it allows the search to move to a new state that is ‘‘worse’’ than the present one in the initial high-temperature stage. It is this mechanism that prevents the search from being trapped in a local maximum. The probability of accepting bad solutions decreases as temperature is decreased according to the Boltzmann distribution and the algorithm ﬁnally converges to the optimum solution.

Fig. 6illustrates the ﬂowchart of SA. For the problem of maximizing the array cost function, the array is initially set to be the URA with microphone positions xi. The corresponding beam pattern bi and cost function Qi are calculated. The

microphone surface is partitioned into m n divisions in a rectangular grid. The localized regions and the associated grid points are deﬁned in the same way as the IBMC. Accordingly, each microphone can be assigned to any position within the localized region in the simulation. The initial temperature Ti, the ﬁnal temperature Tf, and the annealing factor a are selected accordingly. A typical value of a is in the range of 0.8 and 0.99. Initially, set xopt¼xi, bopt¼biand Qopt¼Qi. Next, M microphone positions x are tentatively assigned. Each microphone is randomly assigned to one of the grid points with respect to the localized region. The beam pattern b and the cost function Q are evaluated for a new microphone positions x. Calculate the difference between the present and the optimal cost function,

D

Q ¼ Q Qopt: (10)

No No

Yes Yes

Set m×n grid on the microphone surface. Initial array = URA, microphone positions x_i , beam pattern bi , cost function Qi.

M localized regions are divided on the microphone surface.

xopt= xi, bopt= bi, Qopt= Qi

Randomly allocate M microphones to the grid

points in the localized regions.

Compute beam pattern b and the cost function Q

Q>Qopt Qopt=Q xopt=x bopt=b I + 1 I >I IBMC Stop

(7)

If T 4 Tf and

D

Q 4 0, replace the optimized solutions xopt, bopt and Qopt with the new solutions x, b and Q. Otherwise, if

D

Qr0, evaluate the following probability function:

Pð

D

Q ; TÞ ¼ eDQ =T_: ₍₁₁₎ The above probability will be compared with a random number 0r

g

r1 generated subject to the uniform distribution. A tentative solution is accepted when the probability function P is greater than the random number

g

; otherwise, the solution is rejected. Namely,

Pð

D

Q ; TÞ 4

g

; accepted Pð

D

Q ; TÞo

g

; rejected : (

(12)

Note that the larger the cost function difference

D

Q or the higher the temperature T, the higher is the probability to accept a worse solution.

As the search proceeds, the temperature is decreased according to an exponential annealing schedule that begins at some initial temperature T0and decreases the temperature in steps:

Tk þ 1¼a Tk; (13)

Choose initial temperature Ti, final temperature Tf , and

annealing coefficient a

Randomly allocate M microphones to nine grid points for

far-field arrays or to any grid points in the localized regions for near-field arrays

No No P > γ Yes No Yes Yes a×T

Compute beam pattern b and

the cost function Q

T ≤ Tf ΔQ> 0

Stop

Set m×n grid on the microphone surface. Initial array = URA, microphone positions x_i , beam pattern bi , cost function Qi.

M localized regions are divided on the microphone surface.

T = a×Ti, xopt= xi, bopt= bi, Qopt = Qi

Qopt= Q

xopt = x

bopt= b

(8)

where 0oao1 is the annealing coefﬁcient. The annealing process will be terminated if the temperature is lower than a preset ﬁnal temperature Tf. As the annealing process proceeds and T decreases, the probability of accepting a bad move

becomes increasingly small until it ﬁnally settles to a stable solution.

5. Numerical simulations

In this section, array optimization is conducted in the context of far-ﬁeld and near-ﬁeld imaging. Simulations with and without the intra-block (IB) constraint are carried out. The MC and SA algorithms are exploited to optimize microphone deployment with no IB constraint. On the other hand, the SA, IBMC and a combined SA-IBMC algorithm with 5 6 URA and random arrays with 30 microphones as the initial settings are employed to optimize microphone deployment with the IB constraint. Numerical simulations are all based on ten randomly generated starting points.

5.1. Optimizing far-ﬁeld array deployment by the MC and SA simulation

The MC and SA algorithm are exploited to optimize the microphone deployment. Initially, m n divisions (m =24 and n = 30) of a square grid are set up on the microphone surface, as shown inFig. 7(a). Each side of the square grid measures 0.1 m. The source frequency was considered as f= 1.7 kHz and the speed of sound as c =340 m/s, yielding the wave number k ¼ 2

p

f =c ¼ 31:4 m1_{. In addition, a URA of 5 6 (M= 30) deployment with inter-element spacing d= 0.6 m is used as a} benchmark in the following simulations, as shown inFig. 7(a). Its beam pattern calculated by Eq. (4) is shown inFig. 7(b). As expected, the grating lobes are clearly visible because the microphone spacing violates the

l

/2-rule (d= 3

l

at f= 1.7 kHz). The cost function calculated by Eq. (5) is only 1.0261 because of the grating lobes. This prompts the use of random deployment of microphones as follows.

In the MC simulation, the 30 microphones can freely occupy any 30 positions of the 25 31 grid points on array surface. Exhaustive search would require 16428144912combinations for a 30-element array, while only 105 _{iterations are}

carried out using this MC search. The MC search attains the optimal cost function 2.6532 at the 27596th iteration. The learning curve of the MC search is shown in Fig. 8(a). The corresponding microphone positions and beam pattern are shown inFigs. 8(b) and (c), respectively. Apart from the extremely time-consuming MC search, the SA approach is employed next. The annealing parameters of the SA for array deployment are chosen to be Ti= 10 1C, Tf= 1081C and a ¼ 0:95[23,25]. The learning curve of the SA search (405 iterations) is shown in the left portion (denoted as 1st SA) of

Fig. 8(d). The curve ﬂuctuates initially and then converges to a constant value 2.5767 between the 351st and the 405th iteration. The optimal microphone deployment and beam pattern are shown inFigs. 8(e) and (f). In addition to optimizing the microphone positions, optimizing the microphone weights can further improve the value of the cost function.

On the basis of the conﬁguration found previously by the SA, we continue to optimize the weights of microphones again using the SA algorithm. The number of iterations is increased to 1000. Starting from unity weights, the microphone weights are adjusted in each iteration with a random perturbation within the range of 0.1 to 0.1. The learning curve is shown in the right portion (denoted as 2nd SA) ofFig. 8(d). The cost function is further increased to 2.7561 at the 1283rd iteration. The resulting beam pattern is shown inFig. 8(g), where a unique main-lobe is clearly visible.

(9)

5.2. Optimizing far-ﬁeld array deployment by the SA and IBMC simulation

In this section, the SA, IBMC and a combined SA-IBMC algorithm are exploited to optimize microphone deployment with the IB constraint. Both microphone positions and weights are to be optimized using the SA algorithm. Speciﬁcally, the combined SA-IBMC method proceeds with three stages—the 1st SA stage, the IBMC stage, and the 2nd SA stage. The parameters of the two SA stages are identical to those in Section 5.1. The learning curve of the 1st SA stage (405 iterations) is shown in the left portion of Fig. 9(a). The curve ﬂuctuates initially and then converges to a constant value 2.5328 between the 208th and the 405th iteration. The resulting microphone deployment and beam pattern are shown in

Figs. 9(b) and (c). Being able to avoid local minima by accepting ‘‘bad’’ solutions in the initial SA search can be a beneﬁt and a shortcoming as well. A shortcoming of the SA algorithm is that it can miss the optimal solution in the initial stage and converges prematurely to a suboptimal one. A hybrid SA-IBMC approach is used in an attempt to address this problem.

The previous deployment obtained by the SA search is used as the input to the IBMC simulation. The microphone position can be randomly chosen from the nine grid points in the localized region. Each region necessarily contains one and only one microphone. Exhaustive search would require prohibitively 930combinations for a 30-element array, while only 100 iterations are required in the IBMC search. The learning curve of the IBMC (iteration 406–505) is shown inFig. 9(a). By the IBMC search, the cost function is further increased to 2.5465 at the 482nd iteration.Figs. 9(d) and (e) show the optimal

Fig. 8. The far-ﬁeld array optimized using the MC and the SA algorithms without the IB constraint at the frequency 1.7 kHz. For the MC simulation, maximum cost function Q = 2.6532 is attained at the 27596th iteration. For the SA simulation, maximum cost function Q= 2.7561 is attained at the 1283rd iteration. The circle indicates the main-lobe. (a) Learning curve of the MC search, (b) optimal array deployment obtained using the MC search, (c) beam pattern obtained using the MC search, (d) learning curve of the SA search, (e) optimal array deployment obtained using the SA search, (f) beam pattern obtained using the SA search, (g) beam pattern with weights optimized.

(10)

microphone positions and beam pattern obtained at the 482nd iteration. Next, in the 2nd SA stage, the microphone weights are optimized based on the conﬁguration found previously by the SA-IBMC approach. The microphone weights initially set to unity are adjusted in each iteration with a random perturbation within the range of 0.1 to 0.1. The learning curve in 506 iterations is shown inFig. 9(a). The cost function is further increased to 2.6602 at the 1429th iteration. The resulting beam pattern is shown inFig. 9(f), where a unique main-lobe is clearly visible.

Apart from the URA, the random array deployment is also used as the initial setting in the simulation. For brevity, the results of MC, IBMC, SA and SA-IBMC simulations are summarized in Table 1. The highest value of the cost function obtained with these ten randomly generated starting points is regarded as the optimal Q inTable 1. The simulation results obtained with and without the IB constraint are compared in terms of number of iterations and the maximum cost function values. Although the MC approach has reached the highest cost function (Q= 2.6532), it takes 27 596 iterations to achieve this value. By comparing the results of the MC and IBMC (with the URA as the initial setting), we found that the IBMC approach can attain comparable cost function value to the MC approach with far less amount of computation (Q=2.5638 at the 7662nd iteration of IBMC vs. Q= 2.6532 at the 27 596th iteration of MC). In comparison with the results obtained using the SA algorithm with the IB constraint (Q=2.6602 for the URA as the initial setting and Q=2.6573 for a random array as the initial setting), the SA approach with no IB constraint has attained a slightly higher cost function (Q= 2.7561) with comparable computational complexity. It all boils down to the tradeoff between search time and optimality.

Incorporating the IB constraint could potentially have the following beneﬁts. First, the IBMC algorithm is computationally more efﬁcient than the plain MC algorithm because of smaller search areas. Second, in the hybrid SA-IBMC approach, the IB constraint could possibly improve the SA results when the SA algorithm converges prematurely to a suboptimal result. Third, the IB constraint normally results in uniform distributions of microphones. By ‘‘uniform’’, we simply mean that microphones would not concentrate at only a few areas, which should not be confused with the

(11)

Fig. 9. The far-ﬁeld array optimized using the MC algorithm and the combined SA-IBMC algorithm with the IB constraint at the frequency 1.7 kHz. Maximum cost function Q= 2.6602 is attained at the 1429th iteration. The circle indicates the main-lobe. (a) Learning curve, (b) optimal array deployment obtained using the SA search, (c) beam pattern obtained using the SA search, (d) optimal array deployment obtained using the SA-IBMC search, (e) beam pattern obtained using the SA-IBMC search, (f) beam pattern with weights optimized.

(12)

deployment of the constant-spacing uniform arrays. In summary, it is fair to say that the IB constraint signiﬁcantly reduces the computation complexity at the risk of converging to a suboptimal solution which may not be far from the global optimum. This is generally sufﬁcient in practical applications.

Apart from the source frequency of 1.7 kHz, we also run the simulation for the frequencies of 500 Hz and 1 kHz. For brevity, we only summarize the results inTable 2. Random arrays yield unique main-lobe and higher cost function than the URA at 1 kHz. For the lower frequency of 500 Hz, no grating lobes are seen in the beam pattern of URA, while a higher side-lobe level is found in the beam pattern of the random array. This leads to a higher value of cost function for the URA than the random array at low frequencies. In the next section, we will examine whether the same idea of random array applies to near-ﬁeld imaging as well.

5.3. Optimizing near-ﬁeld array deployment

In this simulation, the MC, SA, IBMC algorithms and a combined SA-IBMC algorithm are exploited to optimize microphone deployment for near-ﬁeld imaging. The 5 6 URA and random array with 30 microphones are used as the initial setting for the simulation of near-ﬁeld imaging. With reference to Fig. 1, the array parameters are chosen as: microphone spacing d =0.6 m, focal point spacing df= 0.6 m (3

l

at the source frequency of 1.7 kHz), and distance of

reconstruction L=0.1 m. The focal points are collocated underneath with the microphones on a conformal plane (the number of focal points is equal to the number of microphones). Thus, ﬁner 25 31 grid points are set up on the microphone surface. The microphone is free to move to any grid point inside the assigned localized regions. The dimensions of the localized regions (see Fig. 4(b)) are dmx= dmy=0.6 m and the spacing of the grid points is

dx=dy= 0.1 m. The main-lobes are deﬁned as the interior of the circles centered at each focal point with the radius

rm= 0.02 m (the minimal distance between two adjacent grid points). The iteration limit IMCand IIBMCare chosen to be 5000

and 1000, respectively.Figs. 10(a) and (b) show the cost function history of the MC search and the IBMC search with the URA as the initial setting, respectively. The maximum value of the cost function Q= 27.3 corresponds to the initial URA deployment cost function at the frequency of 1.7 kHz. Apart from the time-consuming MC algorithm and the IBMC algorithm, another attempt was made to find the optimal deployment using the more efficient SA algorithm. The annealing parameters are chosen to be: Ti= 10 1C, Tf=1081C and a ¼ 0:95. After 405 iterations, we obtained the learning curve with the URA as the initial setting, as shown inFig. 10(c). This learning curve converges to a stable value 3.8502 at the low-temperature stage, as a typical SA behavior. This Q value is apparently suboptimal since it is smaller than the initial Q= 27.3 of the URA. The maximum cost function remains the initial cost function (Q= 27.3 at the frequency 1.7 kHz) corresponding to the URA deployment. Finally, the hybrid approach SA-IBMC is employed to optimize the near-field array with the URA as the initial setting.Fig. 10(d) shows the learning curve (405 SA iterations+ 100 IBMC iterations). The maximal cost function

Table 1

The search performance of different optimization methods for far-ﬁeld array deployment with the inter-element spacing d =0.6 m.

Constraint Method Find best Q iterations Best Q (linear)

Without IB (initially random array) MC 27596 2.6532

SA 351 2.5767

SA + w 1283 2.7561

With IB (initially URA) IBMC 7662 2.5638

SA 208 2.5328

SA + IBMC 482 2.5465

SA + IBMC+ w 1429 2.6602

With IB (initially random array) IBMC 23285 2.5617

SA 222 2.5224

SA + IBMC 406 2.5224

SA + IBMC+ w 1352 2.6573

The letter ‘‘w’’ indicates that weight optimization is performed.

Table 2

The comparison of converged cost function Q of the URA and the optimized far-ﬁeld random arrays at three different frequencies.

Array f = 500 Hz f =1 kHz f = 1.7 kHz

URA 4.0216 1.0192 1.0261

Random array (without IB, initially random array) 1.5961 2.5451 2.7561

Random array (with IB, initially URA) 2.5048 2.3324 2.6602

(13)

remains to be that of the initial URA deployment.Table 3summarizes the search performance of near-field optimization methods for three different frequencies 500 Hz, 1 kHz and 1.7 kHz with the URA and the random array as the initial settings. As concluded from the table, even though the other initial setting is used for optimizing the random array, the optimal deployment with the maximal cost function remains to be the URA. These results suggest that the optimal near-field array is the URA when the number of focal points is equal to the number of microphones. The near-field beam patterns defined in Section 3 calculated for the URA inFig. 7(a) at the frequencies 500 Hz, 1 kHz and 1.7 kHz are shown in

Figs. 11(a)–(c). Thirty main-lobes with equal height centered at the focal points are clearly visible. The beam pattern of a non-optimized near-field array at the frequency 1.7 kHz by the SA and SA-IBMC algorithms is shown inFig. 11(d). Clearly, the converged pattern is a dramatic departure from the URA and the associated beam pattern is quite ugly with quite a few ‘‘defocused’’ points. To see the contrast, we also show an example of random deployment that corresponds to the optimized far-field deployment obtained previous using the SA-IBMC method (the near-field cost function).Figs. 12(a)–(c) show the beam patterns of the deployments at the frequencies 500 Hz, 1 kHz and 1.7 kHz. It can be seen from the quite smeared beam patter that many main-lobes are defocused or even missed. These observations lead to a conclusion that contradicts the far-field experience—random deployment degrades the multi-focusing performance and the URA is the optimal array. In addition,Table 4summarizes the performance of different numbers of focal points using the same 5 6 URA with the inter-element spacing d= 0.6 m (corresponding to 3

l

for the frequency 1.7 kHz) at the frequencies 500 Hz,

Fig. 10. The cost function history of near-ﬁeld array optimization with inter-element spacing 0.6 m (3lat the frequency 1.7 kHz). Maximum cost function value is Q= 27.3 for the initial URA deployment. (a) All Q values, including those rejected during the random search by the MC algorithm, (b) all Q values, including those rejected during the random search by the IBMC algorithm, (c) only the legitimate values accepted by the SA algorithm, (d) only the legitimate values accepted by the SA-IBMC algorithm.

(14)

Table 3

The search performance of different optimization methods for near-ﬁeld array deployment with the inter-element spacing d = 0.6 m at three different frequencies.

Method Frequency (Hz) Best Q (linear)

Initially URA Initially random array

MC 500 28.2 8.5 1000 26.9 8.3 1700 27.3 8.2 IBMC 500 28.2 9.1 1000 26.9 9.5 1700 27.3 10 SA 500 28.2 4.2 1000 26.9 5.6 1700 27.3 4.7 SA +IBMC 500 28.2 4.2 1000 26.9 5.6 1700 27.3 4.7

Fig. 11. The beam pattern of the optimized near-ﬁeld array (URA). The inter-element spacing is 0.6 m (3lat the frequency 1.7 kHz). (a) At the frequency 500 Hz, Q= 28.2, (b) at the frequency 1 kHz, Q = 26.9, (c) at the frequency 1.7 kHz, Q =27.3, (d) the beam pattern of a non-optimized near-ﬁeld array at the frequency 1.7 kHz by the SA and SA-IBMC algorithms, Q= 3.8502. The circles indicate the main-lobe. The symbol ‘‘ ’’ indicates the focal points.

(15)

1 kHz and 1.7 kHz. The cost function for 5 6 focal points (the number of microphones) is 27.3 at the frequency 1.7 kHz. It increases to 28.3 when the number of focal points is increased to 9 11. However, for 13 16 focal points the cost function decreases to 16.6, showing that it does not increase indeﬁnitely with the number of microphones.

6. Conclusions

Optimized planar array deployment for source imaging is investigated in this paper. The outcome of the work is twofold. First, global optimization algorithms have been developed to facilitate the search of the optimized microphone

Fig. 12. The beam pattern of a non-optimized near-ﬁeld array. The microphone deployment is identical to the optimized far-ﬁeld array obtained using the SA-IBMC algorithm. (a) At the frequency 500 Hz, Q = 3.5, (b) at the frequency 1 kHz, Q= 3.4, (c) at the frequency 1.7 kHz, Q= 3.9. The circles indicate the main-lobe. The symbol ‘‘ ’’ indicates the focal points.

Table 4

The performance obtained using different numbers of focal points in a 5 6 URA with the inter-element spacing d = 0.6 m at three different frequencies.

Frequency (Hz) Q (linear)

5 6 focal points 9 11 focal points 13 16 focal points

500 28.2 28.8 26.7

1000 26.9 27.2 21.8

(16)

deployment. Second, whether or not randomization is necessary in near-field imaging is explored with the aid of these optimization techniques. The SA algorithm and the combined SA-IBMC algorithm proved effective in finding the optimal deployment for far-field imaging. For far-field array with sparse deployment in which inter-element spacing is large, random deployment with optimal weights is crucial to avoiding grating lobes. As predicted by the conventional wisdom, the optimized random sparse array provides excellent beam pattern with a unique main-lobe.

To answer the question whether the idea of random deployment acquired in far-field imaging applies to near-field imaging as well, simulation employing the MC, IBMC, SA and SA-IBMC optimization techniques were conducted. Due to its complicated near-field and multi-focusing nature, a special kind of beam pattern alongside the cost function is defined. It was observed that all optimization approaches lead to suboptimal solutions with Q values far less than that of the URA. The maximum cost function is always the initial cost function corresponding to the URA deployment. All results of simulations suggest that the optimal near-field array is the URA and random deployment presents no particular benefit in near-field imaging.

As a limitation of the present study, the conclusion above concerning near-ﬁeld imaging is based on the results obtained using rectangular planar arrays. Whether the same conclusion applies to other array geometries remains to be explored in the future.

Acknowledgment

The work was supported by the National Science Council in Taiwan under the project number NSC 96-2628-E-009-028-MY3.

References

[1] M.I. Skolnik, Introduction to Radar Systems, McGraw-Hill, New York, 1980.

[2] A.V. Oppenheim, R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1989. [3] S. Haykin, Array Signal Processing, Prentice Hall, Englewood Cliffs, NJ, 1985.

[4] H.T. Friis, C.B. Feldman, A multiple unit steerable antenna for short-wave reception, Bell System Technical Journal 16 (1937) 337–419. [5] J.J. Christensen, J. Hald, Beamforming, Br ¨uel & Kjær Technical Review No. 1, 2004.

[6] J. Capon, High-resolution frequency-wavenumber spectrum analysis, Proceedings of the IEEE 57 (1969) 1408–1418.

[7] R.O. Schmidt, Multiple emitter location and signal parameter estimation, IEEE Transactions on Antennas and Propagation AP-34 (3) (1986) 276–280. [8] M.R. Bai, J.W. Lee, Industrial noise identiﬁcation by using an acoustic beamforming system, ASME, Journal of Vibration and Acoustics 120 (1998)

426–433.

[9] D. Zha, T. Qiu, Underwater sources location in non-Gaussian impulsive noise environments, Digital Signal Processing 16 (2006) 149–163. [10] U. Benko, J. Petrovcic, D. Juricic, J. Tavcar, J. Rejec, A. Stefanovska, Fault diagnosis of a vacuum cleaner motor by means of sound analysis, Journal of

Sound and Vibration 276 (2004) 781–806.

[11] E.G. Williams, J.D. Maynard, Holographic imaging without the wavelength limit, Physical Review Letters 45 (1980) 554–557.

[12] E.G. Williams, J.D. Maynard, E. Skuarzyk, Sound source reconstruction using a microphone array, Journal of the Acoustical Society of America 68 (1980) 340–344.

[13] M.R. Bai, Application of BEM (boundary element method)-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries, Journal of the Acoustical Society of America 92 (1992) 533–549.

[14] S.F. Wu, On reconstruction of acoustic pressure ﬁelds using the Helmholtz equation least squares method, Journal of the Acoustical Society of America 107 (5) (2000) 2511–2522.

[15] M.R. Bai, J.H. Lin, Source identiﬁcation system based on the time-domain nearﬁeld equivalence source imaging: fundamental theory and implementation, Journal of Sound and Vibration 307 (2007) 202–225.

[16] P.A. Nelson, S.H. Yoon, Estimation of acoustic source strength by inverse methods: part II, experimental investigation of methods for choosing regularization parameters, Journal of Sound and Vibration 233 (2000) 669–705.

[17] S. Christopher, Numerical annealing of low-redundancy linear arrays, IEEE Transactions on Antenna and Propagation 41 (1) (1993) 85–90. [18] N. Metropolis, S. Ulam, The Monte Carlo method, Journal of the American Statistical Association 44 (1949) 335–341.

[19] J.M. Hammersley, D.C. Handscomb, Monte Carlo Methods, Chapman and Hall, London, 1979.

[20] F. Dellaert, D. Fox, W. Burgard, S. Thrun, Monte Carlo localization for mobile robots, IEEE International Conference on Robotics and Automation, An Arbor, 1999, pp. 1322–1328.

[21] A. Baggio, K. Langendoen, Monte Carlo localization for mobile wireless sensor networks, Ad Hoc Networks 6 (2008) 718–733.

[22] G. Marmidis, S. Lazarou, E. Pyrgioti, Optimal placement of wind turbines in a wind park using Monte Carlo simulation, Renewable Energy 33 (2008) 1455–1460.

[23] S. Kirkpatrick, C.D. Gelatt Jr., M.P. Vecchi, Optimization by simulated annealing, Science 220 (4598) (1983) 671–680.

[24] R.T. Eduardo, J.K. Hao, T.J. Jose, An effective two-stage simulated annealing algorithm for the minimum linear arrangement problem, Computer and Operation Research 35 (2008) 3331–3346.

[25] B. Franco, Simulated annealing overview. /http://www.geocities.com/francorbusetti/saweb.pdfS, 2003, last viewed on June 22, 2009.

[26] V. Murino, A. Trucco, C.S. Regazzoni, Synthesis of unequally spaced arrays by simulated annealing, IEEE Transactions on Signal Processing 44 (1) (1996) 119–123.

[27] A. Trucco, Thinning and weighting of large planar arrays by simulated annealing, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 46 (2) (1999) 347–355.

[28] B. Nobel, J.W. Daniel, Applied Linear Algebra, Prentice-Hall, Upper Saddle River, NJ, 2000.

[29] O. Kirkeby, P.A. Nelson, H. Hamada, Fast deconvolution of multichannel systems using regularization, IEEE Transactions on Speech Audio Processing 6 (1998) 189–194.

[30] P.M. Morse, K.U. Ingard, Theoretical Acoustics, Princeton University Press, Princeton, NJ, 1986.