Alix Jeannerot
email: alix.jeannerot@inria.fr

Malcolm Egan
email: malcom.egan@inria.fr

Lélio Chetot
email: lelio.chetot@inria.fr

Jean-Marie Gorce
email: jean-marie.gorce@inria.fr

Mitigating User Identification Errors in Resource Optimization for Grant-Free Random Access

Keywords: random access, resource allocation, user identification, stochastic gradient descent

In grant-free random access, a key question is how devices should utilize resources without coordination. One standard solution to this problem are strategies where devices randomly select time-slots based on an optimized stochastic allocation rule. However, the optimization of this allocation rule requires accurate knowledge of which devices have been active in previous frames. As user identification algorithms are subject to errors, the expected throughput of the optimized allocation can be highly suboptimal. In this paper, we propose algorithms for optimization of device time-slot allocations that mitigate the impact of user identification errors. We show that when the activity distribution with and without errors is known, then our algorithm converges with probability one to a stationary point. When the activity distributions are not available, we introduce new theoretically-motivated heuristics which significantly improve the expected throughput over existing algorithms and approach the performance when errors are not present.

I. INTRODUCTION

Resource allocation is a central problem for wireless access networks, particularly for the massive access systems critical for 6G [START_REF] Lee | 6G massive radio access networks: Key applications, requirements and challenges[END_REF], where a large number of devices aim to transmit data to an access point. In 5G-NR and NB-IoT, frames are decomposed into time-frequency resources. A common problem is therefore to determine which devices should be allocated to a given time-slot within a frame.

In order to avoid high levels of control signaling in massive access systems, it is desirable to provide devices with a timeslot selection policy [START_REF] Paolini | Coded slotted ALOHA: a graphbased method for uncoordinated massive access[END_REF], [START_REF] Kalør | Random access schemes in wireless systems with correlated user activity[END_REF], [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF]. This means that when a device has data to transmit, limited control signaling is required in order to determine which resources the device should use. The selection of a policy, often a variant of ALOHA (e.g., [START_REF] Paolini | Coded slotted ALOHA: a graphbased method for uncoordinated massive access[END_REF], [START_REF] Wieselthier | An exact analysis and performance evaluation of framed ALOHA with capture[END_REF]), is typically based on a criterion such as the expected system throughput, outage probability, or expected sum-rate.

In event-driven IoT systems, devices (e.g., sensors) only require data transmission at irregular frames. The expected throughput therefore depends on the activity distribution, which governs the probability a given device will be active in a frame and can vary device to device. In general, the activity distribution may not be perfectly known, due to imperfect estimation of active devices within each frame. Indeed, the estimation of active devices, also known as user identification, is currently under active investigation [START_REF] Ahn | EP-based joint active user detection and channel estimation for massive machine-type communications[END_REF], [START_REF] Ke | Compressive sensingbased adaptive active user detection and channel estimation: Massive access meets massive MIMO[END_REF].

In this paper, we consider the impact of imperfect estimation of active devices on resource allocation based on the expected throughput. Our focus is on the problem of designing strategies for allocating a single time slot within a frame to an active device [START_REF] Kalør | Random access schemes in wireless systems with correlated user activity[END_REF], [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF], relevant for NB-IoT. In this setting, a key challenge is to account for uncertainties in user activity probabilities, which are ignored in most existing work [START_REF] Paolini | Coded slotted ALOHA: a graphbased method for uncoordinated massive access[END_REF], [START_REF] Kalør | Random access schemes in wireless systems with correlated user activity[END_REF], [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF] on resource allocation and coding for random access. Our main contributions are as follows:

• We demonstrate that the stochastic resource optimization algorithms recently proposed in [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF] can be highly susceptible to user identification errors, leading to suboptimal performance. This is significant as in the absence of user identification errors, these algorithms are optimal in contrast with other heuristic approaches [START_REF] Kalør | Random access schemes in wireless systems with correlated user activity[END_REF] and do not assume perfect knowledge of activity probabilities as in [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF]. • To mitigate the impact of the bias arising from user identification errors, we use importance sampling in a manner analogous to sample selection bias correction [START_REF] Zadrozny | Learning and evaluating classifiers under sample selection bias[END_REF] and bias reduction in private synthetic data [START_REF] Ghalebikesabi | Bias mitigated learning from differentially private synthetic data: A cautionary tale[END_REF]. We prove that by weighting the gradient estimates in the algorithms in [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF], the resource allocation converges to a stationary point with probability one. • In practice, the weight for the gradient estimates requires the evaluation of the true activity distribution and the imperfect activity distribution induced by erroneous user identification. As these distributions may be difficult to obtain, we propose two heuristic weights that are readily available in practical systems, leading to new stochastic resource optimization algorithms. • The proposed algorithms are validated via a numerical study. We show user identification errors can significantly decrease the expected throughput. Realistic scenarios where a low-complexity message passing algorithm [START_REF] Chetot | Joint identification and channel estimation for fault detection in industrial IoT with correlated sensors[END_REF] is utilized to estimate active devices are also presented.

In certain scenarios it is even possible to approach the performance achieved without the presence of user identification errors.

II. SYSTEM MODEL

Consider a time-slotted network consisting of an access point (AP) and N devices equipped with a single antenna and utilizing a common frequency band. Transmissions occur in frames of a fixed duration, where each device n is active with probability p n in each frame. The activity of each device is denoted by X n ∼ Ber(p n ) and the joint activity probability mass function by p X . We assume that the activities of the devices are mutually independent. As such,

p X (x) = Π N n=1 p xn n (1 -p n ) 1-xn , x ∈ {0, 1} N . (1)
In frame t, the activity of all devices is denoted by X t = [X t 1 , . . . , X t N ] T . The activity vector X t is also assumed to be mutually independent of X t ′ , t ′ ̸ = t.

A. Transmission Protocol

The transmission protocol is detailed in Alg. 1. Each frame consists of a preamble and K data slots. The resource allocation is communicated to devices every f frames to limit downlink usage. Sync and Preamble: After receiving the sync signal in Step 1, each active device in the current frame transmits a pilot sequence of length L during the preamble phase in Step 2, denoted by q n ∈ C L . At the end of the preamble, the AP attempts to identify active devices, as detailed in Step 3.

In particular, the received signal y t ∈ C L during the preamble phase is given by:

y t = N n=1 X t n h t n q n + w t 0 (2)
Where h t n ∈ C is the fading coefficient of the n-th device and w t 0 ∈ C L is circular complex Gaussian noise with variance τ 2 w . An efficient method to estimate the set of active devices in Step 3 is provided by approximate message passing algorithms (AMP) within the hybrid-gaussian AMP (HGAMP) family [START_REF] Chetot | Joint identification and channel estimation for fault detection in industrial IoT with correlated sensors[END_REF], [START_REF] Zou | Message passing based joint channel and user activity estimation for uplink grant-free massive MIMO systems with low-precision ADCs[END_REF]. Nevertheless, identification of active devices is not perfect. Denoting the estimated activity vector at time t by Xt , the probability of observing the vector x at the output of the HGAMP algorithm is denoted by Pr( Xt = x).

Data Transmission: Each active device transmits in exactly one of the K data slots. In Step 5, active devices select a slot randomly, governed by an allocation matrix A t with A t ij , i ∈ {1, . . . , N }, j ∈ {1, . . . , K} corresponding to the probability of device i selecting slot j, with k A t nk = 1, n = 1, . . . , N (as each active user transmits exactly once in a frame).

B. Optimization Objective

The main focus of this paper is the optimization of the slot selection matrix A. In the following sections, we drop the notation t when the frame index is not relevent. This problem has recently been considered in [START_REF] Kalør | Random access schemes in wireless systems with correlated user activity[END_REF], [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF] when the activity vector is perfectly known in each frame. We consider optimization of the expected throughput, defined by [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF]:

T (A) = E X [T X (A)], (3)
where T X (A) is defined as

T n (A; x) = K k=1 X n A nk N m=1 m̸ =n (1 -X m A mk ) T X (A) = N n=1 T n (A; X)
The objective in (3) can be interpreted as the expected number of slots in which devices can transmit collision-free.

The resource allocation problem is then to find an allocation matrix A that maximizes T (A). In particular, we seek a solution A * such that

A * = arg max A∈R N ×K + ; k A n,k =1 T (A)
Note that the objective T (A) is non-convex in A. We also highlight that the methods presented in the remaining of this paper are not limited to the maximization of the expected throughput T (A) and are applicable to other objectives; e.g., proportional-fair, expected sum-rate and outage probability [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF].

III. THE IMPACT OF IMPERFECT USER IDENTIFICATION

As observed in [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF], [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF], the problem in (3) can be viewed as a stochastic optimization problem. Under the assumption that the activity vectors X 1 , X 2 , . . . are perfectly known to the AP, a natural solution is stochastic gradient ascent (SGA), which is detailed in Alg. 2. Here,

H = {A ∈ R N ×K + : K k=1 A n,k = 1, n = 1, . . . , N } (4)
and Π H [•] denotes the projection to the closest point in H with respect to the Euclidean norm. The matrix of the gradient estimate g(A; X) consists of elements g(A; X) ql , q ∈ {1, . . . , N }, l ∈ {1, . . . , K} defined by [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF]:

g(A; X) ql = X q N m=1 m̸ =q (1 -X m A ml ) - N n=1 n̸ =q X q X n A nl N m=1 m̸ =n m̸ =q (1 -X m A ml ). (5)
Note that in this case, g(A; X) is an unbiased estimate of ∇T (A) = E X [∇ A T X (A)] due to the absence of user identification errors.

Algorithm 2: Stochastic optimization algorithm when user identification is error-free. 1 Choose initial allocation matrix A 1 ∈ H and step-size sequence {α t } with α t > 0, t = 1, 2, . . .

2 t ← 1. 3 while not converged do 4 Based on X t , compute an unbiased estimate g(A t ; X t ) of ∇ A t T (A t ) 5 A t+1 ← Π H [A t + α t g(A t ; X t )] 6 t ← t + 1 7 end
When the activity vectors are perfectly known, under weak conditions on the step-size sequence (detailed in Theorem 1), the SGA algorithm converges almost surely to a stationary point [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF]. In fact, this convergence result holds even when device activity is correlated and the process X 1 , X 2 , . . . is an irreducible Markov process [5, Theorem 1].

However, with realistic user identification algorithms, errors will sometimes be introduced in the estimation of the activity X t . As shown in Example III.1, the errors introduced by imperfect user identification can strongly affect the network throughput.

Exemple III.1. To illustrate the potential impact of user identification errors on resource allocation, consider the simple network consisting of 4 devices sharing two time slots. The true activity probability of the devices is p = 0.2 0.3 0.6 0.9 . Suppose that, due to non-orthogonal pilots or a very low signal-to-noise ratio (SNR), the user identification algorithm might not be able to differentiate between the first and last users. This means that, with probability ϵ, the user activity vector is in fact drawn from p = 0.9 0.3 0.6 0.2 . Fig. 1 shows the impact of ϵ on the network throughput T (A) after optimizing A. The case ϵ = 0 corresponds to the case of perfect user identification and ϵ = 1 to the case where an error event always occurs. Observe that the errors in the user identification can strongly affect the throughput, especially if error events are quite likely to occur (i.e., ϵ ≈ 1). We also compare to the throughput of a Frame Slotted Aloha [START_REF] Wieselthier | An exact analysis and performance evaluation of framed ALOHA with capture[END_REF] network, i.e., a network that always use a constant matrix of elements A i,j = 1 K , ∀i, j. and lead to a throughput worse than an Aloha network. The blue line corresponds to a perfect identification (ϵ = 0).

IV. A NEW ALGORITHM TO MITIGATE ERRORS

As illustrated in Fig. 1, user identification errors can have a significant impact on resource optimization. In this section, we propose new algorithms in order to mitigate the impact of user identification errors on optimization of the allocation matrix A.

A. The Effect of Errors

In order to develop a new algorithm to mitigate the impact of user identification errors, we first examine the impact of the errors on the SGA algorithm in Alg. 2. To this end, consider the update rule of Alg. 2:

A t+1 = Π H [A t + α t g(A t ; X t )], (6)
Observe that the update rule can be rewritten as

A t+1 = A t + α t g(A t ; X t ) + α t Z t ,
where α t Z t is the projection term which represent the smallest vector (w.r.t. the Euclidean norm) needed to project

A t + α t g(A t ; X t ) into the constraint set H, g(A t ; X t ) is the gradient estimate, E X [g(A t ; X t )] = ∇ A t T (A t
) + β t with β t the bias term. We then have

A t+1 = A t + α t E X [g(A t ; X t )] + g(A t ; X t ) -E X [g(A t ; X t )] =∂M t + α t Z t (7) = A t + α t ∇ A t T (A t ) + β t + ∂M t + α t Z t (8)
with E[∂M t ] = 0. In the case when there is no error in the user identification, g(A t ; X t ) is an unbiased estimate of the gradient ∇ A T (A t ), and hence β t = 0. This unbiased property is critical in order to guarantee convergence of Alg. 2 [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF].

On the other hand, when user identification errors are present, β t ̸ = 0. As a consequence g(A t ; Xt ) is a biased estimate of ∇ A t T (A t ). Hence, there is no guarantee of convergence to a stationary point of (3).

B. A Bias Reduction Method

We have observed that the key cause of the performance reduction in Fig. 1 is the bias in the gradient estimates introduced by imperfect user identification. In order to improve the expected throughput, it is therefore desirable to reduce the bias in g(A t ; X t ).

We propose the following solution, inspired by importance sampling [START_REF] Ghalebikesabi | Bias mitigated learning from differentially private synthetic data: A cautionary tale[END_REF], [START_REF] Zadrozny | Learning and evaluating classifiers under sample selection bias[END_REF], to reduce the bias. Consider the weight function w : x → w(x) for x ∈ {0, 1} N defined by

w(x) = Pr(X = x) Pr( X = x) , ( 9
)
where X is the true activity vector and X is an estimate; e.g., obtained from HGAMP.

A key property of the importance weight w(x) is that for all

A ∈ R N ×M in H, E X [g(A; X)] = x g(A; x)Pr(X = x) = x g(A; x) Pr(X = x) Pr( X = x) Pr( X = x) = E X[w( X)g(A; X)]. (10)
In other words, w( X)g(A; X) is an unbiased estimate of g(A; X), as long as w(x) < ∞ for all x such that Pr(X = x) > 0. As such, as we will rigorously establish shortly, this choice of weight overcomes the key problem preventing convergence of Alg. 2.

Incorporating the weight defined in (9) into the SGA algorithm yields Alg. 3. This algorithm has the following convergence guarantee.

t ; Xt ) of ∇ A t T (A t ) 5 A t+1 = Π H [A t + α t w( Xt )g(A t ; Xt )] 6 t ← t + 1 7 end
Theorem 1. The iterates A t of Algorithm 3, converge almost surely as t → ∞ to a stationary point provided that the step-size sequence {α t } with α t > 0, t = 1, 2, . . . has the following properties:

∞ t=1 α t = ∞, ∞ t=1 (α t ) 2 < ∞
and w(x) < ∞ for all x such that Pr(X = x) > 0.

Proof. Similar to [START_REF] Ahn | EP-based joint active user detection and channel estimation for massive machine-type communications[END_REF], the update rule can be written as

A t+1 = A t + α t E Xt [w( Xt )g(A t ; Xt )] + w( Xt )g(A t ; Xt ) -E Xt [w( Xt )g(A t ; Xt )] =∂M ′t + α t Z t = A t + α t E X t [g(A t ; X t )] + ∂M ′t + α t Z t (11) = A t + α t ∇ A t T (A t ) + β t + ∂M ′t + α t Z t = A t + α t ∇ A t T (A t ) + ∂M ′t + α t Z t ( 12
)
where the step ( 11) is obtained through property:

E X [g(A; X)] = E X[w( X)g(A; X)]. (13)
As such, β t = 0 as the property in [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF] implies that E X [g(A t ; X)] is an unbiased estimate of ∇ A t T (A t ).

To establish convergence, we seek to apply Theorem 5.2.1 of [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]. To do so, it is necessary to check the following assumptions:

• A.5.2.1:

sup t E Xt [w( Xt )g(A t ; Xt )] < ∞, (14)
which holds provided that P( Xt = x t ) ̸ = 0 for all x t . • A.5.2.2: There exists g : R N ×M → R such that E X[w( Xt )g(A t ; Xt )] = E X w( Xt )g(A t ; Xt ) |A t , w( Xt )g(A i ; Xi ), ∀i < t = g(A t ) + β t which holds since β t = 0 and by setting g(A t ) = ∇ A t T (A t ).

• A.5.2.3: ∇ A t T (A) is continuous, which can be immediately verified from the definition of the objective.

• A.5.2.5:

t α t |β t | < ∞, (15)
which holds as β t = 0 ∀t.

In addition, a condition on the constraint set H must be verified. As the same constraint is presented in [START_REF] Zheng | Stochastic resource optimization of random access for transmitters with correlated activation[END_REF] for the case without errors, we refer the reader to the statement and proof of the constraint condition in [4, Theorem 1].

Theorem 1 shows that even when there are errors, Alg. 3 converges to stationary points of the problems in (3). This is in contrast to the standard SGA algorithm in Alg. 2, for which these guarantees only hold when there are no errors.

C. Proposed Algorithm

While the weight in (13) yields convergence of Alg. 3 to a stationary solution, it must first be computed. In practice this can be challenging for two reasons:

(i) the true activity distribution may not be perfectly known or is difficult to sample from; (ii) the distribution of the estimated activity X must be computed.

We propose the following solution to address (i) and (ii). In order to perform user identification, an estimate of Pr(X = x) is required for the prior within the HGAMP algorithm. This prior is typically estimated via the expectation-maximization algorithm coupled with HGAMP [START_REF] Ke | Compressive sensingbased adaptive active user detection and channel estimation: Massive access meets massive MIMO[END_REF]. In the computation of the weight in ( 9), the probability Pr(X = x) can therefore be approximated by the prior used within HGAMP.

The HGAMP algorithm also produces estimates of the channels h t n . As such, samples of X t n and y t in (2) can be simulated using the prior estimate of Pr(X = x). Let ϵ(p) be the error probability of HGAMP for a given true activity p and e ∼ Ber(ϵ(p)), as a consequence, an estimate Pr( X = x) of Pr( X = x|e = 1) (the posterior of HGAMP given that an error is made) and ϵ(p) can be obtained via Monte Carlo estimation. This gives the following approximation for the weighting function:

ŵ1 (x) = Pr(X = x) (1 -ϵ(p))Pr(X = x) + ϵ(p)Pr( X = x) . (16)
Unfortunately, when N is very large, Monte Carlo estimation of Pr( X = x) and ϵ(p) has a high sample complexity. It is therefore desirable to consider further approximations. A convenient choice is

P (x) = (1 -ϵ(p))Pr(X = x) + ϵ(p) 2 N -1 . (17)
Here, (17) is obtained by considering every possible vector to be equally likely whenever an error is made by HGAMP. An approximation of the weight in ( 13) is then given by

ŵ2 (x) = Pr(X = x) P (x) . (18)
If the number of users is large, ϵ(p) can also be approximated to a value close to one (as it will be very likely to miss at least one user).

Accounting for the challenges in computing the weight in (i) and (ii) then leads to Alg. 4 where ŵ is either equal to ŵ1 or ŵ2 . Due to the approximations used in the computation of these weights, the bias in the estimated gradients is not zero, as the importance sampling property does not hold anymore, unlike the weight in [START_REF] Zheng | Stochastic resource allocation for outage minimization in random access with correlated activation[END_REF]. Nevertheless, as we show in the Sec. V, significant improvements in the performance are obtained over the standard SGA algorithm in Alg. 2, which does not account for the presence of errors in user identification.

V. NUMERICAL RESULTS

While we have established that Alg. 3 converges to a stationary point with probability one, this is under the assumption that the weight can be computed with no approximation error. In practice, Alg. 4 must be used with the weight approximated by ŵ1 or ŵ2 . This inevitably leads to bias in gradient estimates and hence a reduction in the final expected throughput. To investigate the reduction compared with the case of perfect knowledge, we perform a numerical study focusing on two scenarios.

In the first scenario, a network consists of N = 10 devices and an AP, where each active device transmits on one of K = 4 slots. Each device is assumed to transmit independently with probability p n , with p = .19 .20 .21 .23 .31 .55 .71 .82 .85 .94

In the preamble phase, user identification errors are assumed to occur with probability ϵ. In the case an error occurs, the estimated activity vector X is assumed to have independent components with Pr( Xn = 1) = 1-p n . This error distribution corresponds to a worst case scenario and is therefore of interest to validate that Alg. 4 can effectively mitigate the impact of the errors on slot allocations.The initial matrix A 1 is randomly selected by setting the elements to A i,j = 1 K + z ij , with z ij ∼ N (0, 0.05) and then projecting the resulting matrix into H.

Fig. 2 plots the expected throughput with optimized allocation matrix A for varying error probability ϵ. Five algorithms are considered: perfect knowledge (Alg. 2) and Aloha (uniform slot allocation for each user, as in III.1) which serve as references; the unweighted case corresponding to Alg. 2 but with a biased estimate of ∇ A T (A) due to the errors; and two curves corresponding to the proposed Alg. 4 with weights ŵ1 and ŵ2 , respectively. The step-size for each algorithm is set to 1 50 and the algorithms run for f = 5000 iterations. The results are obtained by averaging over 10 runs.

Observe that as expected, the expected throughput with perfect knowledge yields the best performance, while there is significant performance degradation in the unweighted case based on Alg. 2. We highlight that the proposed algorithm in Alg. 4 significantly outperforms the unweighted case for both weights ŵ1 and ŵ2 . For most values of ϵ, Alg 4 is also seen to yield performance close to the case of perfect knowledge, which suggests that the impact of user identification errors is largely mitigated. Note also that unless ϵ = 1, all the proposed methods outperform Frame Slotted Aloha.

In the second scenario, Alg. 4 is again compared with the unweighted case, the perfect decoding case and an Aloha system in a network with N = 10 devices, K = 2 slots and activity vector governed by p = .11 .17 .18 .23 .25 .34 .57 .63 .70 .94 . In contrast with the first scenario, the distribution of the estimated activity Pr( X = x) is based on user identification using the HGAMP algorithm which exploits circular complex Gaussian pilot sequences of length L = 5. Note that L < N meaning that active users are likely to mixed by the decoder. The weights ŵ1 and ŵ2 are estimated from 3000 prior simulations. The initial matrix is obtained as previously. Each algorithm is run 15 times for f = 15000 iterations with a constant step-size of α t = 1 20 . HGAMP algorithm often yields poor user identification and there are limited differences in the performance. On the other hand, as the SNR is increased, while there are still errors it can be observed that the proposed Alg. 4 with weight ŵ1 yields a significant performance improvement over the unweighted case. The reason for the poor performance of weight ŵ2 is that the assumption of uniformly distributed activity in the case of an error is a poor approximation, especially if there are many users. It is interesting to note even when the SNR is small, the algorithms are still outperforming Aloha.

VI. CONCLUSION

Stochastic resource allocation in grant-free random access relies on estimates of device activity via user identification algorithms, such as HGAMP. However, in practice errors are introduced which lead to potentially significant performance reductions. In this paper, we have developed an algorithm that can mitigate the impact of errors via importance sampling. Numerical results show significant performance improvements over existing algorithms, in many cases approaching the performance when user identification errors are not present.

Fig. 1 :

 1 Fig.1: Expected throughput with user identification errors for a network of N = 4 devices sharing K = 2 slots. Drawing the sample with probability ϵ from p instead of p (in orange) might introduce a 60% reduction of the throughput (if ϵ = 1) and lead to a throughput worse than an Aloha network. The blue line corresponds to a perfect identification (ϵ = 0).

Fig. 2 :

 2 Fig. 2: Comparison of the different algorithms for varying error probability ϵ and activity governed by (19).

Fig. 3 :

 3 Fig. 3: Expected throughput of the different algorithms as a function of the SNR.

Fig. 3

 3 Fig. 3 plots the performance of each algorithm for varying SNR in the preamble phase. For low values of SNR, the

Algorithm 4 :

 4 Stochastic Optimization Algorithm with Approximate Weight Function 1 Choose initial allocation matrix A 1 ∈ H, and step-size sequence {α t } with α t > 0, t = 1, 2, . . . Based on the estimate Xt , compute a biased estimate g(A t ; Xt ) of ∇ A t T (A t )

	2 t← 1.
	3 while not converged do
	4
	5	A t+1 = Π H [A t + α t ŵ( Xt )g(A t ; Xt )]
	6 7 end t ← t + 1


domain: phys spi

text:
Joannès Chambon

Jérôme Antoni

Simon Bouley

Galerkin equivalent sources

published or not. The documents may come

Galerkin Equivalent Sources Method for sound field reconstruction around diffracting bodies Joannès Chambon, 1, a) Jérôme Antoni, 1 and Simon Bouley. 2 1) Univ Lyon, INSA Lyon, LVA, 25 bis av. Jean Capelle F-69621, Villeurbanne Cedex, France.

2) MicrodB, 28 Chemin du Petit Bois, F-69134, Écully Cedex, France.

(Dated: 3 October 2022)

The rising interest for three-dimensional acoustic imaging requires the improvement of the numerical models describing the propagation between a radiating body and a microphone array. The commonly used free field transfer functions boil down to assume a full acoustic transparency of the radiating object, which in some case may lead to misleading outcomes for their characterization. Among other approaches, Equivalent Sources Methods (ESM) emerged as a convenient and powerful approach to simulate scattered sound fields. In the following paper, an acoustic imaging algorithm named Galerkin ESM where equivalent sources are tailored to concomitantly match with microphone pressures and a Neumann boundary condition is proposed. By means of a projected matrix inversion and a backpropagation of the equivalent sources, Galerkin ESM aims at the direct synthesis of the pressure field around a diffracting body by making the most of an array measurement. This method is compared with two other existing imaging algorithms fueled by both free field and computed transfer functions. The impact of the chosen transfer model is discussed, and Galerkin ESM performances are evaluated on both numerical and experimental test cases.

[https://doi.org(DOI number)]

[XYZ] Pages: 1-15

I. INTRODUCTION
Acoustic propagation models that account for rigid body diffraction can be obtained with a wide range of numerical methods. When it comes to problems of infinite extent, the Boundary Element Method (BEM) emerges as the main option for its strong theoretical foundation and its convenience (see for instance [START_REF] Bonnet | Boundary integral equations methods in solids and fluids[END_REF] or [START_REF] Burton | The application of integral equation methods to the numerical solution of some exterior boundary-value problems[END_REF] for an overview). As it got enhanced over the years to tackle its non-uniqueness issues, BEM became a complete but more complex and time consuming algorithm that requires consequent computational resources. Over and above the implementation of the method itself, BEM also requires a fine mesh management expertise. The use of advanced iterative solvers (see [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] or Fast Multipoles introduced by [START_REF] Rokhlin | Rapid solution of integral equations of scattering theory in two dimensions[END_REF] led to significant speed ups, but at the same time turned BEM into a rather hard to master option.

With a view to provide a more simple and flexible alternative for more simple test cases, [START_REF] Koopmann | A method for computing acoustic fields based on the principle of wave superposition[END_REF] designed and assessed an Equivalent Sources Method (ESM). Its physical principle can be stated concisely: considering a rigid body exposed to an incident acoustic field, a volume interior distribution of elementary sources can be set to offset the boundary condition induced by the incident field on the skin of the object. The main interest of having equivalent sources strictly in a) Also at MicrodB. Author to whom correspondence should be addressed. Electronic mail : joannes.chambon@microdb.fr the interior domain is to avoid numerical singularities issues encountered with BEM. These sources are thus considered as acoustically equivalent to the presence of the rigid body, and can be propagated toward the domain outside the skin to get the scattered sound field. After describing the theoretical background of ESM, [START_REF] Koopmann | A method for computing acoustic fields based on the principle of wave superposition[END_REF] also established its numerical version in which the rigid body and the volume of equivalent sources are respectively modeled by a rudimentary mesh and a discrete set of acoustic monopoles. Since then, ESM has been extensively studied and [START_REF] Lee | Review: The use of Equivalent Source Method in computational acoustics[END_REF] proposed a review article where the key parameters of its proper implementation are identified:

The optimal number equivalent sources was initially discussed by [START_REF] Koopmann | A numerical solution for the general radiation problem based on the combined methods of superposition and singular-value decomposition[END_REF] on a cylindrical test case. Later on, [START_REF] Dunn | Aeroacoustic scattering via the equivalent source method[END_REF] experimentally showed that a ratio of three times less sources than nodes on the mesh provided the best results on practical cases.

Their spatial distribution proved out to be directly linked to the conditioning of ESM matrices. Notably, a slight distinction was pointed out by [START_REF] Leblanc | A wave superposition method based on monopole sources with unique solution for all wave numbers[END_REF] between ESM with regularly distributed equivalent sources (Method of Fundamental Solutions, see [START_REF] Chen | Using the method of fundamental solutions in conjunction with the degenerate kernel in cylindrical acoustic problems[END_REF][START_REF] Kondapalli | Analysis of acoustic scattering in fluids and solids by the method of fundamental solutions[END_REF] and ESM using randomly located sources (Wave Superposition Method, see again [START_REF] Koopmann | A method for computing acoustic fields based on the principle of wave superposition[END_REF]. They showed that the first is a wellposed problem but suffers from non-uniqueness issues at particular frequencies while the latter is more robust but prone to ill-conditioning deficiencies. This can be dealt with using singular value regularization (see [START_REF] Lee | Assessment of timedomain equivalent source method for acoustic scattering[END_REF]. [START_REF] Pavić | An engineering technique for the computation of sound radiation by vibrating bodies using substitute sources[END_REF] also designed a time consuming but accurate algorithm to determine optimised spatial configurations.

Their retreat distance toward the skin of the rigid body is also a key driver for the proper functioning of ESM approaches. The stake is comparable to BEM collocation techniques : equivalent sources far from the mesh struggle to precisely match with complex incident acoustic fields while sources close to the skin are likely to be numerically unstable because of singularities. [START_REF] Bai | On optimal retreat distance for the equivalent source method-based nearfield acoustical holography[END_REF] recently tackled this issue and proposed leads to find the optimal balance. Lastly the nature of the equivalent sources plays a role in the efficiency of the method. ESM can be described in terms of single and double layer potentials (see [START_REF] Wilton | A clarification of nonexistence problems with the superposition method[END_REF], and [START_REF] Jeans | The wave superposition method as a robust technique for computing acoustic fields[END_REF] first investigated the use of dipolar equivalent sources, while [START_REF] Ochmann | The source simulation technique for acoustic radiation problems[END_REF] later on suggested that any type of sources (like spherical harmonics, see [START_REF] Bouchet | Calculation of acoustic radiation using equivalent-sphere methods[END_REF] could be used.

Featuring all these refinements (see [START_REF] Lee | Review: The use of Equivalent Source Method in computational acoustics[END_REF], for a more in-depth review), ESM is a relevant alternative to BEM for its straightforward implementation and its accuracy on a large scope of test cases.

Bearing in mind the advantages cited above, ESM naturally got introduced into acoustic imaging methods over the last years. The generic set-up common to all acoustic imaging problems,

p = Hq, (1)
where a source grid q has to be identified from microphonic pressures p and a collection of acoustic ) for imaging in enclosed space. If they differ in the algorithm used to identify the sources, the common ground to all these references lies in their use of a discrete set of equivalent monopoles (or dipoles, see for example [START_REF] Valdivia | Advanced equivalent source methodologies for near-field acoustic holography[END_REF] to reconstruct a radiated acoustic field from array measurements. However, another common aspect is that they no more involve these equivalent sources in the simulation of a particular boundary condition in contrast to the historical purpose of ESM.

A few tries to align the use of equivalent sources in acoustic imaging with scattering simulation are to be mentioned : Le Magueresse (2016, section 6.2.4) first introduced the idea of using Koopman's ESM to fill the H matrix in Eq.( 1) with coefficients taking into account the diffracting behaviours in the acoustic scene. Later on, Le Magueresse et al. (2020, section 4) also highlighted how the use of such transfer functions can matter to avoid misleading interpretations of source identification on an experimental engine test bench. At the same time, [START_REF] Chambon | 3D Beamforming for wind tunnel applications using ESM based transfer functions[END_REF] used the same approach in a wind tunnel, and showed how ESM could be used for aeroacoustic source identification on a car mirror to improve beamforming maps resolution.

Drawing on this, the connection between Koopman's ESM approach and acoustic imaging remained to be formalized, which is the point addressed in this article. The concept of fitting a set of equivalent sources to conform to a given boundary condition is borrowed from ESM, and included in the acoustic imaging inverse problem. Considering a Neumann boundary condition, this embedding takes the form of an orthogonal projection of the equivalent sources on the kernel of a matrix modeling the impedance of a rigid body. The basic inverse problem in Eq.(1) remains, but finds itself restricted to solutions satisfying an orthogonality constraint similar to those stemming from Galerkin methods. The aim of such an approach is also a crossover between classical ESM and acoustic imaging : the data provided by the phased microphone array is processed to identify virtual equivalent sources, and the latter are to be propagated to synthesize the acoustic fields at areas of interest.

The proposed algorithm, hence named Galerkin ESM, is depicted in this paper. Section II is dedicated to its mathematical formulation, from the statement of the inverse problem to the backpropagation step providing the display of the pressure maps synthesized from equivalent sources. Then, in section III, the emphasis is put on an in-depth study of the orthogonal projector involved in Galerkin ESM. A lead to precondition the kernel matrix with a physical maximization of the radiation efficiency of the sources is exposed. In section IV, Galerkin ESM is put to the test head-to-head against the well-assessed acoustic imaging methods that are Conventional Beamforming (CBF) and iterative Bayesian Focusing (iBF) on a various set of test cases featuring scattered sound propagations. With a full control on the target pressure fields, this benchmark is an opportunity to discuss the importance of supplying imaging algorithms with realistic transfer functions instead of the commonly used free field models. The last section consists in an illustrative application of Galerkin ESM on an academical experimental set of measurements.

Conventions and notations

Throughout this article, vectors are represented by bold lowercase letters and matrices by bold uppercase letters. M H denotes the transpose conjugate of a complex matrix, and M + stands for its Moore Penrose inverse. The orthogonal complement of a given Hilbert space H is noted as H ⊥ . The complementary of a given set E is Ē. The euclidean norm of a vector v is written ∥v∥ and its p-norm ∥v∥ p . The column space of a matrix M is noted as Im M and its kernel as ker M.

Lastly, the convention chosen for phase sign is e -iωt .

II. GALERKIN ESM A. Overall methodology

First is considered a three-dimensional radiating body bounded by a closed surface Γ. The acoustic field is sensed via a set of M microphones located at the positions (r i ) i≤M . The purpose of the method is to provide a description of the exterior acoustic field (both phase and amplitude) at a given pulsation ω by means of equivalent sources placed inside Γ. The latter are to be determined with respect to the measured pressures p(r i ) i≤M on the one hand, and to the impedance condition induced by the behaviour of the body on the other hand.

A discrete formulation of the problem is proposed now. Γ is modeled by a surfacic mesh featuring N nodes of positions (r j ) j≤N and a set of N s equivalent monopolar acoustic sources are introduced at the positions (r k ) k≤Ns in the volume Ω.

Let G ∈ C M ×Ns be the matrix defining the acoustic free-field transfer functions between the equivalent sources and the pressure at the microphone positions, i.e. ∀i ≤ M, ∀l ≤ N s ,

G il = -iωρ e ik∥ri-r l ∥ 4π ∥r i -r l ∥ , r i ∈ Ω, r l ∈ Ω. ( 2
)
Similarly can be defined the matrix G Γ ∈ C N ×Ns describing the source to pressure transfer between the equivalent sources and the nodes on Γ, and T Γ ∈ C N ×Ns the source to normal velocity transfer ∀j ≤ N, ∀l ≤ N s , T Γ,jl = e ik∥rj -r l ∥ 4π ∥r j -r l ∥

(1 -ik ∥r j -r l ∥) cos θ jl , r j ∈ Γ, r j ∈ Ω.

(

) 3
where θ jl denotes the angle between r j -r l and the direction normal to Γ at r j .

Finally the impedance imposed by the radiating body is introduced as

z j = p(rj ) vn(rj ) , r j ∈ Γ j≤N .
After stacking the M pressures p(r i ), the N s equivalent sources flow q(r l ), and the impedances z j into the column vectors p, c and the diagonal matrix Z respectively, the inverse problem under scrutiny boils down to

p = Gc (measurements matching) (ZT Γ -G Γ ) c = 0 (boundary condition). ( 4
)
From a mathematical standpoint, the goal here will be to solve a classical acoustic imaging inverse problem whose solution has to be projected on a null space accounting for the presence of the boundary condition.

Assuming Eq.( 4) solved, the last step of the approach consists in the repropagation of the obtained sources c on any observation region Γ obs thanks to free field transfer matrices defined similarly to Eq.( 2) and (3).

In this paper, the particular case of a perfectly rigid body is investigated for the sake of simplicity. In other words, the admittance on Γ is set to zero which implies

T Γ c = 0 (5)
to account for a Neumann boundary condition, and the overall problem as stated in Eq.( 4) becomes

p = Gc (measurements adequation) T Γ c = 0 (boundary condition). (6) B. T Γ Null space
As the acoustic inverse problem turns into an algebraic projected inversion, the settings of the operating matrices dimensions is a major stake for the practical implementation of the method.

M and N are considered as inputs of the problem: the number of microphones is given by the available hardware, and the number of nodes describing the body is fixed in accordance with the minimum wavelength under scrutiny. Concerning the test cases presented in this paper, a 5 nodes per wavelength condition is matched to ensure a margin [START_REF] Lee | Assessment of timedomain equivalent source method for acoustic scattering[END_REF] showed that ESM approaches could be more resilient than FEM/BEM on this point).

The number of equivalent sources N s is a key parameter of the method and requires further clarification. As it is mandatory that Eq.( 5) admits non zero solutions for the proper operation of the method, the lowest N s is first given by the rank theorem:

N s = Rank (T Γ ) + dim (ker T Γ ) . (7)
As T Γ is of full rank (i.e. Rank (T Γ ) = min (N, N s )), a non-empty kernel implies N s > N and the number of equivalent sources must be greater than the number of nodes supporting the boundary condition on the body.

Under this assumption, the resolution of Eq.( 4) involves the extraction of the kernel of T Γ . The simplest way to do so is based on the classical relationship

ker T Γ = ImT H Γ ⊥ . ( 8
)
From then on, the QR decomposition (see for example [START_REF] Golub | Matrix Computations[END_REF] of

T H Γ is sufficient to sample a K sized matrix B ∈ C Ns×K from the null space of T Γ . Let us first introduce T H Γ = QR =       Q 1,1 • • • Q 1,N Q 1,N +1 • • • Q 1,Ns Q 2,1 • • • Q 2,N Q 2,N +1 • • • Q 2,Ns . . . • • • . . . . . . • • • . . . ImT H Γ Q Ns,1 • • • Q Ns,N null space of T Γ Q Ns,N +1 • • • Q Ns,Ns       R, (9)
where Q ∈ C Ns×Ns is an unitary matrix and R ∈ C Ns×N is an upper triangular matrix.

Equation ( 9) can be rewritten as

T H Γ = Q 1 Q 2 R, (10)
with TΓ) describing the image of T H Γ and Q 2 ∈ C Ns×(Ns-Rank(TΓ)) being a set of orthogonal rows in the kernel of T Γ . From there, any matrix B of the form

Q 1 ∈ C Ns×Rank(
B = Q 2 Λ, Λ ∈ C (Ns-Rank(TΓ))×K (11)
ensures that T Γ B will equal zero. The B matrix is the cornerstone of the proposed Galerkin ESM as it plays two major roles:

Because of its similarity with Q 2 , it can force equivalent sources to respect the Neumann boundary condition stated in Eq.( 5).

The appropriate construction of Λ allows the reduction of the the inverse problem dimension in Eq.(4) from N s to K unknowns.

Leads to define what can be a relevant definition of this matrix is investigated in section III.

C. Practical inversion

The initial system (4) can now be reduced to a single inverse problem

p = Gc = GBd (12)
where the new unknown vector d ∈ C K may be understood as a set of equivalent source coefficients in a basis accounting for the Neumann boundary condition on the body.

Since it is desirable for this new statement of the problem to be correctly conditioned, a parametric condition arises from Eq.( 12) as it should not be underdetermined. This means that the size of the kernel basis K has to be larger than the output dimension M of the matrix GB to invert, which necessarily leads to

M ≤ K ≤ N s -N. ( 13
)
This prerequisite on the number of equivalent sources will be supported by numerical considerations in section IV.

Now that the inverse problem boils down to the formulation Eq.( 12) remains the choice of an inversion algorithm to reach the solution d. A large scope of sophisticated algorithms are available in the literature, (see for example [START_REF] Leclère | A unified formalism for acoustic imaging based on microphone array measurements[END_REF] or Merino-Martínez et al.

(2019) for a review). The shape of the optimal solution d (i.e. for example sparse or with minimum energy) is casedependant and likely to be influenced by the choice of Λ. For that reason, the generic choice would be to make use of an L 2 regularized inversion. This is the one made in the framework of this article in the validation sections, with the regularization parameter set up according to [START_REF] Pereira | Empirical Bayesian regularization of the inverse acoustic problem[END_REF] Bayesian regularization algorithm.

Once c is recovered from Eq.( 12), the last step consists in the backpropagation of the optimal equivalent sources vector ĉ = B d to estimate the radiated pressure in the region of interest within Ω. This is simply achieved with the free field propagator following

p obs = GΩĉ. ( 14
)
For this backpropagation to be meaningful, a noticeable precaution has to be highlighted: since the equivalent sources arising from Galerkin ESM are designed to acoustically model the object defined by Γ, it is mandatory for this region of interest not to intersect Γ. This would indeed boil down to evaluating the pressure radiated by the body inside or on itself, and most likely lead to incoherent results and interpretations. Indeed, like classical ESM or BEM, the Galerkin inverse problem is dedicated to the simulation of the acoustic field outside the source region since the 1/r singularity of the monopoles is unmanageable within Ω. Apart from this area, it is supposed to provide a realistic synthesized acoustic field outside the rigid body based on the microphone array measurements.

In the end, the overall execution of Galerkin ESM is summarized in Algorithm 1.

∩ Γ obs = ∅ G ∈ C M ×Ns ← Eq.(2). TΓ ∈ C N ×Ns ← Eq.(3). Q2, R ← QR decomposition of T H Γ . Λ ∈ C (Ns-Rank(T Γ ))×K ← User defined, see section III. B ← Q2Λ. d ← regularized inversion of p = GBd. p obs ← GΓ obs B d.
From the computational cost point of view, the two most significant steps are the QR decomposition (O N 3 operations using numpy Householder reflectors) and the inversion of GB that will depend on Λ and the chosen inverse methods.

III. KERNEL LAYOUT THROUGH PRINCIPAL SURFACES

When choosing the simplest form for Λ, i.e. a rectangular matrix with unitary entries on its diagonal, there is no genuine reason to expect any dimension reduction from N s -N to K: in that case B is directly sampled from the QR decomposition, which is a pure mathematical operation with no physical meaning involved. The selection of columns from Q 2 is thus bound to be arbitrary and it is shown in the next section that Galerkin ESM is unlikely to provide meaningful results with small values of K in this set up.

As mentioned in the previous section, the construction of Λ thus deserves an in-depth analysis since it impacts both the physical meaning of the equivalent sources model and the numerical cost of the inverse problem. The current section is dedicated to a practical way to define an acoustic operator from which are derived principal radiating surfaces. For the sake of clarity, only one promising option is advocated in this section, but some other leads could be studied as well depending on the application cases.

Sources with maximum radiation efficiency

The proposed option is to tidy up the columns of Q 2 following an order of increasing efficiency in terms of acoustic radiated power. For that purpose, the choice was to test Galerkin ESM featuring a Λ matrix filled with K eigenvectors associated to the largest eigenvalues of the radiation efficiency operator.

Given a closed surface completely covering the rigid body (Γ W ⊂ Ω) and discretized in N W nodes and surface elements, the acoustic power propagated by the equivalent sources c through Γ W is determined by

W = 1 2 ℜe v H n Ap = 1 2 ℜe (T Γ W c) H AG Γw c (15)
where T Γ W , G Γ W are the free field propagation matrices of c on Γ W respectively for normal velocity and pressure, and A ∈ R N W ×N W is a diagonal matrix containing the areas of Γ W surface elements.

From that perspective, it seems relevant to favour distributions of equivalent sources with the best radiation efficiency when reducing the size of the inverse problem. In practice this can be achieved by setting the columns of Λ as

Λ = [λ 1 , . . . , λ K ] ( 16
)
where

λ i = arg max ∥u∥=1 u H λj =0, j≤i-1 1 2 u H ℜe Q H 2 T H Γ W AG Γw Q 2 u. ( 17
)
All in all, the matrix B = Q 2 Λ brings in a combination of a projection on the kernel of T Γ while sorting the equivalent sources distribution in ascending order regarding their radiation efficiency on Γ W . Through this intermediary step, a physical behaviour is coupled to Q 2 and allows to significantly reduce the size of the inverse problem without damaging its accuracy. An illustration of how radiate such principal surfaces and the quantitative gain of the introduction of Λ defined with Eq.( 17) is put forth in section IV E.

Another subsidiary benefit lies in the orthogonality between the columns of B. Considering further applications of Galerkin ESM, one could think of decomposing the overall radiated acoustic power on the principal surfaces according to

W = K i=1 w i , where w i = di 2 Π i , (18)
with Π i being the i th eigenvector of the matrix 17) and di the i th coefficient of d. Such a process is likely to provide some insightful hierarchy in the operating radiative patterns (see Fig. 6,7,8 in the next section for a example on an academic case).

1 2 ℜe Q H 2 T H Γ W AG Γw Q 2 in Eq.(

IV. NUMERICAL VALIDATION

A comparison of Galerkin ESM with other popular acoustic imaging algorithms is proposed in this section. The purpose of the method is to reconstruct acoustic fields outside a rigid surface. To evaluate these methods from that angle, a spherical geometry is chosen: with respect to Fig. 1, Γ is modeled here by a spherical mesh of radius a = 0.3 meters featuring 1280 triangular elements and N = 642 nodes. As they get equally spread over the sphere, this set of nodes allows a 5 elements per wavelength condition up to ka = 15.

A. Spherical Related Transfer Functions

The acoustic field radiated around a rigid box (see for example [START_REF] Bouchet | Calculation of acoustic radiation using equivalent-sphere methods[END_REF][START_REF] Koopmann | A method for computing acoustic fields based on the principle of wave superposition[END_REF] is a common test case for ESM approaches. However, the point of this paper is also to put Galerkin ESM to the test in terms of three-dimensional diffracted source identification based on array imaging. From that perspective, the spherical case offers valuable advantages. First of all, the pressure in the vicinity of a rigid sphere given a monopolar source on its skin is analytically known (Williams (1999, section 8.8)): using the notations of Eq.( 2) and considering the sphere centered on the origin, the acoustic transfer function between a unitary monopolar source located at r l on the spherical surface Γ and a microphone at r i in the outer field is given by

H s il = -iρc 4πa 2 n∈N (2n + 1) h n (k ∥r i ∥) h ′ n (ka) P n r i • r l ∥r i ∥ ∥r l ∥ , ( 19
) where h n , h ′ n and P n respectively refers to spherical Bessel Functions of second kind and order n, its derivative, and Legendre polynomials of order n. [START_REF] Duda | Range dependence of the response of a spherical head model[END_REF] even proposed an accelerated routine to compute high orders of formulation ( 19), and Pereira (2013, chapter 2) empirically linked the maximum order n max to the frequency of interest through

n max > 1.2(ka) 2 + 8 r/a + 1 r/a (20)
to ensure a truncation error lesser than 10 -9 on the sum. Equation ( 19) provides the acoustic field resulting from any combination of monopoles on the surface of the sphere, and thus allows to assess Galerkin ESM on monopolar source identification. An observation surface Γ obs is set up for the display of the acoustic fields propagated by the equivalent sources from Galerkin ESM and the sources identified through classical acoustic imaging methods. It consists in an additional sphere of radius a obs = 4 3 a merged with a cutting plane of size 5a (see Fig. 2). This set-up allows to check discrepancies between the obtained propagated fields and the analytical one on both near and far field.

For the sake of simplicity, the choice was made to use a regular spherical array of radius 16 3 a with M = 250 microphones.

B. Acoustic imaging methods

For each validation case below, Galerkin ESM is benchmarked against two already existing acoustic imaging approaches.

The first one is Conventional Beamforming (CBF) (see [START_REF] Chiariotti | Acoustic beamforming for noise source localization -reviews, methodology and applications[END_REF]. Known to be fast and robust, CBF is theoretically justified for the identification of a single monopolar radiating source in a free field toward the array. Even if its relevance may be questioned given the more complex radiation patterns used here, it remains

with e j ∈ C N being the j th element of the canonical basis.

The second method is iterative Bayesian Focusing (iBF), a current state-of the-art algorithm initially proposed by [START_REF] Antoni | A Bayesian approach to sound source reconstruction : Optimal basis, regularization, and focusing[END_REF]. The full description of this algorithm as used in the scope of this article is to be found in Antoni et al. (2019, Algorithm 1), featuring a generalized multivariate complex Gaussian as prior density function. iBF was assessed as an imaging approach that provides the full cross spectral matrix of the sources, and includes a regularization step to deal with noisy measurements or poorly conditioned FRFs.

In the end, the comparison study features five different approaches to reconstruct the acoustic field:

1. CBF between N monopole sources on the mesh Γ and the array, featuring the free-field propagation of Eq.( 2) in the transfer matrix H of Eq.( 21).

2. CBF involving the analytical transfer Eq.( 19) in H.

With these first two is tested the usability of CBF for sound field synthesis and its ability to handle non free-field FRFs.

3. iBF between N monopole sources on Γ and the array, with free-field FRF and without any prior hypothesis on the sources in terms of regularization.

This one serves at evaluating what a turnkey version of iBF produces.

4. Galerkin ESM as described in section II, with N s = 3N = 1926 and K = N s -N = 1284 in a first stage, to maintain a margin regarding Eq.( 13). No principal surfaces for the validation in comparison with other imaging algorithms to remain general. The equivalent sources are distributed inside the mesh according to classical ESM literature guidelines : half of it are placed on a 85% scaled replica of Γ and the other half randomly located but at least 0.15a m away from the control points (see [START_REF] Leblanc | A wave superposition method based on monopole sources with unique solution for all wave numbers[END_REF].

5. Lastly iBF again, but with the use of the analytical FRF in H and a strong sparsity constraint (L 1 regularization). It should be noted that this last approach differs from the first four as it is the only one that features a suitable a priori on the sources. It aims at displaying the maximum degree of accuracy reachable on the reconstructed scattered field when the inverse method perfectly fits with the ground truth source distribution.

The aim of such a methodology is multiple:

Checking the error induced by approximated FRFs in CBF and iBF when the actual propagation is not free field.

Assessing the ability of CBF and iBF to benefit from the use of a perfectly accurate transfer function.

Evaluating Galerkin ESM comparatively to what is reachable with blind acoustic imaging approaches (i.e. without proper assumptions on the source distribution) and with the perfectly tuned iBF version.

C. Uncorrelated monopolar sources

The first test case consists in the reconstruction of the acoustic field produced by N src = 10 uncorrelated monopoles randomly placed on the yellow sphere (on Fig. 2). The cross spectral matrix (CSM) resulting from their combined radiation at the microphones is computed through

S pp = H s S qq (H s ) H , ( 22
)
where H s is determined according to Eq.( 19) and S qq is set to the identity matrix in order to describe completely uncorrelated sources of unitary strength. Each algorithm presented in sections II and IV B is then fed with this CSM and the sources obtained are finally backpropagated on a circle of radius 1.5a. The directivities resulting from this process are plotted on Fig. 3. This first configuration is supposed to fall in the range of all the methods assessed, but this first benchmark offers interesting preliminary conclusions. Methodwise, it seems that the overall directivity pattern is seized by both CBF, iBF and Galerkin ESM.

0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 •
It appears then that the use of the accurate FRFs leads to an enhancement at high frequencies, dealing with some 5 dB discrepancies for both beamforming and Bayesian Focusing. This item is quite intuitive and shows how the diffraction of the spherical body plays a role in the acoustic transfer and should not be underestimated.

With respect to the classical methods, Galerkin ESM (green dotted line Fig. 3) may be ranked at the second place in terms of accuracy. It outperforms CBF in general, deals with the impact of diffraction in a much better way than iBF with free-field FRFs, but does not compete with the latter when used with spherical transfer functions and induced sparsity.

It may be noticed that at both high and low frequencies, iBF with correct transfer functions (yellow dashed line) perfectly managed to reconstruct the reference directivity of the ten sources. This basically means that the number and the disposition of the microphones allows enough information to fully describe the acoustic field. That being said, it should be concluded that Galerkin ESM is by far the best performing algorithm among the without a priori approaches even if there still are room for improvements.

D. Correlated monopolar sources

The panel of imaging approaches is now put to the test for the recovery of an acoustic field produced by the same ten monopoles, but this time correlated: the CSM and the the reference directivity pattern is obtained in the same way as in the previous section using Eq.( 22), except that S qq is no longer the identity matrix. Instead, it gets defined as

S qq = I Nsrc + ∆ Nsrc + ∆ H Nsrc (23)
where ∆ Nsrc is an upper triangular matrix randomly filled with 0 or 1 entries. The ten sources on the sphere are thus either fully or not correlated one with each other, causing a much more complex and uneven directivity pattern than in the previous paragraph.

Results are plotted on Fig. 4. As expected, CBF loses relevance when dealing with correlated sources and fails at nailing the interference peaks. Regarding iBF, the scattering effect of the sphere becomes of first order at large ka number and the quality of the output is damaged as long as a free field transfer function is used. In that case, levels seem to be underestimated with respect to the reference: inherently to their propagation model, the identified sources are namely interferated through the spherical mesh at the backpropagation step of the process leading to overall lower levels because of the inappropriate L 2 regularization with respect to the source model.

Here again, Galerkin ESM does well as being in the efficiency range of the optimal set-up offered by iBF with spherical related transfer functions. Every interference lobe is well rendered at both high and low frequencies, only remain some local discrepancies with the reference. This uneven aspect may likely result from a non optimal disposition of the equivalent sources, which constitutes the major lead for further improvements of the method.

With a view to exemplify the potential applications of Galerkin ESM, the backpropagation depicted in Eq.( 14) is achieved and displayed on Fig. 5 for correlated monopoles (5 only for the sake of readability). Comparing (a) and (b), this result highlights to what extent the free assumption leads to misleading interpretations on the left side of the sphere. On this test case, it appears that iBF with accurate FRFs and Galerkin ESM only can provide a reliable backpropagated sound field.

All in all, the conclusion regarding the numerical test cases seems consistent. In terms of methods ranking, it can be stated that :

Beamforming is unusable at low frequencies for radiated field reconstruction. The poor spatial resolution on the identified sources leads to over smoothed propagated fields completely blurring directivity patterns. At large ka values, the use of accurate FRFs as steering vectors for CBF slightly improves it but this statement remains especially when dealing with correlated sources.

At high frequencies, i.e. when the scattering behaviour of the sphere is of first order, the identification of the outer pressure field with iBF should be achieved with great precautions. Without any prior knowledge on the sources nature, the generic parametrization of the algorithm with free-field transfer functions and L 2 regularization leads to significant inaccuracies. In practice, this means that iBF should be used for array-based field synthesis only when refined transfer functions and funded assumptions on the sources are available.

Galerkin stands out as a powerful alternative to have a proper simulation of the acoustic field around the rigid sphere without ground truth FRFs. Discrepancies with the reference are observable but local and acceptable. The peaks induced by the scattering presence of the sphere is well modeled by the kernel projection added to precondition the free-field transfer matrices as exposed in section II.

E. Principal Surfaces

Given its reasonable size, the spherical case used throughout this section is ideal to evaluate the contribution of principal surfaces introduced in section III.

A 5/4 scaled version of the mesh plotted in Fig. 2 was chosen as Γ W discretized in N W = N = 642 nodes. The Λ matrix of size N s × N s was filled with columns computed following Eq.( 17) and it is insightful to visually check their radiative pattern as basis vector for the Galerkin ESM inverse problem. For that purpose, the amplitude of the pressure field radiated by the columns of B on Γ obs , i.e.

0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 •
p λi obs = G Γ obs B i (24)
was computed for various i values and displayed in Fig. 6 with Λ being the identity matrix and Fig. 7 with Λ issued from Eq.( 17).

A relevant point arises from the observation of these figures and confirms the reasoning exposed in section IV E: On the one hand, in terms of acoustic radiation, no obvious hierarchy seems to sort the basis vector inside B when Λ is the identity. In other words, every columns of Q 2 share the same radiation efficiency and it is unlikely to sample K < N s -N of them without any loss on the accuracy on the Galerkin ESM outcome.

On the other hand, from Fig. 7 it appears that the power-optimized version of Λ leads to much more workable basis functions, with a clear ranking between the most radiating vector at low i values and the weakest ones at larger indices.

The quantitative assessment of this interpretation was conducted and shown in Fig. 8. The latter exposes the global L 2 relative error on Γ obs between the analytical pressure field and the one propagated by Galerkin ESM sources, with respect to the number of columns K selected from Λ. Simultaneously is plotted the relative acoustic power of the equivalent sources integrated on Γ W for rising truncation values K in the sum of Eq. ( 18).

It can be stated from this graph that the optimized formulation of Λ drastically improves the convergence speed of Galerkin ESM. On this particular test case for example, the -13 dB error ratio is obtained for K ≈ 0.8(N s -N ) = 1027 when B is only sampled from Q 2 while the same accuracy can be achieved with K ≈ 0.4(N s -N ) = 513 with the power oriented version of Λ. This earning represents a significant reduction of the dimension of the inverse problem in Eq.( 12) for practical applications of the method.

The red lines assesses the point made above with Eq. ( 18): when Eq.( 17) is used for the construction of Λ, most of the acoustic energy is concentrated in the first coefficients of d while at least 85% of them are required to reconstruct the same energy levels when Λ is simply the identity matrix.

F. Computational efficiency

Results above on the spherical test case were computed with Python 3.7 on an Intel Core I7 8750H (2.2 Ghz). Averaged on 25 runs of the code, 1.7 seconds per frequency were needed for the QR decomposition, 1.2 for the BΛ calculation in Eq.( 17) and 2.3 for the GB regularized inversion.

Given the results displayed in the current section, the fair comparison would be with iBF combined with realistic FRFs. Considering a generic test case, iBF itself is of the same complexity than the GB inversion but the FRF simulation has to achieved with FEM, BEM or classical ESM and it is likely that in the end the overall performance does not overrun Galerkin ESM.

V. EXPERIMENTAL VALIDATION

This last section lays out an experimental validation of Galerkin ESM. The aim of this final point is to gauge the behaviour of the method dealing with actual microphone array measurements on a large frequency range. The chosen case is a wooden mock-up flush mounted with omnidirectional sources, as described on Fig. 9. The latter is modeled with triangular mesh of N = 2377 vertices.

Its characteristic length equals L = 0.77 meters, and the planar array features M = 36 analogical microphones placed at 12 centimetres from the sources plan. Similarly to what was proposed in section IV, sources A and C were fully correlated while the last one remained uncorrelated. The measurement is a 9.84 seconds record sampled at f s = 25600 Hz.

The session took place in a semi anechoic room with a locally treated ground, and array signals were post processed using the Welch periodogram with a Hann window. The block size was set to 2560 samples with a 50% overlap to reach a 10 Hz frequency resolution up to 12800 Hz.

Since no analytical transfer function is available for such a geometry, the relevance of this last section namely lies in this transfer function computation step: whatever the choice of the numerical method (ESM, BEM, FEM, etc.), computing transfer functions involves numerical re- sources and may induce additional errors for complex meshes. The Galerkin ESM proposed in this paper could be seen as an alternative method that includes scattering effects in the model while bypassing this transfer simulation step.

With a view to recover the sound field propagated by the three sources around the engine, the same methods as in the previous section were applied here: CBF, iBF and Galerkin ESM, with the difference that this time only free field FRFs can be provided to the two first ones. The numerical modelisation of the set up is exposed on Fig. 10.

An illustrative display of the associated synthesized fields is given on Fig. 11, at a frequency at which scattering effects are likely to be significant (kL ratio far greater than one).

Preliminary observations of these acoustic maps are enlightening. As expected CBF is unsuitable to seize the acoustic produced by the correlated sources A and C. iBF and Galerkin ESM are in accordance on the side of the faces with the sources but differ on the part of Γ obs with pressure levels scattered around the mock-up. At first glance it seems that the acoustic shadow zones caused by the presence of the wooden structure find its most faithful depiction with the Galerkin ESM synthesized pressure field.

Another quantitative evaluation on the experimental case is displayed on Fig. 12. It exposes the backpropagated pressure level on reference microphone placed 12cm in the alignment behind the array shown on Fig. 9. For the sake of readability, the beamforming curve was post-processed by propagating only the three local maxima corresponding to the A, B and C points of the CBF map. At low ka numbers, Galerkin ESM all three methods yield a good adequation with the measured level, with a small edge for iBF and Galerkin ESM. However, beamforming rapidly tends to be globally inappropriate as frequency rises and iBF is subject to discrepancies at specific frequency ranges compared to Galerkin ESM. With regard to Fig. 11, this is consistent with the fact that acoustically masked areas in Ω do differ between iBF and Galerkin ESM issued backpropagation because of the misleading free-field propagator. Lastly, this last plot gives a first empirical clue on the maximum frequency reachable with the chosen number of equivalent sources concerning our method (here approximately) 1800 Hz). As for classical ESM, a sensitivity study will be conducted in further studies to define more precise guidelines on that specific issue.

VI. CONCLUSION

Galerkin ESM is introduced as a promising algorithm to reproduce sound fields scattered around diffracting bodies from array measurements. The transcription of a Neumann boundary condition as a matrix kernel inclusion is fully described, and integrated to the acoustic imaging inverse problem. A sub-process dedicated to the reordering of the kernel is also proposed to reduce the dimension of the resulting inverse problem, by maximizing the overall radiation efficiency the equivalent sources.

The method is assessed first on numerical data issued from correlated and uncorrelated monopoles flush mounted on a rigid sphere. Since highly realistic FRF can be analytically computed for such a geometry, this case first underlines to what extent the free field assumption becomes hazardous when it comes to reconstructing pressure fields at high ka ratios. From this first benchmarks, Galerkin ESM turns out to outperform CBF and iBF restricted to the free field model. For the sphere specific case, it however remains coarser than the inverse problem supplied with analytical ground truth FRFs.

Finally on experimental data, Galerkin ESM backpropagated maps provide insightful refinements on the sources directivity, especially compared to what actually used methods are likely to provide. In view of the initial motivation of the paper, i.e. scattered sound field synthesis from microphone array pressures, it appears that the approach pursued here takes on its full relevancy on this type of industrial application. Indeed, a complete set of realistic transfer functions is seldom available or at the expense of costly upstream computations, and it has been shown that Galerkin ESM provides an efficient framework to grasp strongly directive patterns while bypassing the FRF simulation step.

FIG. 1

 1 FIG. 1. (Color online)A set of equivalent monopolar sources with complex amplitudes c (blue dots) is sought after to match both the impedance condition z on the skin of the body and the pressure p measured by a microphone array at a given frequency.

Algorithm 1 :

 1 Step-by-step description of Galerkin ESM Require: Array measures p ∈ C M , N node positions on Γ, Ns equivalent sources positions in Ω, and a region of backpropagation Γ obs featuring N obs points. Ensure: M ≤ Ns -N and Γ

FIG. 2

 2 FIG. 2. (Color online) Geometrical set-up for validation. The yellow surface is the control point support Γ, the transparent grey is Γ obs and the cloud of points around it is the microphone array.

  FIG. 3. (Color online) Directivity pattern of 10 uncorrelated monopoles on a rigid sphere at ka = 1 (top) and ka = 13 (bottom), computed through various imaging algorithms. At low frequency, both Galerkin ESM and iBF are matching with the reference (all 4 plots are overlapped).

  FIG. 4. (Color online) Directivity pattern of 10 correlated monopoles on a rigid sphere at ka = 1 (top) and ka = 13 (bottom), computed through various imaging algorithms. At low frequency, both Galerkin ESM and iBF are matching with the reference (all 4 plots are overlapped).

  FIG. 5. (Color online) Acoustic pressure (amplitude) on Γ obs radiated by 5 randomly correlated monopoles at ka = 13 (with a 20 dB dynamic range). (a) Analytical sound field. (b) Sound field propagated by sources identified with iBF featuring free-field FRF. (c) Sound field propagated by sources identified with iBF featuring analytical FRF. (d) Sound field propagated by sources identified with Galerkin ESM.

  FIG. 6. (Color online) Pressure amplitudes on Γ obs radiated by the i th column of Q2Λ (Λ being the identity matrix). (a) i = 1, (b) i = 10, (c) i = 100, (d) i = 250 and ka = 8.24 . Plots are displayed within a 30 dB dynamic range, and the red dot indicates the maximum level location.

FIG. 10

 10 FIG. 10. (Color online) Numerical set-up for experimental set up. The yellow surface is the control point support Γ, the transparent grey is Γ obs and the 36 black dots represent the microphone array placed at 12 cm from the mock-up.

  FIG. 11. (Color online) Acoustic pressure (amplitude) backpropagated on Γ obs at kL = 21 (relative levels with a 30 dB dynamic range). (a) Sound field propagated by sources identified with CBF featuring free field FRF. (b) Sound field propagated by sources with iBF featuring free field FRF. (c) Sound field propagated by sources identified with Galerkin ESM.

J. Acoust. Soc. Am. / 3 October 2022 Chambon & al

J. Acoust. Soc. Am. / 3 October 2022Chambon & al

ACKNOWLEDGMENTS

The authors would like to thank Thibaut Le Magueresse (from Amiral Technologies, France and formerly from MicrodB) for making experimental results possible, and for his early stage intuitions on involving equivalent sources directly in the imaging inverse problem.