Student Member, IEEE, Jesper Usama Saqib

Member, IEEE Rindom Jensen

A Framework for Spatial Map Generation using Acoustic Echoes for Robotic Platforms

Keywords: robot audition, TOA estimation, DOA estimation, echolocation, localization, mapping

In this work, we present a framework for constructing a spatial map of an indoor environment using the concept of echolocation. More specifically, we propose a non-linear least squares (NLS) estimator which is combined with a spatial filtering technique, e.g., beamforming, to estimate both the time-of-arrival (TOA) and direction-of-arrival (DOA) of the acoustic echoes. The proposed framework is complemented with an echo detector to classify a spurious estimate and an acoustic reflector, i.e., a wall. Based on these estimators, we then propose two algorithms that can complement existing range sensors and aid a robotic platform in acoustic reflector localization and mapping: a singlechannel localization and mapping (ScLAM) and a multi-channel localization and mapping (McLAM). Compared to commonly used sensors, e.g., lidar and cameras, our proposed method can detect transparent surfaces that are typically found in an office environment. To test our algorithms, a proof-of-concept robotic platform was built. According to our evaluation, both qualitative and quantitative experiments reveal that the proposed methods can detect an acoustic reflector up to a distance of 1.5 m at a signal-to-diffuse-noise ratio (SDNR) of 0 dB in a simulated environment and 10 dB in a real environment with an accuracy of 80 %.

I. INTRODUCTION

Robotic platforms, e.g., drones and unmanned ground vehicles (UGVs), has become an essential part of our society. We use them for tasks that are often monotonous and dangerous for human workers to handle. With the advancement of perception technology, i.e., the ability of a robot to perceive its environment, robots are now able to perform complicated tasks that makes them suitable in different sectors, e.g., agriculture [START_REF] Rovira-Más | Augmented perception for agricultural robots navigation[END_REF], construction [START_REF] Melenbrink | On-site autonomous construction robots: Towards unsupervised building[END_REF], supply chain and logistics [START_REF] Lima | An industry 4.0 approach for the robot@factory lite competition[END_REF], hospitals [START_REF] Huang | Mixed path planning for multi-robots in structured hospital environment[END_REF], etc. Within a warehouse setting, robots are often programmed to follow a predefined trajectory within an environment to transport goods. Over time, robots were equipped with proximity sensors, e.g., lidars, cameras, ultrasonic sensors, etc., for navigation, which also led the robots to plan its own path without human intervention, making these robots more autonomous. According to the IEEE Standard for Robot Map Data Representation for Navigation [START_REF]IEEE standard for robot map data representation for navigation[END_REF], one way to effectively navigate an indoor environment is to construct a spatial map of an environment, which is normally done using a very popular framework called Simultaneous Localization and Mapping (SLAM) [START_REF] Durrant-Whyte | Simultaneous localization and mapping: part 1[END_REF]- [START_REF] Deng | SLAM: Depth image information for mapping and inertial navigation system for localization[END_REF]. The advantage of constructing a spatial map of a surrounding could also aid engineers and building planners to do asset maintenance and survey related work. Additionally, SLAM-based robots also aid rescue workers and surveyors to construct spatial maps of unknown environments, e.g. sewers [START_REF] Wanasinghe | Unmanned aerial systems for the oil and gas industry: Overview, applications, and challenges[END_REF], [START_REF] Worley | Acoustic echo-localization for pipe inspection robots[END_REF], underground tunnels, etc. Traditionally, lidar and camera-based technologies are used to provide input data to SLAM algorithms to construct a spatial map of an environment [START_REF] Hahnel | Map building with mobile robots in dynamic environments[END_REF].

However, lidar and camera-based technologies are susceptible to changing light conditions or foggy environments and are also not suitable for detecting transparent surfaces [START_REF] Saqib | A model-based approach to acoustic reflector localization using robotic platform[END_REF]. This makes these technologies unsuitable to accurately generate a spatial map of a typical office environment [START_REF] Wei | Multi-sensor fusion glass detection for robot navigation and mapping[END_REF]. Furthermore, lidar and camera-based technologies has limited field of view (FOV) and thus offers limited coverage when localizing targets around the corner of the room [START_REF] Hui | Visual SLAM based on EKF filtering algorithm from omnidirectional camera[END_REF]. These issues can be resolved by employing sound [START_REF] Okuno | Robot audition: Its rise and perspectives[END_REF]. Sound is used by animals (e.g., bats, dolphins, and rats) in nature for orienting themselves within an environment and hunting prey [START_REF] Eliakim | A fully autonomous terrestrial bat-like acoustic robot[END_REF] by probing the environment. This process is known as echolocation. An advantage of using echolocation for spatial map generation is that it can enable a robotic platform to navigate an environment under poor lighting conditions. Furthermore, compared to camera and lidar-based technologies, microphones are typically cheaper and may offer omnidirectionality. In the past, the concept of echolocation was studied by several researchers to built active SONAR (Sound navigation and ranging) technologies for naval submarines to detect incoming ships and hostile submarines [START_REF] Muir | Underwater acoustics: A brief historical overview through world war 2[END_REF].

The use of SONAR in air-borne applications is a challenging and complicated task but an attempt to study this was proposed in [START_REF] Kleeman | Mobile robot sonar for target localization and classification[END_REF]. The authors utilizes two ultrasonic transmitters/receivers to effectively localize multiple targets up to a distance of 8 m and classification of the targets were done using template matching. Moreover, the authors in [START_REF] Kreković | EchoSLAM: Simultaneous localization and mapping with acoustic echoes[END_REF]- [START_REF] Boutin | A drone can hear the shape of a room[END_REF] has also proposed several techniques that utilizes sound to make a distance estimate of an acoustic reflector. However, these works assume that the time-of-arrival (TOA) information of the acoustic echoes are known prior to estimation. TOA measurement of an acoustic echo is usually extracted from the estimated room impulse response (RIR) using a standard peakpicking approach, but this has proven to be a non-trivial and time consuming process [START_REF] Kelly | Detecting arrivals in room impulse responses with dynamic time warping[END_REF]. In acoustic signal processing, the RIR is the transfer function between the source and the microphone. It has a distinctive characteristic, i.e., it contains two main components: a direct-path plus early reflections and a stochastic long tail representing late reflections that contributes to the reverberation [START_REF] Moschioni | A new method for measurement of early sound reflections in theaters and halls[END_REF]. The direct-path component corresponds to the shortest distance that a sound travels to reach a receiver while the early reflections correspond to the sound bouncing off an acoustic reflector before reaching the receiver as shown in Fig. 1. Within the context of robotic platforms, the individual RIRs have to be estimated as the robot moves within an environment for TOA estimation of the early reflections as it corresponds to an information about the geometry of the room [START_REF] Baba | 3d room geometry inference based on room impulse response stacks[END_REF]. Furthermore, [START_REF] Steckel | BatSLAM: Simultaneous localization and mapping using biomimetic sonar[END_REF] proposed an algorithm called BatSLAM that utilizes a transmitter/receiver in the ultrasonic frequency range to generate a spatial map of an indoor environment. In [START_REF] Eliakim | A fully autonomous terrestrial bat-like acoustic robot[END_REF], the researchers built a robotic platform that navigates an outdoor environment in order to construct a spatial map as well as classify flora using an artificial neural network. However, both these approach utilizes specialized sensors that operates in a specific frequency range of the sound, i.e., the ultrasonic range. Many robotic platforms that exist in the market are, on the other hand, intended for human-robot interaction (HRI), e.g., the NAO robot. These are often equipped with standard microphones and loudspeakers that operates in the audible frequency range. The estimation of TOAs and DOAs from an observed signal in the audible frequency range has been addressed previously in [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF]- [START_REF] Verhaevert | Direction of arrival (DOA) parameter estimation with the SAGE algorithm[END_REF]. However, the presence of the direct-path component within the recording makes it difficult to detect and estimate early reflections as it contains the highest energy as seen in Fig. 1. Therefore, in [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF], [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF], the direct-path component was assumed to be absent (e.g., removed by preprocessing) from the synthetic data recordings before TOA/DOA estimation, but, in reality, this is a non-trivial problem on its own and could thus have a detrimental effect on the estimation. Additionally, the estimators proposed in [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF], [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF] always gives an estimate regardless of whether there actually is an acoustic reflector in the vicinity of the robot, which may lead to a noisy map generation. To mitigate this, we propose introducing a statistical echo detector within our framework [START_REF] Kay | Fundamentals of statistical signal processing[END_REF]. The detector proposed in this paper is a binary classifier, where the statistics of the background noise is used to optimally define a threshold value against which the spurious estimates are differentiated from those corresponding to actual acoustic reflectors. Furthermore, in order to construct a spatial map of an indoor environment, the direction of arrival (DOA) of the acoustic echo is also required. This helps a robot estimate both the distance and the orientation of the acoustic reflectors, i.e., walls. Several techniques on DOA estimation exist in the literature [START_REF] Valin | Localization of simultaneous moving sound sources for mobile robot using a frequency-domain steered beamformer approach[END_REF]- [START_REF] Maher | Audio forensic gunshot analysis and multilateration[END_REF] to estimate the required information, e.g., TOAs and DOAs, but without addressing the influence of the strong direct-path component of the sound source that could have a detrimental impact in TOA and DOA estimation. To address the presence of the direct-path component, we consider a setup consisting of a microphone array, e.g., a uniform circular array (UCA), and a loudspeaker situated at the center of the array as in Fig. 2. This setup could be placed on different kinds of robotic platforms, e.g., drones, UGVs, etc. Based on these conditions, we propose a framework that extends our existing approach [START_REF] Saqib | A model-based approach to acoustic reflector localization using robotic platform[END_REF], [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF] to multi-channel. The framework consist of four main blocks: 1) a spatial filtering block, that utilizes an adaptive beamformer to filter the observed signals 2) a non-linear least squares (NLS) estimator proposed in [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF] to estimate TOA from the filtered signal 3) an echo detector block, that takes the statistics of the background noise into account to optimally decide a threshold for deciding if the said estimates belong to actual acoustic reflectors and 4) a mapping block to utilize the DOA and TOA to generate a spatial map, as shown in Fig. 3. The process involves probing the environment with a known sound which is recorded by a UCA. The recorded audio data is then processed to estimate the DOAs and TOAs of the acoustic reflections. The advantage of the approach proposed in this paper is that it reduces the influence of the direct-path component resulting from the sound source, e.g., loudspeaker, which is achieved by first processing the acoustic echo with a spatial filter and later using the filtered signal for TOA estimation, and, then, classify whether an estimate belongs to an acoustic reflector or empty space. Finally, the estimates classified as belonging to acoustic reflectors, are used for generating a spatial map. Furthermore, the proposed methods are derived in the frequency-domain which provides a decrease in computational load. The proposed method estimate the parameter of interest directly from the observed signals and does not rely on estimating RIRs as the robot moves. Moreover, to test our proposed framework, we built a proof-of-concept robotic platform to generate real data for a multi-channel scenario similar to our earlier work where we introduced a single-channel estimator for acoustic reflector localization [START_REF] Saqib | A model-based approach to acoustic reflector localization using robotic platform[END_REF]. The dataset is made public 1 .

The remainder of the paper is then structured as follows: Section II describe the signal model and problem formulation, Section III describes the non-linear least squares (NLS) estimator, Section IV describe the first multi-channel localization and mapping (McLAM) algorithm while Section V describes single-channel localization and mapping (ScLAM) algorithm. Additionally, the robotic architecture and components description are detailed in Section VI before proceeding to the Section VII. In Section VIII we test the performance of the proposed method on a robotic platform and then we discuss our findings and conclude the paper in Section IX and Section X. Moreover, in Section X, we propose different ways on how this research could be extended.

II. SIGNAL MODEL AND PROBLEM STATEMENT

A. Time-domain model

Consider an array with M microphones recording a sound emitted from a loudspeaker, including its acoustic reflections from walls, etc. The microphones and loudspeaker are collocated and the loudspeaker is assumed to be a point source. We can then formulate a general model for the recorded signal at microphone m, for m = 1, . . . , M at the kth robot position, as

y m,k (n) = (h m,k * s)(n) + v m,k (n), (1)
= x m,k (n) + v m,k (n)
where h m,k (n) is the acoustic impulse response of the room measured from the loudspeaker to microphone m at robot position, w k , for k = 1 . . . K. Moreover, v m,k (n) is the additive background noise, including interfering sources plus the egonoise of the robot at position, w k . The operator * represents the convolution operator, and x m,k (n) = (h m,k * s)(n). In what follows, the background noise is assumed to be white Gaussian noise, but prewhitening techniques could be employed in cases where such assumptions are not met [START_REF] Christensen | Multipitch estimation[END_REF], [START_REF] Jaramillo | A study on how pre-whitening influences fundamental frequency estimation[END_REF]. By decomposing (1) as a sum of its direct-path component and its reflections and expressing the transfer function between the loudspeaker and a microphone in terms of its gain and delay, the signal model in (1) can be written as:

y m,k (n) = ∞ r=1 g m,r,k s(n -τ ref,r,k -η m,r,k ) (2)
+ v m,k (n)
where g m,r,k is the attenuation of the rth reflection from the loudspeaker to the microphone m at position, w k , while τ ref,r,k is the TOA of the reflected sound received at the reference point of the UCA at robot position, w k , while η m,r,k is the time-difference-of-arrival (TDOA) between the reference point and the microphone m. In our definition in (2), the directpath component corresponds to r = 1. The acoustic impulse response has a certain structure and is distinctively described in two parts: the direct-path plus early reflections and late reflections often described as a stochastic and dense tail. This means that we could rewrite (1) as the sum of the first R reflections to facilitate TOA and DOA estimation as shown

y m,k (n) = R r=1 g m,r,k s(n -τ ref,r,k -η m,r,k ) + d m,k (n) + v m,k (n), (3) = x m,k (n) + v m,k (n) (4)
where d m,k (n) is the stochastic and dense tail of the late reflections. Often, we can combine the late reflections, d m,k (n), with the background noise as shown in (4) [START_REF] Braun | Evaluation and comparison of late reverberation power spectral density estimators[END_REF]. If we collect N time samples from each microphone and assume stationarity across those samples, we can vectorize our data and extend our signal model as shown:

y m,k (n) = R r=1 g m,r,k s(n -τ ref,r,k -η m,r,k ) + d m,k (n) + v m,k (n), = x m,k (n) + v m,k (n) (5) = [y m,k (0) y m,k (1) • • • y m,k (N -1)] T , (6)
where the time-stacked probe signal, s(n), early reflections, x m,k (n), and noise, v m,k (n), are defined similarly to y m,k (n).

Hence, the signal formulation above yield an interesting problem to solve, namely, how to estimate τ ref,r,k and η m,r,k of an acoustic reflector that will aid in simultaneously localizing and mapping an indoor environment. However, this requires us to estimate R unknown TOA and M R TDOAs from the observations y m,k (n), at position, w k . If we assume a known array configuration, however, we can reduce the dimensionality of this problem by incorporating the geometry of the loudspeaker and the microphone array.

B. Array Model

The array model can be chosen to be of any geometry but in this paper, we use a uniform circular array (UCA) with a loudspeaker placed at the center of the array. Although any reference point could be chosen to solve the TOA and DOA problems, we assume the center of the UCA to be the reference point. Assuming that the reflectors are in the far-field of the array and given the geometry of the microphones and the loudspeaker where the center of the microphone array is chosen as the reference point, we can then write the TDOAs of the acoustic echoes as follows:

η m,r (ζ r ) = d sin ψ r cos(φ r -β m ) f s c , (7)
where ζ r = [ψ r φ r ] T , and ψ r and φ r are the elevation and azimuth angles, respectively, while d is the radius of the UCA. Furthermore, β m = 2πi M + α is the angular position of the mth element on the UCA circle counted in an anticlockwise manner from the x-axis and α is the offset angle. Moreover, f s is the sampling frequency and c is the speed of sound. The TDOA model in [START_REF] Huang | A quadratic-complexity observability-constrained unscented kalman filter for SLAM[END_REF] can be combined with the observation model in [START_REF]IEEE standard for robot map data representation for navigation[END_REF] to simplify the dimension of the estimation problem from M R to 2R. The problem of interest is thus to estimate the unknown orientation parameters, i.e., ψ r and φ r , and the distance-related parameter, τ r , based on the posed array [START_REF] Huang | A quadratic-complexity observability-constrained unscented kalman filter for SLAM[END_REF] and observation (5) models. Additionally, a classification of the estimates as either belonging to an actual acoustic reflector or an empty space is needed to generate a spatial map of the acoustic reflectors. Finally, the parameter estimates of the acoustic echoes need to be mapped into the acoustic reflector positions based on the robots movement and orientation.

III. NON-LINEAR LEAST SQUARES (NLS) ESTIMATOR

We can resolve the problem of estimating the unknown parameters in [START_REF]IEEE standard for robot map data representation for navigation[END_REF], i.e., τ ref,r,k and η m,r,k , by using a nonlinear least squares (NLS) estimator, which is statistically optimal under the assumed white Gaussian noise conditions. Mathematically, this can be formulated as

{ g k , τ k , ζ k } = arg min g,τ ,ζ, M m=1 y m,k (n) - R r=1 g m,r,k s(n -τ ref,r,k -η m,r,k (ζ r )) 2 2 , (8)
where

τ = τ 1 τ 2 • • • τ R T , (9)
g = g T 1 g T 2 • • • g T R T , (10)
g r = g 1,r g 2,r • • • g M,r T , (11)
ζ = ζ T 1 ζ T 2 • • • ζ T R T , (12)
with x denoting an estimate of x, and x k denoting a parameter x related to the kth robot position. The displacement k of the robot can be estimated using an accelerometer or can be pre-programmed within the robot so that the robot follows a predefined trajectory. We can also solve [START_REF] Deng | SLAM: Depth image information for mapping and inertial navigation system for localization[END_REF] by converting it into frequency domain because 1) it will reduce the computational load when estimating the desired parameters [START_REF] Avargel | System identification in the short-time fourier transform domain with crossband filtering[END_REF], and 2) by working in the frequency domain, we will have the flexibility to work in specific frequency ranges or account for frequency dependency. For instance, if we want to work in the ultrasonic range then we can select and utilize only the frequency bins corresponding to these high frequencies. This may also help us design probe signals that are nonintrusive to human hearing but this is left for future iteration of this research. Using Parseval's theorem and omitting the frequency dependency in the notation, we can transfer (8) to the frequency domain, which yields the following:

{ g k , τ k , ζ k } = arg min g,τ ,ζ J(g, τ , ζ), (13)
where

J(g, τ , ζ) = M m=1 Y m,k - R r=1 g m,r Z(τ r , ζ r ) S 2 , (14)
Y m,k = Y m,k (0) • • • Y m,k (F -1) T , (15)
Z(τ, ζ) = 1 e -j(τ +η(ζ))2π 1 F • • • e -j(τ +η(ζ))2π F -1 F T ( 16
)
with F denoting the number of frequency bins, Y m,k (f ) denoting the DFT of y m,k (n) in frequency bin f . Moreover, S is the DFT vector of s(n) defined similarly to Y m,k . This estimation problem is multidimensional and thus computationally expensive in practice. To minimize the computational complexity, the multidimensional estimator could instead by implemented using various cyclic methods like the RELAX method proposed in [START_REF] Li | Efficient mixed-spectrum estimation with applications to target feature extraction[END_REF] and later used in [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF] to iteratively estimate the values of τ k and g k . In the special case where we are only concerned with estimating one acoustic reflection, and assuming that the direct-path component has been removed via preprocessing, we can set R = 1. Additionally, if we assume that the gain of each microphone is the same, then we can solve [START_REF] Wei | Multi-sensor fusion glass detection for robot navigation and mapping[END_REF] for the gain g k by taking the derivative of the cost function, yielding:

∂J(g k , τ k ) ∂g k = ∂ ∂g k (Y H Y -g k Y H Z(τ k ) -g k Z H (τ k )Y + g 2 k Z H (τ k )Z(τ k )) = -Y H Z(τ k ) -Z H (τ k )Y + 2g k Z H (τ k )Z(τ k ) (17)
where Z(τ k ) = Z(τ k ) S is the frequency domain probe signal delayed by τ k samples. Solving for the linear gain parameter g k gives:

g k = Y H k Z(τ k ) + Z H (τ k )Y k 2Z H (τ k )Z(τ k ) . (18)
By inserting this back into (13), we get

τ k = arg min τ Y k - Y H k Z(τ ) + Z H (τ )Y k 2Z H (τ )Z(τ ) Z(τ ) 2 (19) = arg max τ IR{Y H k Z(τ )} (20)
where the operator IR represents taking the real part of the signal. The expression in [START_REF] Nguyen | Can a robot hear the shape and dimensions of a room?[END_REF] estimates TOA for a single reflector at position, w k . That is, for the special case with one acoustic echo, the NLS estimator in (20) can be interpreted as a cross-correlation based technique, which is widely used within robotics for source localization [START_REF] Hirata | Real-time ultrasonic distance measurements for autonomous mobile robots using cross correlation by single-bit signal processing[END_REF]. Based on the above problem formulation and methods in Section II and Section III, respectively, we propose two algorithm that could aid different robotic platform for spatial map construction: A multichannel localization and mapping (McLAM) algorthm and a singlechannel localization and mapping algorithm (ScLAM).

IV. MULTI-CHANNEL LOCALIZATION AND MAPPING (MCLAM)

When constructing a spatial map of an environment using sound, a robotic platform requires both DOA and TOA information of the acoustic echoes while distinguishing estimates belonging to an acoustic reflector from spurious estimates. Furthermore, this should be carried out under the presence of a strong direct-path component originating from the sound To address these problems, we propose mounting a microphone array on a robotic platform so that both the DOA and TOA of the acoustic echoes could be estimated, while suppressing the direct-path component. The multi-channel localization and mapping (McLAM) architecture has four important components as shown in Fig. 3. First, we introduce a spatial filter, i.e., a beamformer, to determine the DOAs of the acoustic echoes impinging from the reflectors, e.g., walls. Second, we feed the filtered observation into an NLS estimator to find the TOAs of the acoustic echoes. Then, we introduce a binary classifier to distinguish between spurious and real estimates, to exclude spurious estimates in the subsequent mapping of the acoustic reflectors, which constitutes the final block.

A. Spatial Filter block

The DOA information of an acoustic echo can, for example, be determined using the traditional spatial filtering techniques, e.g., beamforming [START_REF] Pan | Performance study of the MVDR beamformer as a function of the source incidence angle[END_REF], as considered in this paper. Later, a TOA estimate technique is applied, so that acoustic echoes corresponding to the distance of acoustic reflectors are estimated. Apart from DOA estimation, the other advantage of using spatial filtering before TOA estimation is that it can suppress the direct-path component that can affect the parameter estimation. Beamforming is based on the spatial weighting of the signals recorded by a microphone array such that the output signal is the weighted summation of all the signals to extract the signal impinging from a particular DOA [START_REF] Habets | Generating nonstationary multisensor signals under a spatial coherence constraint[END_REF]. In this way, we first employ a beamformer for estimating the angle of an acoustic echo using the steered response power approach. Subsequently, the echo is extracted by applying a beamformer steered towards the estimated angle to produce the output signal for the later TOA estimation using an NLS estimator (20) [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF]. We therefore seek to estimate the signal of interest (SOI), while minimizing the influence of the directpath component of the probe sound and other noise sources, e.g., from the rotors of a drone. With this aim, we consider the use of an adaptive beamformer [START_REF] Pan | Performance study of the MVDR beamformer as a function of the source incidence angle[END_REF]. In addition, this idea also builds on the statistical foundation of the EM method [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF], which indicates that this is the optimal way of solving the problem of localizing acoustic reflectors in an indoor environment.

Due to the broadband nature of the signals involved, we implement the beamformer in the frequency domain. Therefore, the observations in (1) were first converted into frequency

Y k = X k + V k = d(ζ r,k )S r,k + U k (21) = Y 1,k (ω) Y 2,k (ω) • • • Y M,k (ω) T ,
where X k and V k is defined similarly to Y k . Moreover, U k contains the remaining R-1 early reflections as well as the late reverberation and background noise, and S r,k is the complex amplitude of the rth reflection at frequency ω. Assuming a UCA with the center of the array chosen as the reference point, the steering vector can be written as follows:

[d(ζ k )] m = e -j ω c d sin(ψ k ) cos(φ-βm) . (22)
Here, ζ k is the look direction of the beamformer. The objective of the beamformer is then to recover the desired signal S r,k given the observation Y k , i.e.,

Y ζ k = w H Y k , (23)
where w ∈ C M and Y ζ k is the recovered signal from the observed signal from direction ζ k at position w k , which should be an estimate of S k . Here, several beamforming filters could be used, while, in this paper, we consider three types of beamformers which, e.g., facilitates a trade off between computational efficiency, estimation accuracy, and direct-path component suppression: 1) the minimum power distortionless response (MPDR) beamformer, 2) the delay-and-sum (DSB) beamformer, and 3) the linearly constrained minimum variance (LCMV) beamformer [START_REF] Frost | An algorithm for linearly constrained adaptive array processing[END_REF]. The MPDR beamformer is derived by minimizing the power of the of the output of the beamformer Y ζ k subject to a distortionless constraint, i.e.,

w MPDR = arg min w H R Y k w (24) subject to w H d(ζ k ) = 1.
The solution to this is then well known to be given by [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF] w where

MPDR = R -1 Y k d(ζ k ) d H (ζ k )R -1 Y k d(ζ k ) , (25)
R Y k = E[Y k Y H k ] is the M × M covariance matrix of the observed signal, E[•]
is the mathematical expectation operator, and w MPDR is the complex weight vector corresponding to the MPDR beamformer. If the observed signal is assumed to be white Gaussian noise, e.g., R Y k = I M , where I M is the M × M identity matrix, the MPDR design resembles the DSB, i.e.,

w DSB = d(ζ k ) M . (26)
Similarly, the LCMV beamformer is derived by extending the MPDR beamformer with additional constraints such that the optimization problem is solved as shown:

w LCMV = arg min w H R Y k w (27) subject to w H D = f T .
Here, D is a matrix containing all the steering vector for the C different constraints in f ∈ R L . In this paper, we choose f as

f = 1 0 • • • 0 T (28)
By utilizing this, we can reject the interference of the directpath component by introducing a null in the direction of the loudspeaker, i.e., the center of the UCA. The solution to the LCMV beamforming problem is

w LCMV = R -1 Y k D[D H R -1 Y k D] -1 f . (29)

B. TOA estimator block

The output of these beamformers are subsequently feed to the non-linear least squares (NLS) Estimator for TOA estimation in [START_REF] Nguyen | Can a robot hear the shape and dimensions of a room?[END_REF]. This estimator is statistically optimal when estimating τ and g for a single reflection while the background noise is white Gaussian. By preprocessing the observation with the adaptive beamformer, this assumption is better met since we can reduce the impact of directional and colored noise [START_REF] Schwartz | Multispeaker LCMV beamformer and postfilter for source separation and noise reduction[END_REF]. The resulting NLS estimator is then given by

{ g k , τ k , ζ k } = arg min g,τ,ζ Y ζ -gZ(τ ) S) 2 , ( 30
)
where Y ζ is the output of the beamformer (23) extracted from direction ζ at a position w k at frequency ω, while is the element-wise multiplication operator. By solving for the linear parameters in [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF], we get the concentrated estimator for the TOA and DOAs:

{ τ k , ζ k } = arg max τ,ζ IR Y H ζ Z(τ ) (31)
where Z(τ ) = Z(τ ) S and the operator IR represents taking the real part of the signal.

C. Echo detector block

If the robotic platform is expected to move autonomously based on echolocation, a problem of great importance is to detect whether the observed signal received by the microphones represent an acoustic reflector, or if it only contains noise, e.g., the ego-noise of the robotic platform. This is because the TOA estimator in [START_REF] Kay | Fundamentals of statistical signal processing[END_REF] provides estimates even when no acoustic reflector is present, which may lead to spurious localization estimates. To prevent false estimation when no acoustic reflector is present, several approaches could be applied including machine learning approaches [START_REF] Liang | Machine learning-based 3-d geometry reconstruction and modeling of aortic valve deformation using 3-d computed tomography images[END_REF], and deep learning [START_REF] Yu | Room geometry estimation from room impulse responses using convolutional neural networks[END_REF] to categorize acoustic reflectors. Another approach could be to include a Generalized Likelihood Ratio Test (GLRT) detector [START_REF] Kay | Fundamentals of statistical signal processing[END_REF] within our framework to distinguish whether the observed signal contains an acoustic reflection or not. Compared to the data-driven machine learning approaches, the GLRT is based on a priori model assumptions, and does thus not require training data. In this paper, we therefore employ the GLRT detection approach as discussed in the following.

If we assume the acoustic reflection to be in the far-field of the array, the decision about whether an observation contains an acoustic reflection can be formulated as a detection problem [START_REF] Kay | Fundamentals of statistical signal processing[END_REF]:

H 0 : y m,k (n) = v m,k (n) (32) H 1 : y m,k (n) = g k s(n -τ m,k ) + v m,k (n), (33)
for m = 1, . . . , M , where H 0 is the null hypothesis referring to a situation when the observation only includes white Gaussian background noise and late reverberation, v m,k (n), with variance σ 2 , while H 1 refers to the situation when the observation includes a reflected version of the known probe signal s(n) in noise. Here, we assume that the direct-path component is absent, i.e., suppressed via preprocessing. The GLRT is then given by

p(y k ; g k , H 1 ) p(y k ; H 0 ) > γ, (34)
y k = y T 1,k (0) • • • y T M,k (0) T , (35)
It can then be shown that, in order to detect if the observation belongs to H 1 , we can use a threshold that depends on the power of the attenuated probe signal, the noise variance, and γ.

If the power, T (y k ), of a matched filtering between the probe signal and the observed signal at the reference microphone exceeds this threshold, we decide H 1 , i.e., if

T (y k ) = y H k (n)H(τ k )s(n) > g k 2 + σ 2 ln γ 2 g k (36)
with

H(τ k ) = D T τ 1,k • • • D T τ M,k T , (37)
= M s(n) 2 , (38)
g k = 2y H k (n)H(τ k )s(n) M s(n) 2 , ( 39
)
where D τ is a cyclic shift register that delays a signal by τ samples. The GLRT derivation is shown in Appendix ??.

D. Mapping block

In this block, the DOA and TOA estimates are used alongside the robot's position within an environment to localize the position of an acoustic reflector. The aspect of the robot's navigation and path planning is beyond the scope of this paper, however, by utilizing common on-board sensors, e.g., Initial Mass Units (IMUs), of the robotic platform, we can estimate the robot's position. By combining this information with the estimates of the acoustic echoes obtained using, e.g., the methods considered in this paper, a spatial map of the environment can be generated for the robotic platform. The resulting spatial map may then enable the robotic platform to plan its path and move autonomously within the environment.

To estimate the position of the acoustic reflector from the estimated TOA, τ k , we assume that the sound propagates in plane waves (i.e., the source is in the far-field of the array). If we assume the speed of sound to be fixed then also assume estimation of a single acoustic reflector then the distance of the acoustic reflector with respect to the robotic platform is estimated as

δ k = c•τ k 2 .
Additionally, the direction of the acoustic reflector at position w k is determined from the DOA estimates ψ and φ. If we assume a 2D scenario, where the reflections and the hardware is located in the same plane, we can utilize the far-field assumption and the choice of our reference point to conduct the mapping as:

p x k = w x k + δ k cos φ k (40) p y k = w y k + δ k sin φ k
The procedure is then to estimate the acoustic reflector positions for each of the known robot positions, w k , along its trajectory. The estimated acoustic reflector positions are then concatenated in the set

P = {p 1 , . . . , p K } with p k = (p x k , p y k ) for k = 1, . . . , K.
The spatial filtering, the TOA estimator, the echo detector and the mapping block are then combined to form the basis of our proposed McLAM method. The algorithm describing the proposed McLAM method is outlined in Algorithm 1. However, in some applications, only one microphone and loudspeaker pair may be available for the mapping. In the following section, we therefore consider, how the hardware directivity properties may be exploited to localize the acoustic reflectors.

V. SINGLE CHANNEL LOCALIZATION AND MAPPING

(SCLAM) In some applications, robotic platforms, e.g, intended for HRI may consist of only a single loudspeaker and microphone. In such a scenario, it is therefore necessary to reduce

(n), y k (n) FFT ---→ S, Y k ; for φ ∈ Φ do
Compute w, e.g., using [START_REF] Baba | 3d room geometry inference based on room impulse response stacks[END_REF];

Y φ,k (ω) = w H Y k ; end { τ k , φ k } = arg max τ,φ IR Y H φ,k Z(τ k ) ; φ k update ----→ DOA; τ k update ----→ TOA;
Apply the echo detector in [START_REF] Christensen | Multipitch estimation[END_REF];

if y H k (n)H(τ k )s(n) > g k 2 + σ 2 ln γ 2 g k then Compute p k using (40); p k update ----→ P; end end
the McLAM algorithm to a single channel localization and mapping (ScLAM) algorithm. However, using such a singlechannel approach has certain limitations, for instance, it cannot generally not be used to estimate the DOA of the acoustic echoes because of the lacking spatial information. Some possible ways of combating this are to, e.g., exploit the movement of the robot [START_REF] Kreković | EchoSLAM: Simultaneous localization and mapping with acoustic echoes[END_REF], or, as considered in this paper, to exploit the directionality of the employed hardware [START_REF] Saqib | A model-based approach to acoustic reflector localization using robotic platform[END_REF].

As in the McLAM, the loudspeaker probes the room with a known sound, s(n), which is recorded by a microphone as the robot moves via positions w k , for k = 1, . . . , K. The NLS estimator described in [START_REF] Nguyen | Can a robot hear the shape and dimensions of a room?[END_REF] estimates τ k for every robot position, w k . Consider the platform moving in a predefined trajectory W = {w 1 , . . . , w K } with w k = (w x k , w y k ), such that the platform moves from w k to w k+1 etc. For every position, w k , the platform will thus probe the environment with s(n) and record the observed signal y k (n). The probed and observed signals are then converted into the frequency domain before passing them to the NLS estimator. In practice, the analysis window for the TOA could be restricted to a search interval from τ min up to τ max samples. This leads to

τ k = arg max τ [τmin;τmax] IR{Y H k Z(τ )} (41)
In ScLAM, the position of the acoustic reflector is then inferred from the estimated TOA, τ k , by exploiting the typical directionality of a loudspeaker. More specifically, we assume the acoustic reflector to be located at the distance correspond-

(n), y k (n) FFT ---→ S, Y k ; τ k = arg max τ k IR{Y H k Z(τ )}; { τ k } update ----→ TOAs;
Apply the echo detector in; [START_REF] Christensen | Multipitch estimation[END_REF];

if y H k (n)H(τ k )s(n) > g k 2 + σ 2 ln γ 2 g k then τ k remove ----→ TOAs;
p k using [START_REF] Frost | An algorithm for linearly constrained adaptive array processing[END_REF] update ----→ P; end end ing to the estimated τ k in the direction of the loudspeaker. Additionally, the direction in which the robot platform is moving, θ rob,k , at position w k , is related to the direction that the loudspeaker is facing, θ l k , by a fixed offset angle, ∆θ, i.e.,

θ l k = θ rob,k + ∆θ. ( 42
)
Based on the above information, the coordinates of the position of the acoustic reflector is then estimated as follows:

p x k = w x k + δ k cos θ l k , (43)
p y k = w y k + δ k sin θ l k .
The resulting ScLAM algorithm is then proposed in Algorithm 2, which can be used used to construct a spatial map of a 2D environment with a single-channel loudspeaker/microphone setup.

VI. ROBOTIC PLATFORM OVERVIEW

The proposed methods discussed in Section IV and Section V were implemented on an embedded platform running a Windows 10 Operating System. The microcomputer used for the proof-of-concept robotic platform is an UDOO x86, which is a single board development platform. On the platform, we use MATLAB to implement the proposed McLAM and ScLAM methods in Algorithm 1 and 2, respectively. Moreover, for multichannel audio data acquisition, Playrec [START_REF] Humphrey | Playrec: Multi-channel matlab audio[END_REF] was used to probe and record the acoustic signals. The base of the robot used for moving the microphone and loudspeaker array as shown in Fig. 4 is a Kobuki (TMR-K01-W1), which is a wheeled platform with on-board sensors such as accelerometer, odometer, etc., for precise control and movement. The Kobuki platform has a built-in microcontroller (Arduino) that can be programmed with a predefined trajectory to conduct experiments. The microphone and loudspeaker array is connected to a Presonus (1818VSL) audio interface, which was subsequently connected to the UDOO x86 microcomputer. The sampling frequency of the audio interface was set to 48, 000 Hz. Furthermore, a pre-calibrated laser range sensor (TFMini micro Lidar), was also attached to an external microcontroller (Arduino Uno) which was then connected to the UDOO microcomputer to receive a ground truth distance value for the experiments. The laser range finder helps in evaluating the performance of the proposed method at varying distances under different noise conditions. The recorded data was processed by the UDOO x86 microcomputer in real-time as the robot was moving along its trajectory. The final assembly is shown in Fig. 5 where the microphone and loudspeaker array is attached on top of the Kobuki base. The microphones are organized as a UCA with a radius of 0.2 m.

VII. SIMULATED RESULTS

In this section, we evaluate the performance of the proposed method presented in the earlier sections. We evaluate the performance of the ScLAM and McLAM using simulation data and later, implement the methods on a proof-of-concept hardware platform that was built to test the proposed method in a lab setting. In the first experiments, the performance of the considered TOA/DOA estimators in terms of their accuracy were evaluated and compared against existing methods under different background noise levels. Similarly, in the second experiment, we evaluate the TOA/DOA accuracy of the proposed methods against varying distances from the acoustic reflector. The simulated experiments were conducted using the room impulse response generator [START_REF] Habets | Room impulse response generator[END_REF]. The dimension of the simulated room was set to 8 × 6 × 5 m., the reverberation time (T 60 ) was set to 0.6 s, while the speed of sound was fixed at 343 m/s. The loudspeaker was positioned at the center of an UCA with a radius of 0.2 m and M = 6 microphones. A white Gaussian noise sequence was used as the known probe signal, s(n), consisting of 1, 500 samples from a Gaussian distribution. Using such a broadband signal minimizes the effect of spatial aliasing [START_REF] Dmochowski | On spatial aliasing in microphone arrays[END_REF] and was also used in [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF] to simplify the EM estimator. However, any type of known broadband signal could be used to probe the environment, e.g., a chirp signal or a maximum length sequence (MLS) [START_REF] Florencio | Maximum a posteriori estimation of room impulse responses[END_REF]. Additionally, we have zero-padded the probe signal to get a total length of 20, 000 samples so that we get a longer analysis window which will ensure that all the reflections are captured in the observed signal. The sampling frequency f s was set to 48, 000 Hz. The background noise for the evaluation was composed of three components: a cylindrical diffuse noise e m,k , the sensor noise, f m,k , and an interfering source, i m,k , e.g., external and directional noise source. The diffuse cylindrical noise was generated using the method in [START_REF] Habets | Generating nonstationary multisensor signals under a spatial coherence constraint[END_REF] with the rotor noise of a drone from the DREGON database [START_REF] Strauss | DREGON: Dataset and methods for UAV-embedded sound source localization[END_REF]. The audio file used to generate the cylindrical noise has a rotor speed of 70 revolutions per second (RPS). The thermal sensor noise was simulated as a white Gaussian noise while the interfering source is modelled as a point source. These noises were then added to the observed probe signal before

y m,k (n) = x m,k (n) + v m,k (n), ( 44
) = x m,k (n) + e m,k (n) + f m,k (n) + i m,k (n). ( 45
)
The noise was added to achieve certain signal-to-diffuse noise ratios (SDNR's), signal-to-sensor noise ratios (SSNR's), and signal-to-inteference noise ratios (SINR's). These are defined, for the microphones m = 1, . . . , M , as

SDNR m = σ 2 xm σ 2 em , (46)
SSNR m = σ 2 xm σ 2 fm , ( 47
)
SINR m = σ 2 xm σ 2 im , ( 48
)
where σ 2 y denotes the variance, σ 2 y = E[y 2 (n)] of a zero-mean signal y(n). In the following experiments, we then compare our proposed method with existing TOA/DOA methods found in the literature. This includes, the multi-channel expectationmaximization method (EM-MC) method proposed in [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF] and the common approach to extracting TOAs from the estimated RIR using dual-channel method [START_REF] Herlufsen | Dual channel FFT analysis (part I)[END_REF] through the peak-picking approach (RIR-PP). This is done by computing H(f ) = Y (f )/S(f ) and then taking the inverse DFT to get h = F -1 { H(f )}. These methods were compared with different variations of the proposed beamforming and NLSbased approach, utilizing DS (DS-NLS), MPDR (MPDR-NLS), and LCMV (LCMV-NLS) beamforming, respectively.

Although the proposed methods can be extended and applied to 3D scenarios, we focus on the construction of 2D maps in our experiments and therefore set ψ = 0. The generalization to 3D is left for future research. In contrast to earlier works in [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF], [START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF], the direct-path component is accounted for and thus included within the simulations. Within the experiments, we assume that the robotic platform is closer to one acoustic reflector. Therefore, we choose R = 1 to estimate the TOA and the DOA of the nearby acoustic reflector. In order to estimate multiple reflections R > 1, we can adopt several iterative methods, such as, RELAX and EM method [START_REF] Li | Efficient mixed-spectrum estimation with applications to target feature extraction[END_REF], [START_REF] Feder | Parameter estimation of superimposed signals using the em algorithm[END_REF]

A. Implementation of the proposed DOA estimator

To implement the beamformers, we use the overlap-add technique [START_REF] Pan | Performance study of the MVDR beamformer as a function of the source incidence angle[END_REF]. The output of the microphone was divided into overlapping frames with a frame width of 960 samples (20 ms with a sampling rate of 48 kHz) with a window overlap of 50 %. Later, each frame is multiplied with a Hanning window. These frames are then transformed using a short-time Fourier transform (STFT). For each frequency bin, a beamformer was designed and applied to the received signals Y k . Furthermore, for each sub-band, the observed signal covariance matrix, needed in forming the MPDR and LCMV beamformers, is estimated as

R Y k = 1 T T -1 t=0 Y k Y H k . ( 49
)
Moreover, to make the beamformers robust against, e.g., miscalibration and reverberation, we regularize the covariance matrix of the observed signal as in [START_REF] Benesty | Noise reduction in speech processing[END_REF]

R Y k = (1 -β)R Y k + β Tr{R Y k }I M ( 50
)
where β is the regularization parameter, Tr(•) is the trace of a matrix, and I M is the M ×M identity matrix. When evaluating the performance of our estimator, a value of β = 0.1 was selected for the MPDR beamformer. The noise covariance matrix, R Y k , in ( 49) is then replaced by regularized noise variance matrix [START_REF] Florencio | Maximum a posteriori estimation of room impulse responses[END_REF], R Y k . For the LCMV beamformer, we added an additional regularization using γ to mitigate poor matrix conditioning for certain constraint and frequency combinations. This was done as

w LCMV = R -1 Y k D[A(γ)] -1 f , where A(γ) = (1 -γ)D H R -1 Y k D + γ Tr{D H R -1 Y k D}I M . ( 51
)
Values of γ = 0.1 and γ = 1 were selected empirically and used in the simulations. To initiate the method, we probe the environment with a known sound. The observed signals recorded by the microphone array were first processed to determine the DOA of the acoustic echoes. To estimate the DOA and the TOA of the acoustic echoes, a uniform grid of DOAs over the interval [0°; 360°] and a uniform grid of TOAs corresponding to a distance interval from 0.5 m up to 3 m were considered. The estimators were then evaluated over these grids of candidate DOAs and TOAs. The reason for selecting 0.5 m as the lower bound was done to search for acoustic echoes that are outside the UCA which has a radius of 0.2 m, and so that the direct-path component is not included within the estimation window. Moreover, the upper bound of 3 m were selected because the performance of the proposed method degrades after 3 m according to [START_REF] Saqib | Estimation of acoustic echoes using expectation-maximization methods[END_REF].

B. Comparison of the proposed methods

In our first experiment, we compare our proposed method with the existing TOA/DOA methods. We compare the proposed methods against different SDNRs while placing the setup at a distance of 1 m close to an acoustic reflector. The performance of the proposed methods are shown in Fig. 6. The accuracy is defined as the percentage of estimated TOAs that are within ±10 % of the true TOA/DOA parameter of the first order acoustic echo computed using the image-source method [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF]. This was measured for different SDNRs while the SSNR was fixed to 40 dB and the interfering source was absent in this experiment. For each SDNR values the accuracy was measure over 50 Monte-Carlo simulations. As seen in Fig. 6, the proposed methods, MPDR-NLS and LCMV-NLS, outperforms the existing TOA/DOA methods, EM-MC and RIR-PP, for SDNR levels greater than -10 dB. The DSB-NLS method offers similar performance to EM-MC both in terms of TOA and DOA estimation for most SDNRs as seen in Fig. 6(b).

C. Evaluation of the proposed method in the presence of a point source interference

In this experiment, we investigate a scenario where the robot is placed within an environment in the presence of an external interfering source, e.g, a human-speaker, machinery, a radio, etc. In such a scenario, the proposed method will be affected from external elements present in the environment. Therefore, the objective of this experiment is to evaluate both the TOA and the DOA performance of the proposed method against different SINR values. The interfering source was modelled as a point source for this experiment. More specifically, within this experiment, the robotic platform was placed close to an acoustic reflector at a position, [1, 3, 2.5] m within an environment of dimension 8 × 6 × 5 m 3 . Furthermore, the external interfering point source was positioned at a location [2, 1, 2.5] m such that the acoustic reflector is at a fixed angle of 180°while the point source is placed at an angle (a) TOA evaluation of 300°with respect to the robotic platform. The performance is shown in Fig. 7. The SINR level selected for this experiment are within the interval [-40; 40] dB while SDNR and SSNR were both set to 40 dB. Moreover, some additional consideration was taken into account when modelling the interfering point source. For instance, if a human talker is considered as a point source, then it is natural for the human to move within the environment. To model this, the position of the point source was randomize in both the x-axis and y-axis.

The interval selected to model point source movement for both x-axis and y-axis are [1; 3] m and [1; 2] m, respectively. As seen in Fig. 7, the TOA of the MPDR-NLS and LCMV-NLS offer more robustness at low SINR compared to EM-MC, RIR-PP and DSB-NLS. Similar performance is seen in the DOA estimation. The accuracy is defined similarly to the previous method with a tolerance of ±10% of true TOA and DOA.

D. Evaluation of proposed methods against distance

In this experiment, we consider a scenario where the robotic platform is placed closer to an acoustic reflector and its distance with respect to the acoustic reflector was changed after every 50 iteration. With this setup, the performance of the proposed method and existing methods over distance interval [0.8; 2.2] m was investigated. Here, the SDNR and SSNR values were set to 40 dB while the interfering source was absent. As seen in Fig. 8, the MPDR-NLS, and the LCMV-NLS variants outperform other methods in terms of TOA estimation and accurately estimate the DOA of the acoustic reflector as it can detect an acoustic reflector up to a distance of around 2 m. This is because at larger distance the acoustic echoes loses its energy quadratically due to inverse square law.

E. Visualizing Acoustic Echoes

Microphone array imaging has been around for quite sometime and are used in aviation [START_REF] Legg | Automatic 3d scanning surface generation for microphone array acoustic imaging[END_REF] for structural analysis as well as to study low frequencies [START_REF] Williams | Sound source reconstructions using a microphone array[END_REF]. Similarly, our proposed method could also be used to generate an acoustic image of acoustic echoes which could aid researcher to analyse the direction and distance of acoustic reflectors or be used as input data for the development of deep learning based methods.

To generate an acoustic image using the methods proposed in this paper, we consider estimation of the reflector in 2D only, i.e., we only estimate φ. For each beamformer we considered a grid of candidate steering angles with a resolution of 4°in the interval [0°; 360°]. The output of the beamformer is then passed to the NLS estimator in [START_REF] Nguyen | Can a robot hear the shape and dimensions of a room?[END_REF], which then estimates τ from candidate grid of delays in [τ min ; τ max ]. The resulting 2D cost functions are shown in Fig. 9 including the direct-path component. As seen, the cost function of the MPDR beamformer Fig. 9(a) shows a peak at times and angles corresponding to the TOAs and the DOAs at which the beamformer received the acoustic echo, these regions are marked by a red circle. In comparison, the DSB cost function in Fig. 9(b) is very noisy, despite evaluating under the high SINR of 40 dB, which makes it difficult to extract the true TOAs and DOAs. This is partly caused by the presence of the direct-path component and the ego-noise, which cannot be sufficiently suppressed by the DSB.

LCMV-NLS γ = 1.0 LiDAR [m] µ [m] RMSE [m] µ [m] RMSE [m] µ [m] RMSE [m] µ [m] RMSE [m] µ [m] RMSE [m] 1.

F. Computational cost

The computational cost of the proposed methods were measured using MATLAB's built-in function timeit. These were tested on a standard desktop computer running a Microsoft Windows 10 operating system with an Intel Core i7 CPU with a 3.40 GHz processing speed and 16 GB of RAM. A Monte Carlo Simulation of 50 trials were conducted and an average time was calculated. The measured computational time of EM-MC, RIR-PP, LCMV-NLS and MPDR-NLS are 63.25 s, 0.024, 59.75, and 60.65 s, respectively, for R = 1 for SINR = 40 dB. The proposed algorithms are computationally expensive when implemented within a robotic platform compared to lidar technologies. To address these issues, one tweak that would allow faster computation on the robot is to probe the environment and use echo detector first to determine whether the robot is closer to an acoustic reflector before proceeding with the proposed DOA/TOA estimates. This accelerate processing and prevents the robot from estimating the parameters when not in the presence of an acoustic reflector. Moreover, tracking, e.g., in the form of gradient searches, may be employed instead of performing a full grid search for every new robot position.

VIII. EXPERIMENTS USING PROOF-OF-CONCEPT ROBOTIC PLATFORM

In this section, we evaluate the performance of the proposed algorithm (McSLAM) using a robotic platform under different SINRs and distances. The objective of these experiments are to compare our simulated data with the real data to test the performance of the proposed method in real scenarios. Two sets of experiments were conducted using the proof-of-concept robotic platform described earlier in Section VII. The first set of experiments were performed under different SINR and distances while the second set of experiments were performed as a qualitative test to show the mapping ability of the robotic platform while comparing the MPDR-NLS algorithm against the lidar data (ground truth). The data is also summarized in Table I andTable II. A. Evaluation of the proposed method against different SINRs and distances Similar to the experiments performed in Section VII-B, the proof-of-concept robotic platform was placed against an acoustic reflector within Aalborg University's Sound Lab. In the first part of the experiment, the robotic platform was placed at varying distances while the SINR value was set fixed to 40 dB. The platform was placed at an interval of [1, 1.5, 2] m. At each distance, the robotic platform probed the environment with a known sound and 50 samples were collected at each distances. The TOA/DOA obtained from the robotic platform are shown in Fig. 10. As seen in the figures, the proposed McLAM algorithm gives an accurate TOA estimate up to a distance of 1.5 m for all combination of spatial filter. The accuracy is defined as the number of estimates that are ±10% of the true TOAs obtained from the lidar data. DOA accuracy is defined similarly.

The next experiment was performed to test the proposed method against different SINR values of the environment. The SINR value of the environment was changed by using a separate loudspeaker playing an audio file from YouTube called Cocktail party 2 . The loudspeaker was placed 6.3 m away from the robotic platform while the robotic platform was fixed at a distance of 1 m away from the acoustic reflector. The SINR of the environment was estimated by dividing the variance of the probed signal, σ 2

x , with the variance of the background noise, σ 2 v . The background noise v(n) was recorded by the robot before probing the environment. By tuning the volume of the loudspeaker, we then select 5 SINR values, [0, [START_REF] Worley | Acoustic echo-localization for pipe inspection robots[END_REF][START_REF] Nguyen | Can a robot hear the shape and dimensions of a room?[END_REF][START_REF] Saqib | Sound-based distance estimation for indoor navigation in the presence of ego noise[END_REF][START_REF] Li | Efficient mixed-spectrum estimation with applications to target feature extraction[END_REF] dB. The results for this experiment is shown in Fig. 11. Here, we see that the proposed MPDR-NLS is robust under low SINR value of 10 dB for both TOA and DOA estimation with 80% accuracy. The changes seen in this experiments are discussed in Section IX.

B. Application Examples

Two qualitative experiments were performed to test the performance of the proposed method (MPDR-NLS) in constructing a spatial map of an indoor environments. Two environments were selected to perform this task: 1) a typical office environment with a glass partition and 2) Aalborg University's Sound Lab. These experiments are similar to the one performed in our earlier work with ScLAM [START_REF] Saqib | A model-based approach to acoustic reflector localization using robotic platform[END_REF]. In the first experiment, the McLAM algorithm was used to move within an office environment in a predefined trajectory (straight line). The objective of this experiment was to compare the proposed method against lidar, e.g., in detecting a glass surface. The robot moves a distance of 0.5 m and stops momentarily to probe the environment with a known sound before moving to a new location. The robot repeats this process for k = 1, . . . K, positions. The results are shown in Fig. 12. As seen from the experiment, the proposed method is capable of detecting a glass surface compared to the commonly used lidar sensor. This shows that the proposed method is suitable for constructing a spatial map of a typical office environment. https://youtu.be/IKB3Qiglyro

In the second experiment, Algorithm 1 was used within Aalborg University's Sound Lab, which has a dimension of 5.4 × 6.38 × 4.05 m 3 , to construct a spatial map. The objective of this experiment was to move the robot in a more elaborate path within a 3D space such that the robot encounters acoustic reflectors as well as empty space along its trajectory. This was done to construct a spatial map of an enclosed environment and also to test the echo detector method presented in Section IV-C. To accomplish this task, the room was divided in to a grid of 20 square boxes, each box has a size of 1 m 2 so that we can ensure that the robot moves along its predefined trajectory and robot's location with respect to the acoustic reflector is always known. Autonomous navigation is also possible but this would require additional on-board sensors, e.g., using Initial Mass Unit (IMU), odometer, gyroscope, etc., to estimate the robot's current position which can then be combined with our estimates to generate a spatial map. As the robot moves within the square grids and follows a predefined trajectory as shown in Fig. 13(b), the robot probes the environment with a known sound. The recorded sound is spatially filtered using MPDR beamformer which is later feed to a NLS estimator for TOA estimation. Later, the estimated data is passed to a echo-detector to determine whether it belongs to an acoustic reflector or is an spurious estimate. Finally, the estimated data are combined with the trajectory of the robotic platform to localize acoustic reflectors. As seen in Fig. 13(b), if the robot moves without the echo detector then it will estimate spurious estimates even when the robot is away from any reflecting surfaces. However, these spurious estimates are removed when echo detector is applied as seen in Fig. 13(c) IX. DISCUSSION In the experimental section, the performance of the different methods were evaluated in both simulated and practical environment. According to the simulation results in Fig. 6 and Fig. 8, the MPDR-NLS, LCMV-NLS γ = 0.1 and LCMV-NLS γ = 1 methods detects an acoustic reflector up to a distance of around 2 m under the SDNR of -10 dB. Similarly, in practical scenario, the methods could detect an acoustic reflector up to a distance of around 1.5 m as seen in Fig. 10. However, in Fig. 11, only the MPDR-NLS is seen to provide good accuracy at low SINR for TOA/DOA estimation while LCMV-NLS γ = 0.1 is the second best choice for TOA estimation compared to its other variant LCMV-NLS γ = 1 which performs less then EM-MC and DSB-NLS when evaluating under different SINR. From these experiments, we can deduce that the MPDR-NLS estimator provides better performance compared to other methods. The results from practical experiments are also detailed in Table . II. The RMSE of all beamformer variants are robust when the distance of the acoustic reflector is less than 1.5 m with respect to the robot. At higher distances, the RMSE increases while the RMSE decreases with higher SDNR values. One noticeable difference we see between simulated and practical evaluation is the low accuracy in practical experiment. There could be various reasons for a lower accuracy in real scenarios, for instance, in our proposed method, we do not take sensor calibration into account as well as sensor drift due to an increase in thermal temperature that can arise when the sensors are operating for a long period of time [START_REF] Kim | Acoustical-thermal noise in a capacitive mems microphone[END_REF].

Additionally, in the simulation, we see that the presence of a point interfering source does not limit the proposed methods' robustness as seen in Fig. 7. This goes to show that an MPDR-NLS and both implementations of the LCMV-NLS method are effective in localizing the acoustic reflector of the environment when compared to other beamformer variants in the presence of an interfering source. In the qualitative experiments, we exploit robot's movement to construct a spatial map of an environment. Here, we compare the current technologies, e.g., lidar, with the proposed McLAM algorithm. As seen in Fig. 13, our proposed method successfully construct a spatial map of an environment. However, one obvious limitation of the proposed method is that lidar accurate distance measurements over longer distance. This is because sound intensity decreases quadratically over distance due to the inverse squares law. One major advantage of our proposed method, on the other hand, is that it provides is that it can be used to detect transparent surfaces as seen in Fig. 12 that are typically found in an office environment, hence our proposed method could compliment existing technologies for spatial map generation. Additionally, we also test our echo-detector in our qualitative experiment. As seen in Fig. 13(b), without echo-detector enabled, spurious estimates are seen when the robotic platform is at an empty space. However, as seen in Fig. 13(c), with the echo-detector enabled, the spurious estimates are removed.

Moreover, both variant of LCMV beamformers behave differently using real data and offers lower performance compared to the simulated results. This could be due to mismatch of the microphones/loudspeaker positions in the array that leads, in which case the null constraint of the LCMV beamformers are not aligned with the direct-path component. Moreover, the LCMV beamformer implicitly assumes the loudspeaker to be a point source, which will not hold for larger loudspeakers located close to the array in practice. In our future work, we plan to incorporate these inaccuracies within our models and methods to improve their robustness.

X. CONCLUSION

In this paper, we proposed a non-traditional method of constructing a spatial map of an indoor environment using the concept of echolocation. Instead of working in the ultrasonic range, a novelty in this paper is that our proposed method could work in the audible frequency range and any kind of broadband probe signal could be used for DOA/TOA estimation. In addition to this, our proposed algorithms could be used on existing robotic platform, e.g., NAO robots, as these consist of microphones and loudspeaker that work in audible frequency range. As seen from the experimental results, our proposed framework could utilize different beamformers for DOA estimation, which could be combined with our NLS estimator for TOA estimation and an echo detector, to construct a spatial map. One obvious advantage of our framework is that each modules in Fig. 3 could be separately improved over time in order to increase the performance of the acoustic echo localization. Our proposed methods can detect acoustic reflectors up to a distance of 1.5 m at an SINR of 40 dB and robustly estimate TOA at an SINR of 10 dB with 80% accuracy in realistic scenario. Furthermore, an echo detector is proposed based on the statistics of the background noise that could aid a robot in classifying estimates stemming from acoustic reflectors from spurious estimates. This can help a robot map and environment to facilitate its autonomous planning and movement. In the qualitative experiments, we see that, compared to the commonly used lidar technology, our proposed method can detect transparent surfaces as seen in Fig. 12 and it can also construct a spatial map of an indoor environment as seen in Fig. 13.

In a future iteration of this research, we aim to include a loudspeaker's directivity and transfer function within the signal model. This will enable our algorithms to work more efficiently and help us understand and develop sophisticated sound propagation model that could enable even more accurate construction of spatial maps in an indoor environment. Additionally, we aim to decrease the computation load of the proposed method to make it run faster on embedded devices.

Fig. 1 :

 1 Fig.1: An example illustrating synthetic Room Impulse Response (RIR)

Fig. 2 :

 2 Fig. 2: Example of a uniform circular array with six microphones.

Fig. 3 :

 3 Fig. 3: Proposed McLAM Architecture

Fig. 4 :

 4 Fig. 4: An overview of components required to built a multichannel robotic platform used for this research

Fig. 5 :

 5 Fig. 5: The multi-channel proof of concept robotic platform.

Algorithm 1 :

 1 Proposed method McLAM. Input : Trajectory W = {(w x1 , w y1 ), . . . , (w x K , w y K )}; Output: Reflector position estimates P = {(p x1 , p y1 ), . . . , (p x K , p y K )}; Initialization: P = {}, DOA = {}, TOA, = {}, Φ = [0 • ; 360 • ]; for k = 1, . . . , K do Probe the environment with s(n); Record echoes in y k ; Transform signals to frequency domain s

Fig. 6 :

 6 Fig. 6: Comparison of the proposed method against state-of-the-art

Fig. 7 :Fig. 8 :

 78 Fig. 7: Evaluation of proposed McLAM method and state-of-the-art against different SINR

Fig. 9 :

 9 Fig. 9: Acoustic image of MPDR-NLS and DSB-NLS beamformer

Fig. 10 :

 10 Fig. 10: Evaluation of proposed McLAM method and state-of-the-art against varying Distances

Fig. 11 :

 11 Fig. 11: Evaluation of the proposed method using proof of concept robotic platform against different SDNR LIDAR = 1 m SINR = 0 dB SINR = 10 dB SINR = 20 dB SINR = 30 dB SINR = 40 dB Methodsµ [m] RMS error [m] µ [m] RMS error [m] µ [m] RMS error [m] µ [m] RMS error [m] µ [m] RMS error [m] MPDR-NLS1.0558 0.2843 0.9797 0.0992 0.9718 0.0281 0.9890 0.0118 0.9861 0.0138 DSB-NLS 1.0231 0.2802 0.8647 0.8896 0.1103 0.1174 0.9075 0.0924 0.9861 0.0138 EM-MC 1.0229 0.2803 0.8647 0.1485 0.8908 0.1094 0.9075 0.0924 1.0040 0.0039 LCMV-NLS γ = 0.1 0.9899 0.2603 0.8758 0.1647 1.0387 0.0556 1.0647 0.0647 0.9861 0.0138 LCMV-NLS γ = 1 1.0325 0.2819 0.8813 0.1409 0.7996 0.2134 0.8084 0.2042 1.0254 0.0254

  (a) Layout of office with glass surfaces.

  2D map of an office with glass surfaces.

Fig. 12 :

 12 Fig. 12: Detecting of glass surface at Aalborg University.

  2D map of the Sound Lab without GLRT.

  2D map of the Sound Lab with GLRT.

Fig. 13 :

 13 Fig. 13: Generating a spatial map of the Sound Lab.

Algorithm 2 :

 2 Proposed method ScLAM. input : Trajectory W = {(w x1 , w y1 ), . . . , (w x K , w y K )}, Initialization P = {}, TOA, = {}; output : Reflector position estimates P = {(p x1 , p y1 ), . . . , (p x K , p y K )}; for k = 1, . . . , w k do

	Acquire direction of robot movement: θ r,k ;
	Acquire direction of loudspeaker: θ l,k ;
	Probe the environment with s(n) ;
	Record echo: y k ;
	Transform signals to frequency domain
	s

TABLE I :

 I Evaluation of the proposed McLAM against ground truth and SDNRs

	SINR = 40 dB	MPDR-NLS	DSB-NLS	EM-MC	LCMV-NLS γ = 0.1

TABLE II :

 II Evaluation of the proposed McLAM against ground truth and Distances

The dataset used for simulation and evaluation can be found at http://homes.create.aau.dk//ussa/journal/index.php


John M Dudley

Goëry Genty

Frederic Dias

Bertrand Kibler

Nail Akhmediev

The dynamics of a developing CW supercontinuum: analytical predictions and experiments

Keywords:

We show that the development of the supercontinuum spectrum in the quasi-CW regime can be interpreted analytically in terms of Akhmediev Breathers. Theory and experiment are in excellent agreement.

Introduction

Despite intensive research over three decades, there are surprisingly-few analytic solutions describing nonlinear pulse evolution in optical fibers. Certainly there are exceptions such as the well-known soliton and self-similar solutions [START_REF] Akhmediev | Nonlinear Fiber Optics[END_REF][START_REF] Dudley | Self Similarity in Ultrafast Nonlinear Optics[END_REF], but complex processes such as supercontinuum generation have generally been considered to require numerical approaches [START_REF] Dudley | Supercontinuum Generation in Photonic Crystal Fiber[END_REF]. In this paper, however, we report on the use of the "Akhmediev Breather" to describe the development of the supercontinuum spectrum when seeded from spontaneous noise.

Akhmediev Breather dynamics and the onset of supercontinuum generation

The Akhmediev Breather (AB) is an exact analytic solution describing the evolution with distance of a modulated continuous wave with distance in an optical fiber [START_REF] Akhmediev | Modulation instability and periodic solutions of the nonlinear Schrodinger equation[END_REF]. In fact, the AB solution describes the dynamics of induced modulation instability (MI), and it is surprising that, although these results are nearly 25 years old, they have largely been overlooked in studies of MI and its applications. The AB evolution of a CW wave with small initial modulation yields a train of ultrashort pulses that undergoes an initial phase of temporal compression followed by subsequent temporal broadening. At the point of maximum temporal compression, the individual subpulses have maximum amplitude and minimum temporal width given by: mod 0 mod

(1 4 ) 2 cos( ) (0 , ) . 2 cos( ) 1 aa T Az T P aT Z Z (1)
Here 2a = [1 -(Z mod / Z c ) 2 ] with Z c 2 J P 0 / |E 2 |with Z mod and P 0 the modulation frequency and peak power of the initial wave, and J and E 2 <0 the fiber nonlinearity and group velocity dispersion. The corresponding spectrum consists of discrete frequency sideband modes with separation Z mod , and intensities that decrease following a geometric progression [START_REF] Akhmediev | Modulation instability and periodic solutions of the nonlinear Schrodinger equation[END_REF]. For the particular case of the maximally-compressed solution of Eq. ( 1), the spectral mode intensities are: S 0 = P 0 (2 1/2 -1) 2 for the pump component, and S n = 2 P 0 (2 1/2 -1) 2|n| (where n = ±1, ±2, ±3, … ) for the sidebands. Thus the pump and sideband intensities follow the relative progression {I 0 , I 1 , I 2 , I 3 , I 4 …} = {1, 2, 0.3431, 0.0589, 0.0101…} so that there is a 3 dB increase from the pump (n=0) to the first sideband (n=±1) and then a constant decrease of 20 log 10 (2 1/2 -1) = -7.66 dB between subsequent sidebands. This geometric progression describing the decrease in sideband amplitudes yields a characteristic triangular shape in the wings of the spectrum when plotted semi-logarithmically.

To discuss how this result can be used to interpret the onset dynamics of quasi-CW supercontinuum generation, we begin by showing experimental results: Fig. 1(a) shows spectral measurements at the output of 3.9 m of a highly nonlinear photonic crystal fiber with zero-dispersion wavelength at 780 nm for various peak powers using 1 ns pulses at 1064 nm. Experiment (top) is compared with generalised nonlinear Schrödinger equation (GNLSE) simulations (bottom). The agreement between experiment and simulation is very good and illustrate how the spectral characteristics vary significantly with power. Firstly, at the lowest power of 26 W, the spectrum consists of small number of distinct sidebands separated by the frequency of peak MI gain. As the power increases to 43 W, we enter an "extended MI" regime where we can still resolve sidebands close to the pump, but we also see the development of continuous low amplitude wings that appear triangular when plotted with a semi-logarithmic y-axis.

The spectral width in this regime approaches 60 nm at the -40 dB level, and we note the first sign of a very low amplitude Raman peak around 1120 nm. Finally at 98 W, we see a significant increase in bandwidth to ~400 nm at the -40 dB level, and the spectrum assumes characteristics of a "fully-developed" supercontinuum [START_REF] Dudley | Supercontinuum Generation in Photonic Crystal Fiber[END_REF]. We next focus on the characteristics of the extended MI spectrum at 43 W peak power. This is a regime where we are in a "NLSE dynamics" regime, before the onset of higher-order effects. To examine how the AB theory can be applied to interpret the properties of the broadened spectrum seen at 43 W, Fig. 2 The simulations use fiber parameters corresponding to experiment, and assume a 43 W peak power CW field with one photon per mode noise. The grayscale plots show results of a single numerical realization illustrating the field evolution over 3.9 m. The spectral bandwidth first reaches its maximum extension at a distance of 3.5 m where we extract the corresponding temporal profiles. The extracted temporal profile shows a high contrast pulse train, but the broadband noise that seeds the MI results in significant variation in the structure of the individual temporal peaks across the extracted profile. Nonetheless, the highest amplitude peaks in the temporal profile are fitted very well by the maximally-compressed AB sub-pulses calculated from Eq. ( 1) using Z mod corresponding to the peak of the MI gain (f mod = 1.32 THz for our parameters.) This agreement is shown explicitly in the expanded view of the shaded section in Fig. 2(b). In fact, analysis of the highest-amplitude peaks at other regions of the profile (not shown) reveals similarly good agreement with the calculated AB solution, and we have confirmed that simulations using different noise seeds reveal the same characteristics.

The results in Fig. 2 suggest that the extended MI spectrum contains a dominant contribution from spontaneously-generated AB sub-pulses generated at the peak of the MI gain. In fact, this is also confirmed using experimental spectral measurements. Specifically, we compare the form of the measured spectrum with the analytic geometric progression predicted for the maximally-compressed AB as shown in Fig. 3. The figure compares numerical and analytic results with the experimental results of Fig. 1 for P 0 = 43 W. The numerical results show averaged multi-shot simulations using the full GNLSE with a 1 ns pulse input field (blue, short dashes) and using only the NLSE with a CW input field (red, long dashes). It is clear that both GNLSE and NLSE simulations produce similar results in agreement with experiment, with the GNLSE of course yielding better agreement in predicting the onset of a small Raman peak on the long wavelength edge around 1120 nm.

The most significant result, however, is the agreement between experiment and the calculated amplitudes of the discrete frequency modes associated with the analytic form of the maximally-compressed AB. Because spontaneous MI is seeded by broadband noise, the spectral structure is not discrete, but rather consists of a continuous span of frequencies. Nonetheless, as with the time domain characteristics, we expect the spectral characteristics to contain a dominant component from the maximally-compressed AB solution calculated at peak MI gain. Specifically, normalising to the n = ± 1 sideband amplitudes seen in experiment, the decay of spectral intensity with frequency in experiment and simulations is reproduced very well by the analytic geometric progression of the maximally-compressed AB breather. Fig. 3 Comparison between experiments (solid black line), numerical simulations using the full GNLSE (blue dashed line), numerical simulations using the NLSE only (red dashed line), and the calculated spectrum of the maximallycompressed AB (green lines from zero).

Conclusions

These studies of the temporal and spectral properties of spontaneous MI have shown that the characteristics of both the modulated temporal profile and the associated spectrum can be explained in terms of the development of high amplitude Akhmediev Breather sub-pulses. The analytic form of the AB spectrum has been shown to be in very good quantitative agreement with the wings of the spontaneous MI spectrum observed in experiments studying the onset of quasi-CW supercontinuum generation. These results represent an important step in applying the analytic formalism of breather evolution to the interpretation of supercontinuum generation. Moreover, in the case of induced MI from dual frequency fields, the AB theory would be expected to also be applied to fields such as optimization of pulse train and frequency comb generation from multiwave mixing processes [START_REF] Barviau | Generation from continuous waves of frequency combs with large overall bandwidth and tunable central wavelength[END_REF].

Fig. 1 .

 1 Fig. 1. Experimental (top) and simulation (bottom) results for 1 ns pulses at 1064 nm injected into highly nonlinear PCF at peak powers as shown. Simulation results are averaged and convolved with a resolution function matching the bandwidth of the spectrum analyzer used.

  presents simulation results showing in more detail how the (a) spectral and (b) temporal properties develop with distance.

Fig. 2 .

 2 Fig. 2. Single shot (a) spectral and (b) temporal evolution. The shaded region in the extracted temporal trace is shown in detail in the rightmost figure, comparing simulations (solid line) with the AB solution at a modulation frequency corresponding to peak MI gain (dashed line).


J M Dudley
email: john.dudley@univ-fcomte.fr

B Kibler

G Genty

S Coen

P Kinsler

From Supercontinuum Generation to Carrier Shocks: Extreme Nonlinear Propagation in Photonic Crystal Fiber

We review supercontinuum generation in photonic crystal fiber and discuss the underlying spectral broadening processes from the femtosecond to the continuous wave regime. We also describe a new propagation model that integrates carrier dynamics.

Introduction

Since its first observation by Ranka et al. in 2000, supercontinuum (SC) generation in photonic crystal fiber (PCF) has been intensively studied over a wide range of experimental conditions [START_REF] Ranka | Visible Continuum Generation in Air-Silica Microstructure Optical Fibers with Anomalous Dispersion at 800 nm[END_REF][START_REF] Dudley | Supercontinuum Generation in Photonic Crystal Fiber[END_REF]. Nonlinear spectral broadening from self-phase modulation (SPM) plays a central role in SC generation, and its interaction with the novel fiber dispersion results in rich dynamics that are manifested in different ways depending on the particular parameter regime considered. Despite the complexity of these dynamics, numerical modeling has now resulted in a comprehensive physical understanding of the SC broadening mechanisms, and quantitative agreement between simulation and experiment has been obtained even for spectral bandwidths around an octave [START_REF] Corwin | Fundamental Noise Limitations to Supercontinuum Generation in Microstructure Fiber[END_REF].

In this paper, we will present a synthetic review of the physics of SC generation in PCF as it is currently understood. In addition, we will discuss recent theoretical developments considering interaction regimes where the χ (3) nonlinearity must be considered in a way that goes beyond a simple description involving only SPM. A particular result is the presentation of a new nonlinear propagation equation that extends the powerful concept of the pulse envelope to model dynamical effects on the underlying ultrafast carrier, and we use this model to examine the consequences for nonlinear propagation effects in PCF. The key features of our approach are a generalized nonlinear polarization term including third harmonic generation (THG), and the integration of carrier phase evolution effects onto the complex pulse envelope itself. Our extended envelope simulations are compared with the numerical integration of Maxwell's equations, and excellent quantitative agreement is obtained.

Generalized Nonlinear Envelope Equation

Our development follows the Blow-Wood time-domain formalism [START_REF] Blow | Theoretical description of transient stimulated Raman scattering in optical fibers[END_REF], but includes the full field χ (3) nonlinearity, and uses a rigorous approach to forward and backward envelope factorization. Subject only to two approximations of (i) unidirectional propagation and (ii) a perturbative nonlinearity, we obtain a first-order propagation equation with no constraint on the temporal structure of the propagating envelope. This allows us to lift "slowly varying envelope" and "slowly evolving wave" approximations that have previously limited the application of envelope-approaches to sub-cycle regimes [START_REF] Brabec | Nonlinear optical pulse propagation in the single-cycle regime[END_REF]. Specifically, writing the field E(z,t) = Re[A(z,t) exp(-iω 0 t)], the propagation equation for A(z,t) in the laboratory frame can be expressed as:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + = ∂ ∂ - ∂ ∂ - ≥ + ∑ ) , , ( 3
1 ) 1 ( 1 ! 2 2 2 0 1 0 A t z g f A A e A f t i i t A k i z A R t i R shock k k k k k ω τ γ β (1)
where, on the right hand side, the term proportional to |A| 2 A describes SPM, the term proportional to A 3 describes THG, f R represents the fractional Raman contribution to the nonlinearity and g(z,t,A) is a generalized Raman nonlinear response function [START_REF] Karasawa | Comparison between theory and experiment of nonlinear propagation for a-few-cycle and ultrabroadband optical pulses in a fused silica fiber[END_REF]. This generalized envelope propagation equation allows for the accurate simulation of extreme nonlinear propagation phenomena whose study has previously been accessible only through the numerical solution of Maxwell's equations.

To illustrate the success of Eq. ( 1) in modeling ultra-broadband nonlinear dynamics, we first consider the propagation of 5 fs (FWHM) sech 2 pulses of peak power 100 kW in an ideal dispersionless χ (3) medium with nonlinearity γ = 0.68 W -1 m -1 . The input wavelength is 830 nm and we consider propagation over 3.71 μm. Fig. 1(a) shows the output temporal characteristics obtained using numerical simulations based on Eq. ( 1), with the exploded view clearly illustrating the onset of a shock on the optical carrier. This shock arises from the presence of multiple higher-order harmonics, and these are apparent in the field spectrum that is also shown. Significantly, we have compared the temporal and spectral characteristics obtained from our generalized envelope Eq. ( 1) with simulations that directly integrate Maxwell's equations using the pseudospectral spatial domain (PSSD) method [START_REF] Tyrrell | Pseudospectral spatial-domain: A new method for nonlinear pulse propagation in the fewcycle regime with arbitrary dispersion[END_REF]. These results are shown as the open circles in Fig. 1(a). We stress the remarkable quantitative agreement that is obtained, even for this extreme case where spectral broadening extends over 10 times the carrier frequency, and the carrier exhibits shock dynamics on a sub-10 attosecond timescale. Additional simulations show similar agreement with increased propagation, and in the presence of dispersion.

These results confirm the validity of our generalized envelope model to describe ultra-broadband and sub-cycle propagation. We have thus used numerical simulations based on Eq. ( 1) to consider a fully realistic case in a 1 mm long silica tapered PCF nanowire of diameter 600 nm. With parameters as above, but also including dispersion, Raman scattering, and the frequency dependence of the guided mode area, the results are shown in Fig. 1(b). The SC in this case spans over more than 2 octaves, significantly exceeding the range that has been successfully modeled in previous envelope-based numerical studies.

Fig. 1 (a) Shock formation after 3.71 μm propagation in a dispersionless χ (3) medium. The left subfigure shows the temporal field obtained from the solution from the generalized nonlinear envelope equation, together with an exploded view comparing these envelope simulation results (solid line) with those obtained from Maxwell's equations (circles). The right subfigure compares the corresponding spectrum obtained using the generalized nonlinear envelope equation (solid line) with that from Maxwell's equations (circles). (b) Temporal and spectral characteristics obtained using Eq. ( 1) to model SC generation in a 1mm nanowire of 600 nm diameter including dispersion, Raman scattering and modal effective area variation.

Conclusions

These results extend the use of an envelope description of pulse propagation into the sub-cycle regime, allowing the convenient modeling of carrier shocks and multi-octave SC generation. Significantly, additional simulations confirm the validity of this approach at all intensity levels less than the damage threshold of fused silica, and thus these results are particularly relevant to studying high intensity pulse propagation in PCF nanowires. We also suggest that these results could motivate a new generation of experiments studying nonlinear carrier shaping effects.