Direction-of-arrival estimation for acoustic signals based on direction-dependent parameter tuning of a bioinspired binaural coupling system

Bioinspired methods for sound source localization offer opportunities for resource reduction as well as concurrent performance improvement in contrast to conventional techniques. Usually, sound source localization requires a large number of microphones arranged in irregular geometries, and thus has high resource requirements in terms of space and data processing. Motivated by biology and using digital signal processing methods, an approach that adapts the coupled hearing system of the fly Ormia ochracea with a minimally distant two-microphone array is presented. Despite its physiology, the fly is able to overcome physical limitations in localizing low-frequency sound sources. By exploiting the filtering effect of the coupling system, the direction-of-arrival of the sound is determined with two microphones at an intermediate distance of 0.06 m. For conventional beamforming algorithms, these physical limitations would result in degraded localization performance. In this work, the bioinspired coupling system is analyzed and subsequently parameterized direction-sensitive for different directions of incidence of the sound. For the parameterization, an optimization method is presented which can be adopted for excitation with plane as well as spherical sound wave propagation. Finally, the method was assessed using simulated and measured data. For 90% of the simulated scenarios, the correct direction of incidence could be determined with an accuracy of less than 1∘ despite the use of a minimal distant two-microphone array. The experiments with measured data also resulted in a correct determination of the direction of incidence, which qualifies the bioinspired method for practical use in digital hardware systems.


Introduction
The localization of sound sources and the analysis of sound characteristics has applications in different scenarios e.g. in manufacturing industry, and health care. Sound source localization is used for quality control, predictive maintenance, and condition monitoring of mechanical systems/devices based on their emitted sound [1][2][3]. Sound analysis focuses on typical noise characteristics in the emitted sound of wear parts to optimize maintenance work/intervals [4,5]. An associated noise reduction can be part of preventive health measures [6,7]. Acoustic cameras are the respective technical tools to visualize the locations of sound sources [8,9] based on beamforming algorithms [10][11][12][13]. However, small distances between sensors in these cameras lead to small time differences of the arriving sound signals, which reduce the localization performance for low-frequency sounds. Due to the correspondingly large wavelength for low-frequency sounds, sufficiently large sensor distances are required which in principle increases the size of acoustic cameras [14]. Furthermore, an important factor in sound source localization is the number of microphones. The motivation is to minimize their number without losing localization performance. For this purpose, investigations with a minimum number of two microphones have already been conducted using state-of-the-art techniques such as interaural time difference (ITD) and interaural intensity difference (IID) [15][16][17][18].
The approach in this study incorporates insight from the biological example, the fly Ormia ochracea (O. ochracea). Ormia is able to determine the direction of incoming sound waves, although its interaural distance is comparatively small with respect to the detection task [19,20]. By studying the physiological characteristics of O. ochracea, Miles et al have shown that the mechanical structure of the fly's acoustic organs (tympanal organs) is the reason for its accurate localization performance despite the small sensor distance [19]. In contrast to human ears, there is no physical separation, but a mechanical coupling that connects the fly's acoustic organs and increases ITD and IID [21,22]. Due to this mechanical coupling, the parasitoid fly is able to accurately locate its cricket host, even though the tympanal distance of O. ochracea is 1.2 mm, while the wavelength of the cricket's call is about 0.07 m [22] at a peak frequency of 4.8 kHz [19]. In conventional sound source localization algorithms, this mismatch would result in reduced precision [23].
Miles et al studied the mechanical coupling between the tympanal pits inside the fly's head. It consists of an intertympanal bridge connecting two rigid bars at the pivot point through a coupling spring (k c ) and a damper (c c ). All moving elements of the system have the effective mass (m) which is assumed to be concentrated at each end of the intertympanal bridge [19]. The springs (k) and the dampers (c) at the extreme ends of the bridge represent the dynamic properties of the tympanal membranes, bulbae acustica and the surrounding structures [19]. The according system model is given in figure 1.
The system theoretical behaviors as well as the compatibility of miniaturization of the coupling system of Ormia were used in the optimization of antenna arrays [24][25][26][27], in the design of microphones using a built-in analog coupling system [28][29][30], in localizing underwater targets [22] and to estimate the direction-of-arrival (DOA) based on microphone arrays [23,31,32]. In [33], an approach to estimate the most sensitive frequency and the most sensitive angle of incidence of an incoming plane sound wave were considered. For the coupling parameters determined by Miles et al and Ormia's tympanal distance, this resulted in a frequency of 8341.9 Hz and an angle of incidence of 40.48 • . In particular, the ability to localize low-frequency sound sources with a minimally distant two-microphone array characterizes In addition, each mass is connected to a fixed reference surface by another spring k and damper c each. The parameters k, kc, c and cc are spring and damper constants. F1(t) and F2(t) represent the forces of the acoustic pressure field that is applied to the moving masses within the tympanal pits. l1(t) and l2(t) are the translational vibration displacements of the masses when the coupling system is excited [20]. the fly O. ochracea. Therefore, one goal of this work is to parameterize the bioinspired coupling system also for other angles of incidence, other excitation frequencies, as well as other microphone distances, and transfer this extended concept to a digital system.
In this work, first, an analytical description of the coupling system is given for the determination of the most sensitive angle of incidence. Second, a variational approach for a setup with two microphones is proposed that allows to fit the coupling parameters of the system model for different receiving directions at given excitation frequency and sound source distance. Third, an optimization approach for the variational parameters is developed that satisfies constraints to ensure comparable sensitivity for different parameter sets. Fourth, in addition to the usual simplification for plane waves, the whole approach is further developed for the more general case of spherical wave propagation. Finally, the functionality and accuracy of the new localization method are proven based on computer simulation and real experiments.

Bioinspired system model
A translational circuit model as shown in figure 1 is derived from the rotational model commonly used in literature as a basis for an analytical description of the principle of coupling two sound signals. The two sound pressures can be regarded as two sinusoidally varying point forces F 1 (t) and F 2 (t), each acting on a separate mass, representing the ends of the intertympanal bridge [20]. The system responses with displacements l 1 (t) and l 2 (t) of the two masses [19]. The forces are the inputs x i (t) and the displacements are the outputs y i (t) (with i ∈ {1, 2}) in a general description of the system model of the coupling system. According to the physiological structure of O. ochracea, this system model is a two-input, twooutput filter system [32,34].

Lumped mechanical system model of coupled tympanal organs
The lumped mechanical system model is based on a system of second order differential equations which can be derived from figure 1 (see also [19]): Here the inputs are written as vector x(t) and the outputs as vector y(t). Solving equation (1) using Laplace transform [34] leads to and as transfer functions for the two-input, two-output system with For this description approach, Y 1 (s) and Y 2 (s) each consists of two transfer functions, where P(s) is a fourth-order function. In [35], partial fraction decomposition is used to minimize the structure of the system. This leads to a linear combination of two second-order transfer functions with coupling parameters only in the denominator: Closer inspection of equations (7) and (8) shows, that the expressions only differ in the sign before the last term. In consideration of this, a signal processing model can be built with a cross-connection at the input and a second one before the output of the system, as shown in figure 2. Substituting s with jω in equations (7) and (8), results in two transfer functions depending on ω, where ω = 2π · f is the temporal angular frequency with f as the frequency of the propagating sound wave [10]: Block diagram of the coupling system. Input signals X1(s) and X2(s) are cross-connected, followed by the transfer functions H1(s) and H2(s). A second cross-connection is established in the signal path before the outputs Y1(s) and Y2(s).

Signal definition for the excitation of two microphones with plane waves
To define the input of the coupling system, a description of the pressure signals measured at the microphones is required. Those are defined in the time domain and consider the amplitude and the phase delay of a plane wave. In this case, it is assumed that the wave propagates from a source at infinite distance and that the amplitudes at both microphones are equal [10]. For the definition of the pressure signals p 1 (t) and p 2 (t) measured by the microphones m 1 and m 2 , a time delay τ must first be calculated. This is used to represent the different arrival times of the sound wave at the microphones. The time delay depends on the angle of incidence Θ in measured with respect to the array normal and the distance d between the microphones, as shown in the figure 3. Considering the propagation speed of the sound v = 343 m s −1 , the time delay is determined by Here the center point is in the middle between the microphones, so the sound wave reaches m 1 by τ 2 earlier than the center point. Microphone m 2 , on the other hand, is reached by the sound wave τ 2 later relative to the center point. In the time domain, p 1 (t) and p 2 (t) are shifted by half of the time delay calculated with equation (11) and have different signs: Fourier transform is applied to equations (12) and (13) to obtain the frequency domain input signals: with P p,i (ω) being the Fourier transform of p i (t). The signals have the same spectrum and differ only in their phase shifts X p,1 (ω, τ (Θ in )) and X p,2 (ω, τ (Θ in )).
Together with equations (9) and (10), the frequency response can be calculated to analyze the system behavior: showing that the coupling system here is less sensitive [33]. Moreover, at a frequency of 8000 Hz, a nearly 180 • phase rotation across the angle of incidence is most pronounced.

Frequency response of coupling system for different angles of incidence
Considering the standard parameters determined by Miles et al [19,20] as an example, given as k = 0.576 N m −1 , k c = 5.18 N m −1 , c = 1.15 × 10 −5 Ns m −1 , c c = 2.88 × 10 −5 Ns m −1 , and m = 2.88 × 10 −10 kg, the frequency response can be calculated when the coupling system is excited by different X p,1 (ω, τ (Θ in )) and X p,2 (ω, τ (Θ in )). Therefore, Θ in is varied in the range from 0 • to 90 • , while the temporal frequency remains unchanged. In a next step, X p,1 (ω, τ (Θ in )) and X p,2 (ω, τ (Θ in )) can be substituted into equations (9) and 10. In [33], the amplitude difference in [dB] and the phase difference in [ • ] between Y p,1 (ω, τ (Θ in )) and Y p,2 (ω, τ (Θ in )) are calculated over the entire interval of Θ in : Figure 4 shows the frequency response for plane wave excitation with frequencies of 5000 Hz, 6000 Hz, 7000 Hz, and 8000 Hz. It can be seen that the amplitude difference has a maximum that represents the most sensitive angle of incidence at a certain frequency and sensor spacing. At a frequency of 8000 Hz, this maximum is much more pronounced than at the other frequencies and a hill like characteristic is observed. Furthermore, it can be seen that the phase rotation is larger as the hill becomes narrower.
Hou et al [33,35] have shown an approach to determine the most sensitive angle of incidence. Considering this, it is possible to change the characteristic and especially the most sensitive angle of incidence using the coupling parameters.

Determination of the most sensitive angle of incidence to plane wave excitation
To determine the most sensitive angle of incidence to plane wave excitation, X p,1 (ω, τ (Θ in )) and X p,2 (ω, τ (Θ in )) must be inserted in the equations (16) and (17). The obtained transfer functions depend on ω and τ . Furthermore, Euler's formula is used to replace the exponential function by trigonometric functions. Thus, an analytical description for the upper and lower signal path (see figure 2) is obtained: Note that Y p,1 (ω, τ (Θ in )) and Y p,2 (ω, τ (Θ in )) both consist of two components, denoted R p,1 (ω, τ (Θ in )) and R p,2 (ω, τ (Θ in )), multiplied by the spectrum P(ω) of the input signal. Thus, Y p,1 (ω, τ (Θ in )) and Y p,2 (ω, τ (Θ in )) can be expressed as linear combination of R p,1 (ω, τ (Θ in )) and R p,2 (ω, τ (Θ in )) which differ by the sign of R p,2 (ω, τ (Θ in )). Each is corresponding to one of the natural modes of the coupling system [19]. R p,1 (ω, τ (Θ in )) belongs to a mode indicating a rocking motion and means that the ends of the tympanal bridge move in anti-phase. In contrast, R p,2 (ω, τ (Θ in )) considers a mode with translational motion characterized by in-phase movement of the two ends of the bridge [20]. The two terms can each be split into a real and an imaginary part, so that the imaginary unit disappears from the denominator: The use of equations (22) and (23) facilitates the maximization of the amplitude difference according to equation (18), neglecting the logarithm and using a limit value consideration. For this, it is assumed that the angle of incidence Θ in approaches the most sensitive angle of incidence Θ msa : when the denominator approaches zero (29) One solution to this is that R p,1 (ω, τ (Θ msa )) and R p,2 (ω, τ (Θ msa )) are equal, which implies that the absolute values and phases of both would also be equal. For phase equality this means: Figure 5. For illustration, the frequency response of the coupling system over the entire interval of Θ in is used as already known from figure 4 with an input frequency of 8000 Hz. In addition to the amplitude difference and the phase difference, the amplitude responses of Rp,1(ω, τ (Θ in )) and Rp,2(ω, τ (Θ in )) are depicted. At the maximum amplitude difference |Rp,1(ω, τ (Θ in = Θmsa))| and |Rp,2(ω, τ (Θ in = Θmsa))| are equal and the phase difference reaches the largest rate of change. This point represents the most sensitive angle of incidence.
Note that after the last transformation equation (33) no longer depends on the angle of incidence Θ in or on the most sensitive angle of incidence Θ msa , but only on the angular frequency ω (i.e. on the frequency f, respectively). However, since an expression for the most sensitive angle of incidence Θ msa is wanted, the requirement for phase equality can be discarded. Therefore, the equality of absolute values is considered instead: As can be seen in figure 5 for |R p,1 (ω, Θ in )| = |R p,2 (ω, Θ in )|, the amplitude difference is maximal at this point as well as the rate of change of the phase difference between Y p,1 (ω, Θ in ) and Y p,2 (ω, Θ in ). This point represents the most sensitive angle of incidence. Removing the square root and joining the terms results in This is simplified and written with the tangent function as Considering equation (11) and inserting it in the argument of the tangent function ω τ 2 , a term is obtained that contains Θ in and f as source-specific parameters: In addition, d depends on the distances of the microphones in the array and v on the propagation speed of the surrounding medium. Inserting the wavelength λ = v f in equation (38) results in: To obtain the most sensitive angle of incidence Θ msa for a parameter combination of the coupling system, equation (39) with Θ in = Θ msa is inserted in equation (37): This is finally solved for Θ msa : This expression is used as the basis for calculating the most sensitive angle of incidence for plane wave excitation while using a specific combination of coupling parameters.

Determination of the most sensitive angle of incidence for spherical wave excitation
The simplification of a plane wavefront is based on the assumption that the source is located far away from the sensor (far field). However, the more general assumption is spherical wave propagation in the near field. In addition to changing phase delays, the amplitude also decreases inversely proportional to the radial distance [36]. This more complex description of wave propagation is based on focusing on a point source at finite distance and assuming spherically propagating sound waves [10]. The corresponding setup is shown in figure 6. The finite distance is given as the radius r 0 of a semicircle whose center is in the middle of the microphone array. By targeting a focal point x f on the semicircle and calculating its radial distance r i to microphone m i , the measured signal of a spherical wave in the time domain can be defined as The corresponding Fourier transform is Figure 6. Linear array with two microphones m1 and m2 placed at a distance d. A spherical wave excitation is assumed. The sound source is located on a semicircle with radius r0. The center of the semicircle is located in the middle between the microphones. The angle of incidence Θ in is measured from the normal of the microphone plane at its center point. The radial distances from the sound source to the microphones are r1 and r2.
with the signal spectrum P(ω) and X s,i (ω, r i ) as steering vector of the spherically propagating waves [10]. Setting r i v = τ i , the inputs of the coupling system are obtained as Considering now the transfer functions from equations (9) and (10), X s,1 (ω, r 1 ) and X s,2 (ω, r 2 ) can again be inserted: Y s,2 (ω, r 1 , r 2 ) = 0.5 · P(ω) · (X s,1 (ω, r 1 ) + X s,2 (ω, As mentioned in section 2.1, P(ω) can be neglected for further calculations, since only the components containing X s,i (ω, r i ) are used. The same transformation steps as in section 2.2 are performed so that R s,1 (ω, r 1 , r 2 ) and R s,2 (ω, r 1 , r 2 ) are obtained. Using Euler's formula, this can again be written with trigonometric functions including the imaginary unit: Considering equations (24)- (27), these are substituted into R s,1 (ω, r 1 , r 2 ) and R s,2 (ω, r 1 , r 2 ) so that numerators and denominators of all terms can be shortened. In addition, R s,1 (ω, r 1 , r 2 ) and R s,2 (ω, r 1 , r 2 ) are simplified and split into real and imaginary parts because the absolute value is calculated in the next step and the imaginary unit must be removed from the denominator: ( 1 r1 cos(ωτ 1 ) − 1 r2 cos(ωτ 2 ))C C 2 + D 2 As already shown in equation (34), the absolute values of R s,1 (ω, r 1 , r 2 ) and R s,2 (ω, r 1 , r 2 ) are equated. The equation is rearranged such that all coupling parameters are on one side: ⇐⇒ r 2 1 − 2r 1 r 2 cos(ω(τ 1 − τ 2 )) + r 2 2 r 2 1 + 2r 1 r 2 cos(ω(τ 1 − τ 2 )) + r 2 For later use, this is written using two functions as: with the vector of coupling parameters ⃗ φ = [k, k c , c, c c , m] ⊤ . If the parameters in ⃗ φ satisfy equation (55) for selected ω, r 1 and r 2 , then the coupling system is sensitive to a matching excitation. Equation (54) in this form still depends on the radial distances r 1 and r 2 of the microphones from the sound source (see figure 6). In the next step, r 1 and r 2 are replaced by the angle of incidence Θ in , the radius r 0 and the fixed microphone distance d. This makes the system description in the near field comparable to that in the far field. On closer inspection of figure 6, the law of cosines from Euclidean geometry can be used to determine the radial distances r 1 and r 2 . Therefore two triangles are formed -i.e. one consisting of r 1 , r 0 and d/2 and the second of r 2 , r 0 and d/2: The expressions for r 1 and r 2 can be substituted in equation (41) respectively in (55) which is then written as: where the left side can be regarded as a target function which has to be satisfied by the right side which represents the actual coupling system.

Bioinspired sound source localization approach based on direction optimized system parameters
For the purpose of estimating the DOA, the coupling system should be sensitively parameterized for each possible direction of incidence, similar to the delayand-sum beamformer, which individually delays the measured pressure signals [10]. Therefore, the derived expressions for the most sensitive angle of incidence were used to optimize parameter sets for different target angles Θ trg,i representing the respective direction of incidence. The target angle was varied in the range of 0.25 • to 90.0 • in steps of 0.25 • : where i ∈ {1, 2, . . . , 360} was the target angle index for iteration during the optimization process. For this purpose, a minimum search was performed for a nonlinear constraint function, which yielded the optimal parameter vector function Since the minimum search for target angles smaller than Θ trg,1 = 0.25 • did not yield an optimal parameter vector, this was chosen as the minimum during optimization. In the following, the procedure is described for both plane and spherical wave excitation, using the expressions for the most sensitive angle of incidence as derived in sections 2.2 and 2.3.

Optimization for plane wave excitation
For plane wave excitation, equation (41) was used to calculate the most sensitive angle of incidence Θ msa in the minimum search for the difference between Θ msa and the respective target angle of incidence Θ trg,i . The minimum search returned the argument, in this case the optimal parameter vector function ⃗ φ opt (Θ trg,i ): subject to lower and upper boundaries and an additional expression ceq as nonlinear equality constraint: Considering figure 4, the equality argument in equation (63) was necessary to ensure that the system responded with the same amplitude difference for each parameter set when excited from the direction of incidence matching the parameter set. For this, the equations (20) and (21) were used to calculate the amplitude difference (which is represented by the fraction in equation (63) after the logarithm has been applied). Interaural time delay τ (Θ trg,i ) was calculated as described in section 2.2. The active-set algorithm [37,38] in Matlab R2021b (The Mathworks Inc. Natick, MA, USA) was used as optimizer with a step tolerance of 1 × 10 −12 . A maximum number of 8000 function evaluations was chosen arbitrarily after conducting pre-tests. For a larger number of function evaluations the error reduction per function evaluation was below 1 × 10 −6 . As mentioned, the nonlinear equality constraint in equation (63) was used to drive the optimizer towards the same maximum value (in [db]) in each run to keep optimization results comparable. For this, a threshold I thresh = 30 dB was set to ensure a narrow maximum for the amplitude difference (cmp. figure 5). If a larger I thresh had been chosen instead, the amplitude difference, which here is equivalent to the localization performance for this direction of incidence, would have been higher. In this case, however, the condition for Θ trg,i at the beginning (i.e. smaller than 3 • ) and at the end (i.e. greater than 80 • ) of the interval (see equation (59)) could not have been satisfied. This also means that comparability between the possible directions of incidence could not have been ensured. As a consequence, the parameters could have been optimized for the corresponding angle of incidence, but a different set of parameters would have resulted in a higher amplitude difference, even if it would not have been optimal for the exact angle of incidence. Through trial and error, the value of I thresh = 30 dB emerged as a good compromise between the inherently contradictory demands of a value that could be achieved over the entire working range of Θ trg,i while at the same time producing a narrow maximum for the amplitude differences.
In order to determine the lower and upper limits of the parameters to be optimized, the physical and mechanical boundary conditions had to be considered, respectively. For this reason, no negative values are allowed. Even if the optimizer seems to have found a minimum, negative values are unsuitable for the following sound source localization. A suitable procedure for defining the boundaries and start parameters will be described below.

Optimization for spherical wave excitation
In the case of spherical wave excitation, equation (58) was considered to minimize the difference between Ψ ω, Θ trg,i , r 0 and χ(ω, ⃗ φ(Θ trg,i )). Substituting it into the optimization specification, the optimized arguments ⃗ φ opt (Θ trg,i , r 0 ) for all Θ trg,i and r 0 were obtained: In addition, for constraint ceq in equation (66), the equations (46) and (47) were used taking into account the excitation by spherical waves as described in section 2.3. Otherwise, the same procedure as for the optimization for plane wave excitation was applied.

Exemplary optimization for fixed angular frequency and microphone distance
The goal in this exemplary optimization was to optimize the parameters of the coupling system χ(ω, ⃗ φ) such that equation (58) was satisfied. This means that χ(ω, ⃗ φ) -i.e. the right side of equation (58) -followed the course of the target function Ψ(ω, Θ trg,i , r 0 ) -i.e. the left side of equation (58). In figure 7, the target function Ψ(ω, Θ trg,i , r 0 ) was plotted as an example with ω = 2π · 1 kHz and a microphone spacing of d = 0.06 m for varying r 0 . The curves exhibit a sigmoid shape. For small r 0 the separate curves are distinguishable from each other, while for growing r 0 they overlap and converge against the target function for plane wave propagation as given in equation (40). From the figure 7, it can be seen that the curve of Ψ(ω, Θ trg,i , r 0 = 0.2 m) is noticeably different from the others, having a larger maximum at Θ trg = 90 • and a larger slope. It is important to mention that the nonlinear equality constraint must still be taken into account, which is why optimization must be performed separately for each r 0 . When choosing I thresh = 30 dB, it is possible to optimize the parameters with a high accuracy. However, this requires an increased optimization effort. Therefore, this was done only for r 0 = 0.2 m. For all other r 0 , an error for the nonlinear equality constraint of ±0.25 dB was allowed during optimization. This facilitated optimization and kept the loss of localization power low, i.e. the deviations from a peak amplitude difference of 30 dB were so small that comparability between the possible incident directions was still ensured. Moreover, the optimized parameters for r 0 = 0.2 m could be used as initial parameters for the optimization with r 0 > 0.2 m. The optimization was performed iteratively over all Θ trg,i with fixed ω and d and desired r 0 , starting from a parameter vector ⃗ φ 0 (Θ trg,i ). The aim was to obtain an almost continuous curve for the optimized parameters of the coupling system χ(ω, ⃗ φ), which approximately corresponds to the shape of the target functions Ψ(ω, Θ trg,i , r 0 ) shown in figure 7. This can be achieved if ⃗ φ 0 (Θ trg,i ) also follows the sigmoidal shape, which leads to the fact that also the optimized parameters adopt to the sigmoidal shape. For this purpose, Ψ(ω, Θ trg,i , r 0 ) was normalized to its maximum value assumed for Θ trg,360 = 90 • and multiplied by an initial parameter vector ⃗ φ init : During optimization, it was found that the parameters for k, k c , c c , and m do not become minimal at small Θ trg,i . For this reason, an offset was applied to the associated ⃗ φ 0 (Θ trg,i ) over the entire angular range. In the case of k, k c , and c c , the offset is half of ⃗ φ init . For m, the offset is 1 × 10 −6 kg, since this parameter was already in a very small range of values anyway. Finally, it had been shown that the parameter c does not converge to a curve with a sigmoid shape. In an unlimited optimization run, the parameter c followed some kind of sinusoidal function. It could be fitted by the sum of two sine functions according to the following equation that is obtained by using the Curve Fitting Toolbox in Matlab R2021b: As a result, for a radial distance of r 0 = 0.2 m, a well optimizable starting vector was now given. Using the boundaries as given in table 1 for optimization, it was possible to optimize a nearly continuous parameter curve that provided the desired incidence angles and peak amplitude differences. The values in table 1 were multiplied by ⃗ φ 0 (Θ trg,i ) to obtain lower (lb) and upper (ub) bounds for each Θ trg,i . The range of target angles was divided into two intervals, because at small target angles the parameters c and m become very small relative to the other parameters, so that especially the lower bound must be chosen smaller. The first interval, from 0.25 • to 9.75 • (i.e. for all i < 40) and second, for all larger target angles (i.e. for all i ⩾ 40). The optimized parameter curve is shown in figure 8.

Optimization for other radial distances
As described, the optimization of the start vectors was time-consuming. Therefore, to save time, the optimization process was organized in two steps, which included an initial optimization run and additional runs to minimize the optimization error. In the initial run, the results of the optimization with r 0 > 0.2 m were used as starting vectors (see also figure 7 for reference). The corresponding parameter vector could then be used as the start vector for the optimization, but for a different Θ trg . For all different target angles the changing nonlinear constraint had to be considered in the optimization. As before, the boundaries could be calculated as multiples of ⃗ φ 0 (Θ trg,i ) and can be found in table 2. In the improving runs, the optimization was repeated, using the optimized parameters of the previous step as the starting vector. Repeating the optimization only serves to further minimize the optimization error. The parameter curves for r 0 = 0.5 m, r 0 = 1.0 m and r 0 = 2.0 m are shown in figure 8.

Results
For a test, two simulated microphone signals x 1 (t) and x 2 (t) were used as inputs for the model of the coupling system, each signal in chunks of 4096 samples. For less than 4096 samples, the localization performance was lower because the modeled coupling system was not in a steady state. The output of the coupling system consisted of two time signals y 1 (t) and y 2 (t) with modified amplitude and phase shift, again handled in chunks of 4096 samples. To obtain a scalar value, the root mean square (RMS) of each signal was calculated for each chunk. The RMS of y 1 (t) was divided by the RMS of y 2 (t) and the result was expressed in decibels. This procedure was performed for all Θ trg,i to be tested with the set of parameters optimized as described in the sections 3.2.1 and 3.2.2 and which are shown in figure 8. The simulation model was implemented in Simulink (Matlab R2021b). It was conducted with a sampling frequency of 51.2 kHz using the ode14x fixed-step solver. Here, only positive angles of incidence are considered because the distribution of negative angles is identical to this. In (c), the peak amplitude difference can be seen for each simulated angle of incidence (positive only). In (d), the calculated 3 dB beam widths are plotted over all angles of incidence (positive only).

Simulation results
The simulation was conducted as described above, using simulated microphone data without noise from varying angles of incidence. For the signals, a frequency of f = 1 kHz was chosen and spherical wave propagation was assumed. The microphones were placed at a distance of d = 0.06 m. The angle of incidence Θ in was varied in a range of 1 • to 90 • in steps of 1 • . As already described in section 3, the radial distance r 0 of the sound source was varied over r 0 = {0.2 m, 0.5 m, 1.0 m, 2.0 m}. For all radial distances r 0 chosen, r0 d < 2.7, so that the error in the far-field approximation is less than 0.5 • [39,40]. Figure 9(a) shows the localization results for an angle of incidence of based on the calculated amplitude difference. For all simulated radial distances the selected angles of incidence were determined, i.e. that the direction of incidence of the sound source could be correctly localized. A positive amplitude difference represents a positive angle of incidence, while a negative amplitude difference indicates a negative angle of incidence. In figure 9(b) a histogram of the absolute error (AE) between the angle of incidence as simulated and the target angle with the highest amplitude difference is shown. For 90% of the simulations performed, the absolute error was between 0 • and 1 • . The resolution of tested ⃗ φ opt (Θ trg,i ) was the same as the step size for parameter optimization in quarters of a degree increments.
In figure 9(c), the compliance with the nonlinear equality constraint (see equations (63) and (66)) is quantified. The peak amplitude differences over all r 0 and all simulated angles of incidence Θ in are shown. It can be seen, for r 0 = 0.2 m the peak amplitude difference during optimization was 30 dB. In comparison, for the other r 0 an error of ±0.25 dB was allowed, which is also visible in the simulation results. In addition, the 3 dB beam width was calculated for all Θ in and plotted for varying r 0 in figure 9(d). As in figure 9(a), it can be seen that the 3 dB beam width increases with a growing angle of incidence. Simulations were also performed with angles of incidence Θ in > +90 • or Θ in < −90 • , i.e. the source was located in the area behind the microphone array. Due to the symmetrical design of the coupling system and the chosen propagation characteristic of the sound waves, which depends on the distance between the source and the microphone (see section 2), the coupling system behaves the same for a source in the backfield as its counterpart in front of the microphone array.

Experimental results
Based on the optimization results, the coupling system was tested with real microphone data. The experiments were performed in a soundproofed room protected from external noise and with a reduced reverberation time, which according to the simulation is documented as T 60 = 0.177 s at a frequency of 1000 Hz. Measurements were carried out with two free-field array microphones (GRAS 40PH-10 CCP, GRAS Sound & Vibration, Holte, Danmark) which were fixed in a 3D-printed mount at an inter-microphone distance of 0.06 m. The analog microphones were connected to a sound and vibration measurement module (NI-9234 in a cDAQ-9185 chassis, National Instruments Corp., Austin, Texas). This setup allowed for a sampling frequency of 51. During the experiment, the point source was placed at a radial distance of 0.5 m from the microphone array origin. The angle of incidence to the array was chosen to be 25 • . The experiment ran for 1.5 s. Over the entire duration, the microphone signals were recorded, but only after 20 000 recorded sampling steps, the sine tone was output via the point source. Subsequently, the measured data were highpass filtered to remove frequencies below 500 Hz from the signals. The filtered data was then fed into the inputs of the coupling system (x 1 (t) and x 2 (t)). The measured signals are shown in figures 10(a) and (b) and have a SNR of 54.4 dB and 52.2 dB, respectively. The peak amplitude difference for each angle of incidence is plotted in figure 10(c). The angle of incidence can be determined with a deviation of 0.25 • . The peak amplitude difference of 17.28 dB is lower than for the simulated ideal data (30 dB) and for the simulated noisy data (29.72 dB). The simulated noisy data had the above given SNR of 54.4 dB and 52.2 dB, respectively. For the results shown in figure 11, the point source was placed at a radial distance of 1.0 m from the origin of the microphone array. Angle of incidence and recording time remained unchanged. The measured signals have a SNR of 44.1 dB and 45.0 dB, respectively. In this experiment, an angle of incidence of 26 • was determined (deviation of 1 • ). The peak amplitude difference was 24.17 dB for the measured data and 28.39 dB for the simulated noisy data (SNR of 44.1 dB and 45.0 dB). Due to this reduced peak amplitude differences (17.28 dB and 24.17 dB, respectively) and the resulting flatter curves, the 3 dB beam widths were wider. This was 14.5 • and 6.5 • for the experiments with measured data, respectively.

Discussion
In this work, a bioinspired method of sound source localization is presented which is based on the auditory coupling system of the fly O. ochracea [19]. The method is based on optimization of the system's coupling parameters according to Hou et al [33], but for specific directions of incidence of the sound waves. The optimization was conducted under simplified far field assumptions (plane wave propagation) as well as for near field assumptions (spherical wave propagation) while the latter is the more general case which also allows for amplitude attenuation over distance. Like the fly, this model also uses two sound sensors (microphones). The goal of the parameter optimization was to tune the coupling system for a specific target angle. This resulted in a different set of coupling parameters for each specific target angle in a given interval. In this process, the optimization was carried out in such a way that the model was forced to respond with the same value (amplitude difference of the outputs equals to 30 dB) whenever the model was excited from the direction it was optimized for. To reach this, the optimizer was provided with a nonlinear equality constraint. The optimization process was carried out for sound signals of 1 kHz and a microphone distance of 0.06 m. The radial distance between sound source and microphones was varied in the interval r 0 = {0.2 m, 0.5 m, 1.0 m, 2.0 m}.
In order to evaluate the optimized parameters for their respective direction of incidence, the method was tested on both simulated and physically measured data. In the tests with simulated microphone signals, the coupling system responded with high accuracy (see section 4.1), i.e. 90% of all possible directions of incidence could be determined with an accuracy of less than 1 • (see e.g. figure 9(b)). In comparison, in [41] a sound source localization system for mobile robots is presented that computes an angular error of less than 15 • for 92% of the tested data from a room with noise and reverberation. Four microphones in a cross arrangement are used for this purpose. The opposing microphones have a distance of 0.40 m to each other.
The method consistently responded with the optimized peak amplitude difference of 30 dB for the simulated angle of incidence of the sound wave. In the experiments with the physically measured data, the correct angle of incidence was detected with a lower peak amplitude difference (see section 4.2). During the experiments, it was found that the relative amplitude difference between the measured microphone signals deviates from those modelled according to the theoretical assumptions made in section 2.3 for spherical wave propagation and on which the derivation for the most sensitive angle of incidence is based. Finally, this deviation also causes the reduced peak amplitude difference, which is below 30 dB. Considering figures 10(c) and 11(c) the white noise that is added to the simulated ideal data is not the reason for the discrepancy of the amplitude differences. Both simulated and experimental results showed a small 3 dB beam width. Under the same conditions (two microphones at an inter-microphone distance of 0.06 m), a standard delay-and-sum beamforming algorithm could not even reach values below 3 dB outside the direction of incidence.
In future work, the chosen microphone distance and excitation frequency could be changed and the coupling system optimized accordingly to be sensitive to a wider frequency band. An optimization to the most sensitive frequency could be used for this purpose [33]. The presented method also offers the possibility to be extended into a third spatial dimension by placing more microphones on a twodimensional array. This would result in a higher number of coupled microphone pairs that could be integrated into the bioinspired method. In comparison, a method was developed in [17] to perform 3D sound source localization based on ITD signals of a selfrotating two-microphone array. That method also achieved a high accuracy and was successful in 3D localization of wide-band speech sources and white noise sound sources. The microphones in that setup were separated by a distance of 0.2 m.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.