Mohamed Lotfi Derouiche
email: mohamedlotfi.derouiche@enit.rnu.tn

Soufiene Bouallègue
email: soufiene.bouallegue@issig.rnu.tn

Joseph Haggège
email: joseph.haggege@enit.rnu.tn

Guillaume Sandou
email: guillaume.sandou@centralesupelec.fr

LabVIEW Perturbed Particle Swarm Optimization Based Approach for Model Predictive Control Tuning

Keywords: Model Predictive Control, parameters tuning problem, perturbed Particle Swarm Optimization, Genetic Algorithm, LabVIEW implementation, MAGLEV system

In this paper, a new Model Predictive Controller (MPC) parameters tuning strategy is proposed using a LabVIEW-based perturbed Particle Swarm Algorithm (pPSA). This original LabVIEW implementation of this metaheuristic algorithm is firstly validated on some test functions in order to show its efficiency and validity. The optimization results are compared with the standard PSO approach. The parameters tuning problem, i.e. the weighting factors on the output error and input increments of the MPC algorithm, is then formulated and systematically solved, using the proposed LabVIEW pPSA algorithm. The case of a Magnetic Levitation (MAGLEV) system is investigated to illustrate the robustness and superiority of the proposed pPSA-based tuning MPC approach. All obtained simulation results, as well as the statistical analysis tests for the formulated control problem with and without constraints, are discussed and compared with the Genetic Algorithm Optimization (GAO)-based technique in order to improve the effectiveness of the proposed pPSA-based MPC tuning methodology.

INTRODUCTION

Being part of robust and optimal control techniques, the Model Predictive Control (MPC) approach is almost pervasive in various engineering applications [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]; [START_REF] Rawlings | Model Predictive Control: Theory and Design[END_REF]; Bao-Cang (2010); [START_REF] Rossiter | Model Based Predictive Control: A practical Approach[END_REF]. The first generation of predictive control algorithms tried to solve constrained control problems manipulating several variables in the petrochemical field in the late 80s. MPC applications are present in many systems where the PID controller reaches its limits. It covers both the Single-Input Single-Output (SISO) and Multi-Inputs Multi-Outputs (MIMO) problems for linear and nonlinear systems.

However, for model predictive controller design, the relationship between the adjustment parameters and physical phenomena is shady which make tuning task tedious, error-prone and out of reach of non-expert users. Therefore, this type of tuning still presents a difficulty for manufacturers who do not have particular expertise. Well tuning MPC parameters is still challenging in control engineering field. This challenge leads the MPC tuning to be the subject of several research projects. In [START_REF] Lino | A new input-output constrained model-predictive control with frequency-domain tuning technique and its application to an ethylene plant[END_REF], authors proposed a tuning method taking into account the robust stability based on the frequency response analysis. [START_REF] Rowe | Tuning mpc using h-infinity loop shaping[END_REF] proposed a tuning method using the H ∞ loop-shaping concept. A heuristic method based on setting rules established by experts has been proposed in [START_REF] Drogies | Predictive control: propositions for the design methodology[END_REF]. The methods of C. Rowe and Y. Lino require complicated calculation procedures and keep the need for qualified or rather expert user. In practice, the method presented by S. Drogies must be modified for each specific problem. Thus, none of the above mentioned methods has led to a systematic tuning procedure allowing a user to well design an MPC controller without resort to an expert as it is the case for PID controllers.

The MPC parameters tuning problem can be formulated as an optimization program. Many approaches are proposed in the literature for dealing with hard choice of MPC parameters. See [START_REF] Sandou | Nonlinear Model Predictive Control, chapter Particle Swarm Optimization Based NMPC: An Application to District Heating Networks[END_REF]; [START_REF] Sandou | Metaheuristic strategy for the hierarchical predictive control of large scale energy networks[END_REF]. Knowing that the predictive control is treated as a black box, classic gradient-based optimization techniques will not be able to solve this problem. However, metaheuristics algorithms [START_REF] Yang | Engineering optimization: an introduction with metaheuristic applications[END_REF]; [START_REF] Dréo | Metaheuristics for Hard Optimization: Methods and Case Studies[END_REF]; [START_REF] Siarry | Advances in Metaheuristics for Hard Optimization[END_REF]; [START_REF] Gendreau | Handbook of Metaheuristics[END_REF], which have reached a remarkable maturity, especially the Particle Swarm Optimization technique, present a solution to such problems. Recently, [START_REF] Suzuki | Model predictive control using PSO[END_REF] introduced an automatic tuning of predictive controller parameters using the standard PSO technique. Newly, [START_REF] Xinchao | A perturbed particle swarm algorithm for numerical optimization[END_REF] introduced a new "perturbed PSO" algorithm, denoted as pPSA. This metaheuristic comes to refine the standard PSO technique, which after a premature convergence, is often dropped into local minima. Based on the principle of the possibility measure, the pPSA pursues to update the optimum position of the particle swarm noted gbest using the "possibility at gbest" strategy that seeks to model the lack of information regarding the true optimality of gbest. The pPSA is proposed so as to escape from the local optimal trap and the rapid loss of diversity during the optimization process. The new particle updating strategy aims to promote exploration capabilities in the search space. In this paper, a LabVIEW implementation of this efficient metaheuristic algorithm is proposed for using in MPC parameters tuning problem. An application to a position control of the suspended sphere in known Feedback MAGLEV system is investigated in order to show the effectiveness of our proposed MPC tuning approach.

The remainder of this paper is organized as follows. In Section 2, the MPC parameters tuning problem is introduced and formulated as an optimization program to be solved later using the proposed metaheuristic-based technique. A brief overview on MPC approach is firstly described in this paragraph. Section 3 presents the proposed perturbed PSO-based algorithm as well as its original implementation under graphical LabVIEW tool. This algorithm is validated through an optimization benchmark functions from the literature. An application to the MPC tuning parameters for a MAGLEV system using GA and the proposed pPSA is investigated in Section 4. All simulation results, obtained for the unconstrained and constrained formulations of such a problem, are presented and discussed. Section 5 concludes this paper.

MODEL PREDICTIVE CONTROL PARAMETERS TUNING PROBLEM

Basic Concepts of MPC Design

The MPC design approach focuses on constructing and optimizing feedback controllers at each discrete-time instant, see [START_REF] Rawlings | Model Predictive Control: Theory and Design[END_REF]; Bao-Cang (2010); [START_REF] Rossiter | Model Based Predictive Control: A practical Approach[END_REF], as shown in Fig. 1. Such predictive controllers can adjust the control action before a change in the output setpoint actually occurs. The ability of handling operational constraints in an explicit manner is one of the main reasons for the popularity and success of MPC strategies in various industrial applications, as shown in [START_REF] Sandou | Nonlinear Model Predictive Control, chapter Particle Swarm Optimization Based NMPC: An Application to District Heating Networks[END_REF]; [START_REF] Sandou | Metaheuristic strategy for the hierarchical predictive control of large scale energy networks[END_REF].

As depicted in Fig. 1, the main elements of the discretetime model predictive control are the plant input, the controlled output and the reference trajectory which are denoted by u ∈ R, y ∈ R, and r ∈ R, respectively. The plant model determines the predicted plant outputs on the prediction horizon, denoted N p . The optimization algorithm is aimed at determining the control sequence given by {u (k -1 + i) , i = 1, 2, . . . , N c } for the control horizon, denoted N c . Only the first element u * (k) of the optimized control sequence is applied to the plant and the control input is updated at each sampling instant. The optimization process is repeated at the next sampling time, on the basis

(k -1 + i) = u (k + N c -1) for N c < i ≤ N p .
A model describing the behavior of the controlled plant is required in MPC framework. It is assumed that such a plant is represented by a discrete-time model of the form:

x (k + 1) = A A Ax (k) + B B Bu (k) y (k) = C C Cx (k) (1)
This model is used to compute system predictions over a finite prediction horizon of N p samples. So, at every sampling time and for a specified prediction and control horizons, the MPC controller attempts to minimize the following cost function:

J (k) = Np i=1 ê (k + i|k) T Q Q Qê (k + i|k) + Nc-1 i=0 ∆u T (k + i|k) R R R∆u (k + i|k) (2)
where ê (k + i|k) = ŷ (k + i|k) -r (k + i|k), ŷ (k + i|k), r (k + i|k) and ∆u (k + i|k) are the predicted plant output, the output setpoint profile and the predicted increment of change in control action at time k + i, given all measurements up to and including those at sampling-time k, respectively.

Q Q Q = Q Q Q T > 0 and R R R = R R R T > 0 are
weighting matrices on the output error and input increments.

Minimizing the cost function ( 2) is usually subject to the operational constraints on the control action, its rate of change and plant output signals. Thus, limitations of these entities can be defined, respectively, as given by Eqs. ( 3), ( 4) and ( 5):

u min ≤ u (k) ≤ u max (3) ∆u min ≤ ∆u (k) ≤ ∆u max (4) y min ≤ y (k) ≤ y max (5)

MPC Tuning Parameters Problem Formulation

A formulation of an optimization based tuning problem of parameters in MPC paradigm is given in the followings. Since prediction and control horizons, as well as operational constraints of MPC problem are fixed, weightings matrices Q Q Q and R R R of Eq. ( 2) are tuned in this study. The concept of the proposed pPSA based MPC tuning procedure is depicted in Fig. 2.

In this paper, we consider the control of SISO systems so that Q Q Q and R R R weighting matrices become scalars to be tuned within a time-domain framework. So, the cost function, which penalizes tracking errors and to be minimized, is given as follows:

J 1 (τ, Q Q Q, R R R) = +∞ τ =0 [y (τ ) -r (τ )] 2 (6)
where τ is the time expressed as number of samples from the beginning of simulation until achieving the steadystate domain.

PROPOSED PERTURBED PSO-BASED APPROACH FOR MPC TUNING

Perturbed PSO Algorithm

The perturbed particle swarm optimizer pPSA maintains a swarm of n P ART particles in the search space with dimension D. Each particle in the swarm is characterized by its current position x x x i t = x i,1 t , x i,2 t , . . . , x i,D t , the previous best position p p p i t = p i,1 t , p i,2 t , . . . , p i,D t , and

the velocity v v v i t = v i,1 t , v i,2 t , . . . , v i,D t at iteration t ∈ [[1, n GEN ]].
Where n GEN is the number of generations.

The trajectory of each particle is updated according to its own flying experience and that of the best particle in the swarm as given by the followings motion equations:

x x x i t+1 = x x x i t + v v v i t+1 (7) v v v i t+1 = wv v v i t + c 1 r i 1,t p p p i t -x i t + c 2 r i 2,t p p pg t -x x x i t ( 8
) where w represents the inertia factor, c 1 and c 2 are the cognitive and the social scaling factors respectively, r i 1,t and r i 2,t are random numbers uniformly distributed in the interval [[0, 1]], and p p p i t is the best previously obtained position of the i th particle.

In contrast to the standard PSO algorithm, the gbest p p pg t in Equation ( 8) of pPSA, denoted as p-gbest for "possibly at gbest", is characterized by a normal distribution as follows:

p p pg t = N (p p p g t , σ)

where σ represents the degree of uncertainty about the optimality of the gbest p p p g t of the standard PSO algorithm as shown in [START_REF] Xinchao | A perturbed particle swarm algorithm for numerical optimization[END_REF]. So, for perturbed particle updating strategy, σ is modeled as some non-increasing function of the number of cycles t, called max-min model and given by Xinchao (2010) as follows:

σ (t) = σ max , t < α × n GEN σ min , otherwise (10)
where σ max , σ min and α are manually fixed parameters.

Two others non-increasing models (linear and random) are proposed in [START_REF] Xinchao | A perturbed particle swarm algorithm for numerical optimization[END_REF] and can be used instead Eq. ( 10). On the other hand, the function of p-gbest is to encourage the particles to explore a region beyond that defined by the search trajectory. It provides a simple and efficient exploration at the early stage when α is large and encourages local fine-tuning at the latter stage when α is small.

Finally, a pseudo code of this swarm algorithm is given, for a minimization problem, as follows:

•

Step 1: Initialize a population of n P ART particles having random positions and velocities on D dimensions of the search space; • Step 2: At every iteration t and for each particle x x x i t , evaluate the considered optimization fitness function on the D decision variables;

• Step 3: Compare particle's fitness evaluation J i 1,t = J 1 x x x i t with its pbest i t = J 1 p p p i t . If J i 1,t ≤ pbest i t then pbest i t = J i 1,t and p p p i t = x x x i t ; • Step 4: Identify the particle in the neighborhood with the best success so far and assign its position to the global best variable noted gbest; • Step 5: Change the velocity and position of the particle according to the motion equations ( 7) and ( 8); • Step 6: If the termination criterion is met (satisfied fitness or max of iterations), the algorithm terminates with the solution x x x * = arg min

x x x i t J 1 x x x i t , ∀i, t . Otherwise, go to Step 2.

Implementation of the pPSO Algorithm

In order to validate the LabVIEW implementation of perturbed PSO algorithm, a benchmark optimization with 8 functions are adopted (see Appendix A.). In order to have a fait comparison between PSO and pPSA in Xinchao (2010) and our LabVIEW coded pPSA, the parameters settings are as follows: dimension of search space D = 30, number of particles n P ART = 30, number of generations n GEN = 2000, cognitive and social coefficients c 1 = 0.5 and c 2 = 0.3, inertia weight w = 0.9 and σ max = 0.15, σ min = 0.001, α = 0.5.

All the statistical results over 30 independent runs of the implemented algorithms, as well as those published in [START_REF] Xinchao | A perturbed particle swarm algorithm for numerical optimization[END_REF], are summarized in Tables 1 and2. From these numerical optimization results, it is clear that the pPSA greatly outperforms the standard PSO algorithm. The good performance, in terms of better exploring ability for promising area and exploiting for locating the optima, is carried out.

APPLICATION TO THE MAGLEV SYSTEM CONTROL

Plant Description and Modeling

The Feedback Magnetic Levitation system (MAGLEV 33-006) of Fig. 3 is used as a process example to validate the pPSA-based approach for MPC tuning parameters.

In required range of operation, the distance h of the suspended sphere is given by the infrared photo-sensor voltage y as follows : y = γh + y 0 (11) where γ is a positive gain depending on the position sensor, and y 0 is the offset voltage such that y ∈ [-2V, +2V ].

The coil current is regulated by an inner control loop within the driver block. Its characteristic is linearly related to the input voltage u as follows, neglecting its high frequency dynamics: i = ρu + i 0 (12) where ρ > 0 is the coil resistor and i 0 > 0 is the offset value of current. The working excursion of u is limited between -3V, corresponding to a null coil current, and +5V that defines the saturation value as given in Feedback; [START_REF] Santos | Robust model predictive control for a magnetic levitation system employing linear matrix inequalities[END_REF]. For this process, the predictive control of the suspended sphere position is investigated. The dynamics of the vertical movement of such a sphere is modeled as follows:

m d 2 h dt 2 = mg a -K i 2 h 2 (13)
where K is an electromechanical conversion gain depending on the MAGLEV system, m is the mass of the sphere, g a is the acceleration of gravity, and i is the coil current.

According to the given sensor and current driver characteristics ( 11) and ( 12), the equation ( 13) can be re-written as:

m d 2 y dt 2 = γmg a -K (ρu + i 0 ) 2 γ 3 (y -y 0 ) 2 (14) Taking x = [ y ẏ ]
T as state vector, the following state space representation of the studied system is obtained:

   ẋ1 = x 2 ẋ2 = γg a -K (ρu + i 0 ) 2 γ 3 m (x 1 -y 0 ) 2 (15)
Finally, while using the values for the physical model parameters, a linear discrete-time model of the MAGLEV plant, where a zero-order-hold with 5 ms sampling period was adopted at the input of the system as in [START_REF] Santos | Robust model predictive control for a magnetic levitation system employing linear matrix inequalities[END_REF], is given by the following state-space representation:

   x k+1 = 1.0108 0.0050 4.3185 1.0108 x k + -0.0142 -5.6779 u k y k = [ 1 0 ] x k (16)

Simulation Results and Discussion

As described in section 2.2, the optimization is performed on a SISO system using the unconstrained control problem (6). For this design, we use the control parameters of Table 3. The simulation results for the pPSA and GA algorithms are summurized in Table 4. Recall that the GA algorithm is also implemented and executed under LabVIEW environment, see W. [START_REF] Golebiowski | Waptia: Genetic optimization algorithm[END_REF] and National Instruments (2009).

Table 3. GA and pPSA control parameters. The controlled position of the suspended sphere as well as the control input voltage of the MAGLEV current driver are shown in Fig. 4 and Fig. 5, respectively. The tracking performances of the controlled system are satisfactory but a high amplitude of the control action is noted and exceeds the saturation value, see Fig. 5. In an experimental point of view, this control voltage amplitude can not be applied to the real MAGLEV system which have an upper limit of +5V. This then requires the reformulation of the control problem to incorporate this kind of operational constraints in the MPC synthesis procedure.

In order to control the real MAGLEV system and to show the effectiveness of the proposed pPSA based MPC approach, operational constraint on upper limit of the control action signal u max = 1V is considered. The problem can be formulated using as a constrained optimization problem:

       J 2 (τ, Q Q Q, R R R) = +∞ τ =0 [y (τ ) -r (τ )] 2 s.t : max τ ≥0 |u (τ )| ≤ u max (17)
PSO being an unconstrained optimization algorithm, the objective function, to be minimized using the LabVIEW coded pPSO algorithm, is augmented with the external static penalty technique, that is proportional to the degree of constraint infeasibility, and becomes as follows:

J 2 (τ, Q Q Q, R R R) = +∞ τ =0 [y (τ ) -r (τ )] 2 +exp   1000 max τ ≥0 |u (τ )| -u max u max   (18)
The LabVIEW implementation for the constrained optimization problem (18) leads to the simulation results of Fig. 6 and Fig. 7 as well as the results in Table 5. For this design, both metaheuristic algorithms improve good performances of tuned MPC controller in terms of stability, tracking dynamic behavior and handling operational constraints. However, the superiority of our proposed LabVIEW pPSA based approach is shown mainly in terms of the fastness convergence and the simplicity of implementation.

It is important to highlight that in the MPC parameters Q Q Q and R R R optimization, only the ratio of these two coefficients is most important. The same setting is obtained if we divide by 2 these two coefficients. This is shown in the comparison of the LabVIEW implemented pPSA and GA, i.e., 0.898/0.226 = 3.97 and 59.36/14.95 = 3.97. As depicted in Table 5, this is the same tuning for the two proposed methods. The main observation is that these are the same two settings found by two different metaheuristic algorithms, so there is a greater probability that the obtained solution is the global optimum of the MPC tuning optimization problem. The superiority of the proposed pPSA technique is shown mainly at the remarkable fastness convergence in solving the MPC tuning problem as well as its simple software implementation and algorithm parameters choice.

CONCLUSION

In this paper, a new perturbed PSO based approach for MPC parameters tuning have been proposed and successfully applied to the position control of a MAGLEV system. The proposed pPSA algorithm is implemented under LabVIEW graphical environment and validated through a benchmark of test functions from the literature. All optimization results are compared with the standard PSO and statistical analysis are carried out in order to show the validity of such an implementation. The MPC tuning problem, such as the choice of weighting factors on the output error and input increments signals, is after that formulated as an optimization problem, with and without operational constraints, and solved using the GA (WAPTIA) and the proposed LabVIEW-based pPSA algorithm for the considered MAGLEV plant.

Forthcoming works deal with final implementation of the control laws on the real MAGLEV benchmark available in our laboratory. Increasing the size of the optimizationbased tuning problem, i.e. while considering full weighting matrices, is also investigated. x 2 i -10 cos (2πx i ) + 10 [-5.12, 5.12] n 0

f 8 = 1 4000 n i=1 x 2 i - n i=1 cos x i √ i
+ 1 [-600, 600] n 0

Fig. 1 .

 1 Fig. 1. Model Predictive control structure. of the measured (or estimated from the available output y (k)) state x (k). This sequence minimizes the specified cost function (2) in the MPC design formalism, subject to problem constraints on the input, change in input and output of the plant. The optimization algorithm assumes that u (k -1 + i) = u (k + N c -1) for N c < i ≤ N p .

Fig. 2 .

 2 Fig. 2. Proposed pPSA based MPC tuning approach.

Fig. 3 .

 3 Fig. 3. Feedback Magnetic Levitation system.

  Fig. 4. pPSA-based control of the MAGLEV suspended sphere position: unconstrained optimization case.

  Fig. 6. pPSA-based control of the MAGLEV suspended sphere position: constrained optimization case.

Table 1 .

 1 Numerical optimization results over 30 runs of standard PSO from[START_REF] Xinchao | A perturbed particle swarm algorithm for numerical optimization[END_REF]. Where Best, Mean, Median and STD mean the best, average, median result and the standard deviation over 30 runs.

	Function Best	Mean	Median	STD
	f1	4.39	42.37	27.02	49.76
	f2	2.03	19.95	15.37	18.04
	f3	1.23e+03 6.77e+03 4.88e+03 5.37e+03
	f4	1.47e+02 2.81e+03 9.98e+02 5.23e+03
	f5	1.67e+02 5.62e+02 5.04e+02 2.94e+02
	f6	3.81e-01	8.20e-01	7.2e-01	3.79e-01
	f7	39.50	84.85	82.76	22.71
	f8	1.03	1.55	1.27	7.58e-01

Table 2 .

 2 Numerical optimization results over 30 runs of pPSA under LabVIEW.

	Function Best	Mean	Median	STD
	f1	5.72e-06	0.00	7.32e-06 8.15e-07
	f2	1.19e-02	2.4e-01 5.62e-02 5.56e-01
	f3	1.66e-03	1.90	6.55e-01 2.79
	f4	1.73e+01 62.40	28.10	57.10
	f5	0.00	24.10	12.00	56.10
	f6	1.93e-11	0.00	6.42e-11 1.88e-11
	f7	53.72	81.22	82.58	15.56
	f8	3.24e-07	0.00	7.40e-03 2.33e-02

Table 4 .

 4 Performances comparison of GA and pPSA for Problem (6).

		pPSA	GA (WAPTIA)
	weight factor Q Q Q	36.8780	0.6050
	weight factor R R R	0.0004	7.00e-06
	gbest value	0.1889	0.1887
	umax	26.3265	25.9861
	computation time 12min, 22sec 2h, 29min, 52sec

Table 5 .

 5 Performances comparison of GA and pPSA for Problem (18).

		pPSA	GA (WAPTIA)
	weight factor Q Q Q	14.9501	0.226115
	weight factor R R R	59.3676	0.898016
	gbest value	2.85314	2.85314
	umax	0.995273	0.995252
	computation time 18min, 47sec 2h, 21min, 23sec


B Naylor

M Brewczyk

M Gajda

O Gorceix

E Maréchal

L Vernac

B Laburthe-Tolra

Competition between Bose Einstein Condensation and spin dynamics

Keywords: numbers: 03.75.Mn, 05.30.Jp, 67.85.-d, 05.70.Ln

We study the impact of spin-exchange collisions on the dynamics of Bose-Einstein condensation, by rapidly cooling a chromium multi-component Bose gas. Despite relatively strong spin-dependent interactions, the critical temperature for Bose-Einstein condensation is reached before the spindegrees of freedom fully thermalize. The increase in density due to Bose-Einstein condensation then triggers spin dynamics, hampering the formation of condensates in spin excited states. Small metastable spinor condensates are nevertheless produced, and manifest strong spin fluctuations.

Dilute quantum gases are especially suited for the investigation of non-equilibrium dynamics in closed or open quantum systems, for example associated to the physics of thermalization [1], prethermalization [2], or localization [3]. In particular, they provide a platform to study the kinetics of Bose-Einstein condensation. Soon after the first Bose-Einstein condensates (BECs) were obtained, it was for example possible to investigate how the BEC nucleates [4,5]. More recently, experiments performing a temperature quench below the superfluid transition investigated the dynamics of spontaneous symmetry breaking [6] and revealed the production of longlived topological defects [7]. The aim of this work is to extend the dynamical studies of Bose-Einstein condensation to the case of a multi-component Bose gas, in order to establish the mechanisms to reach both superfluid and magnetic orders. While these orders are intrinsically connected due to Bose stimulation [8,9] (which contrasts with the case of Fermi fluids [START_REF] Coleman | Introduction to many body physics[END_REF]), it was predicted that strong spin-dependent interactions induce spin ordering at a finite temperature above the BEC transition [START_REF] Natu | [END_REF].

We find that the dynamics of Bose-Einstein condensation is drastically modified due to spin-changing collisions arising from relatively strong spin-dependent interactions. Thermalization of the spin degrees of freedom is influenced by the occurrence of BEC, and in turns influences which multi-component BECs can be produced. Our experiment also demonstrates the difficulty to thermalize the spin degrees of freedom, which has a strong impact on the spin distribution of the BECs, and on their lifetime. This is of particular relevance for large spin atoms, and most notably for strongly magnetic atoms such as Cr [12,13], Er [14], and Dy [15] for which dipolar relaxation strongly limits the lifetime of multicomponent gases [16].

We induce fast evaporative cooling of a multicomponent s=3 chromium thermal cloud by lowering the depth of a spin-insensitive optical dipole trap (ODT) (see Fig 1a). When the gas is only slightly depolarized, the thermal gas of the most populated, lowest energy, m s = -3 state rapidly saturates (i.e. reaches the maximal number of atoms in motional excited states allowed by Bose statistics [17]) and a BEC is produced for this spin state. Saturation is also reached for the thermal gas in the second-to-lowest energy state m s = -2. However, surprisingly, this state fails to condense and the BEC remains fully magnetized. In contrast, when the experiment is performed with an initially more depolarized thermal gas, spinor (i.e. multi-component) condensates are obtained, although they remain very small (see Fig. 1) and show strong spin fluctuations. Comparison with numerical simulations based on the classical field approximation (CFA) [18] reveals that the difficulty to obtain a multi-component BEC is due to spin exchange collisions, which rapidly empty the condensates in spin excited states by populating spin states for which the thermal gas is not yet saturated. There is an intriguing interplay between condensation and spin dynamics, as the large increase in density associated to BEC triggers fast spin dynamics which in turn tends to deplete the BEC in spin excited states. The observed spin fluctuations in the BEC are ascribed to a combined effect of phase fluctuations due to symmetry-breaking at the BEC transition, and spin dynamics.

To prepare an incoherent spin mixture of thermal gases, we start from a thermal gas of 2.10 4 52 Cr atoms, at T = 1.1 × T c = 440 ± 20nK, polarized in the Zeeman state m s = -3. We adiabatically reduce the magnetic field B so that the Zeeman energy is of the same order as the thermal kinetic energy. Depolarization of the cloud is driven by magnetization-changing collisions associated to dipole-dipole interactions [19,20]. We obtain a gas of longitudinal magnetization M = -2.50±0.25, with M ≡ s i=-s in i (n i is the relative population in Zeeman state m s = i). We then reduce the trap depth by applying an approximately linear ramp to the ODT laser intensity in a time t S . This results in fast forced evaporative cooling of all Zeeman states (which we refer to as "shock cooling", see Fig. 1a). We study spin dynamics and condensation dynamics by measuring both the spin and momentum distributions as a function of the time t after the beginning of the evaporation ramp. This measurement is performed by switching off all trapping lights, and applying an average magnetic field gradient of 3.5 G.cm -1 during a 6 ms time of flight to perform a Stern-Gerlach separation of the free-falling atoms.

Fig. 1c shows a typical absorption picture. It reveals a BEC in m s = -3, and a thermal gas in spin excited states. We extract the number of thermal and condensed atoms of each spin state through bi-modal fits accounting for Bose statistics. We plot in Fig. 2 the thermal atom numbers as well as the condensate fractions in m s = (-3, -2) as a function of time t for a shock cooling time t S = 500 ms. We found similar results for t S = 250 ms and t S = 1 s.

The gray area in Fig. 2 highlights a relatively long cooling time during which the m s = -3 and m s = -2 gases hold approximately the same number of thermal atoms, and there is a BEC in the lowest state m s = -3 but not in m s = -2. This phenomenon is surprising, because it shows that the m s = -2 component fails to undergo Bose-Einstein condensation although its thermal gas is saturated. Indeed, m s = -2 and m s = -3 thermal atoms have the same measured mechanical temperature (within our 5% experimental uncertainty), experience the same trapping potential, and both interact through the S = 6 molecular potential with the existing m s = -3 BEC. Therefore the m s = -2 cloud has the same thermodynamical properties than the m s = -3 thermal gas, and like this latter spin component, should condense for further cooling [17]. However, BEC does not occur in this state, until t ≤ 700 ms. Only for t ≥ 700ms do we distinguish a very small BEC also in m s = -2. This demonstrates that a BEC in m s = -2 hardly forms, although the thermal gas is saturated and cooling proceeds.

To interpret these observations, we stress that magnetization-changing collisions occur on a larger time scale than shock cooling dynamics and can be neglected, contrarily to [19]. Here, spin dynamics is almost entirely controlled by spin exchange interactions at constant magnetization driven by spin dependent contact interactions [9]. A key point is that for an incoherent mixture the spin dynamics rate γ k,l i,j for the spin changing collision (m s = i, m s = j) → (m s = k, m s = l) is set by the density of the cloud through γ k,l i,j = nσ k,l i,j v with n the atomic density, v the average relative atomic velocity and σ k,l i,j the relevant cross section within Born approximation [21]. This rate is extremely sensitive to the presence of a BEC (which enhances n). Therefore, the emergence of a BEC in a spin-excited state should trigger faster spin dynamics. In addition to these dissipative spin-exchange processes, BEC can also trigger coherent spin oscillations due to forward scattering, with a typical rate Γ k,l i,j = 4πh m n S a S i, j|S S|k, l where the sum is on even molecular potentials S, with associated scattering length a S . Our interpretation for the absence of a BEC in the state m s = -2 is thus that a large BEC cannot form in this state because fast spin-exchange processes (-2, -2) → (-1, -3) deplete the BEC as soon as it is produced. Thus spin dynamics and condensation dynamics are strongly intertwined.

To check this interpretation, we have performed numerical simulations using the Gross-Pitaevskii (GP) equation and the classical field approximation to describe thermal states. According to CFA, the GP equation determines the evolution of the classical field which is a complex function carrying the information on both the condensed and thermal atoms [18,22,23]. The initial classical field corresponds to 13.10 3 Cr atoms at the critical temperature of about 400 nK and with the experimental Zeeman distribution. To describe such a sample we follow the prescription given in [23]. Evaporative cooling is mimicked by adding a purely imaginary potential to the GP equation at the edge of the numerical grid.

Our simulations confirm the existence of a saturated m s = -2 gas and the absence of a condensate in this state (see Fig. 3). To evaluate the impact of spin-exchange processes on the dynamics of condensation, we have reproduced these simulations assuming a 4 = a 6 . In this case, the rates associated to spin-exchange processes (-2, -2) → (-1, -3), γ -1,-3 -2,-2 and Γ -1,-3 -2,-2 , which respectively scale as (a 6a 4 ) 2 and (a 6a 4 ), both vanish. As shown in Fig. 3, a BEC then forms in the spin excited state m s = -2. This confirms the crucial role of spindependent interactions in the dynamics of Bose-Einstein condensation.

It is interesting to face our observations with the accepted scenario for the thermodynamics of noninteracting multi-component Bose gases at fixed magnetization [24]. In this picture, a BEC polarized in the most populated state forms below a first critical temperature; all the other thermal spin states saturate simultaneously and condense below a second critical temperature [24]. In our situation, our observations indicate that the external degrees of freedom have reached an equilibrium, at an effective temperature which we find identical for all spin components. However, although the thermal clouds of the two lowest spin components are saturated, the other thermal clouds are not saturated. This is in profound contradiction with the prediction of Bose thermodynamics, and shows that the spin degrees of freedom in our experiment remain out of equilibrium.

This lack of thermal equilibrium for the spin degrees of freedom results from the fact that spin exchange processes for the thermal gas are slow in regards of con-FIG. 3: Numerical results. Evolution of the condensate fractions for different values of a4 (blue diamonds: ms = -3; red circles: ms = -2). Filled markers correspond to the experimental case: a4 = 64 aB and a6 = 102.5 aB [16] where aB is the Bohr radius. Empty markers correspond to simulations where a4 was set equal to a6 to suppress spin dynamics. A significative BEC fraction in ms = -2 is then obtained. densation dynamics. For example, the rate of the dominant spin exchange term, averaged over density,

γ 1 = n0σ -1,-3 -2,-2 v 2 √ 2
, with n 0 the peak atomic density, is typically 3 s -1 for a thermal gas at T C . A much longer timescale would therefore be necessary in order to reach spin equilibrium. This rate is slow compared to typical thermalization rates of the mechanical degrees of freedom e.g.

γ 2 = n0σ -2,-2 -2,-2 v 2 √ 2
≈ 40 s -1 . γ 2 >> γ 1 insures that the mechanical degrees of freedom thermalize faster than the spin degree of freedom, and that a small m s = -2 BEC can in principle be formed. However, once the m s = -2 BEC is formed, the rates associated to (-2, -2) → (-3, -1) collisions rise to typically γ 1,BEC ≈15 s -1 and Γ -1,-3 -2,-2 ≈ 100 s -1 (for 500 atoms in the condensate). Spin exchange collisions then deplete the m s = -2 BEC as fast as it is created and a multicomponent BEC cannot be sustained due to the lack of saturation of the m s = -1 thermal gas.

Under our experimental conditions, non-saturated spin-excited states thus act as a reservoir into which population may be dumped, thus preventing BEC but in the stretched state, the only collisionally stable one. The situation bears some analogies to the condensation of magnons [START_REF] Fang | Stamper-Kurn[END_REF]26] and polaritons, where BEC is obtained in the lowest momentum state by collisions of higher states in the lower polariton branch [27]. Like for polaritons, it is likely that spin-exchange interactions are increased by Bose-stimulation due to the pre-existing m s = -3 condensate.

To produce multi-component condensates during evaporation, we performed a second series of shock cooling experiments, with a lower gas magnetization M = -2.00±0.25 (where the uncertainty is associated to detection noise). The initial m s = -3 thermal gas is now depo- larized using a radio-frequency pulse. After decoherence of the spin components, this leads to a thermal incoherent mixture with initial fractional population in m s = -3, -2, -1 and 0 states approximately (31%,40%,21%,6%) with a relative uncertainty of 10%. When shock cooling is performed fast compared to γ 1,BEC , we observe the production of very small multi-component condensates in all three lowest energy states (as illustrated in Fig 1 and4).

Our numerical simulations show that spin-dynamics again a very profound influence on the dynamics of condensation. In practice, spin excited states with m s > 0 are not saturated. Therefore, spin dynamics tends to populate these non saturated states and empty the condensates, which thus remain small and short-lived. An important observation is that the spin distribution of the obtained multi-component BECs shows strong fluctuations (see Fig. 4a) compared to the thermal fraction. We interpret this feature in the following way. Bosecondensation of the different spin-excited states introduce a spontaneous symmetry breaking as the phase of each condensate is chosen randomly. This spontaneous symmetry breaking, already observed in [START_REF] Fang | Stamper-Kurn[END_REF], can also be interpreted as the production of fragmented BECs [28,29]. As spin-dynamics is sensitive to the relative phase between the condensates in the various spin states [9,30], we propose that the observed spin fluctuations result from a combined effect of spin dynamics and of the spontaneous symmetry breaking.

We performed numerical simulations to test this scenario. As the CFA does not provide a direct way to provide symmetry breaking at the BEC transition, we chose to apply random relative phases to the wave-functions describing the thermal atoms in different spin compo-nents before condensation. This provides an empirical way to simulate symmetry breaking. We performed a series of numerical simulations for different sets of relative phases between the Zeeman components. We then obtained small condensates with fluctuating magnetization (see Fig. 4b). Furthermore, we also observe that spin and condensation dynamics are also significantly modified by the applied random phases. Due to large computational time for each run, a systematic study of BEC magnetization as a function of initial phases has not yet been performed and remains to be thoroughly investigated. However, while the magnetization fluctuations obtained in the numerical simulations are typically five time smaller than the experimental measurements, these preliminary results thus support the scenario that the combined effect of spontaneous symmetry breaking and spin dynamics lead to the observed spin fluctuations.

As a conclusion, our study reveals a strong interplay between Bose condensation and spin dynamics, which is of particular relevance when spin-dependent and spinindependent interactions take place on a similar timescale (in contrast to previous studies with alkali atoms, see [31]). This interplay can for example result in a delay in obtaining a BEC in spin-excited states or alternatively to the production of weak metastable spinor gases which decay due to spin-exchange collisions. Our results also show that the difficulty to fully thermalize the spin degrees of freedom is a prominent effect to be taken into account for very large spin systems (such as Dy [15] and Er [14]), where all spin states must be saturated for a stable multi-component BEC to be produced. Finally, we point out that when a multi component BEC is dynamically produced, spontaneous symmetry breaking leads to independent phases within the BEC components which triggers spin fluctuation.

FIG. 1

 1 FIG. 1: a) Experimental sequence, showing the reduction in the ODT intensity in a duration tS. An absorption image is taken after a time t and Stern-Gerlach separation. b) Simple cartoon of the evolution of the momentum distributions (p(k)) of atoms in the three lowest spin excited states, illustrating the difficulty of achieving BEC (peak on top of the broad thermal k distribution) in spin excited states due to spin dynamics. Absorption pictures showing: (c) a BEC in ms = -3 and a thermal gas in other spin states for a gas initially prepared with magnetization M = -2.5±0.25; d) a small multi-component BEC for M = -2.00 ± 0.25.

FIG. 2

 2 FIG. 2: a) Number of thermal atoms in ms = -3 (black diamonds) and ms = -2 (red disks) as a function of time t for a shock cooling time tS = 500 ms. b) Corresponding condensate fractions in ms = -3 (black diamonds) and ms = -2 (red disks). The shaded region highlights when both ms = -3 and ms = -2 thermal clouds are saturated, but only ms = -3 atoms condense.

FIG. 4 :

 4 FIG. 4: a) Absorption images after shock cooling experiments performed with tS = 50 ms and an initial magnetization M = -2 ± 0.25. A small BEC is present in the three lowest spin states i.e. ms = -3,-2, and -1. The different images illustrate the fluctuations of magnetization of the condensate fraction. b) Numerical results after evaporation, for three initial relative sets of phases. c) Total condensate fraction of the multi spin component gas for t = tS ms and M = -2±0.25 as a function of tS. We observe small multi-component condensates in the three lowest energy states. The solid line guides the eye.

This work was supported by Conseil Régional d'Ilede-France under DIM Nano-K/IFRAF, Ministère de l'Enseignement Supérieur et de la Recherche within CPER Contract, Université Sorbonne Paris Cité (USPC), and by the Indo-French Centre for the Promotion of Advanced Research-CEFIPRA. M.B. and M.G. acknowledge support of the National Science Center (Poland) Grant No. DEC-2012/04/A/ST2/00090.


Mohammed Alloulah
email: alloulah@outlook.com

Airborne Broadband Ultrasound: Wireless Channel and Motion Tracking Algorithms

Airborne broadband ultrasound (ABU) covers a wide stretch of inaudible acoustic frequencies which have been used to provide scalable (i.e. multiuser) and highly robust ranging and localization of static tags in indoor environments. This is useful for mobile computing applications and robotics. However, because ABU signals are wideband and propagate relatively slowly, they are prone to aggressive Doppler distortion as a result of human-scale movements indoors. Additionally, the ABU channel is characterized by pronounced distance-dependent variations. This paper gives a treatment of the ABU channel both theoretically and experimentally, in order to inform the design of ABU tracking systems. Then, Doppler-tolerant ABU tracking algorithms are developed. Their performance is reported, using a direct-sequence (DS) code-division multiple access (CDMA) signalling approach. Successful motion inference is demonstrated in the ABU band.

I. INTRODUCTION

H IGH-resolution tracking of people and mobile devices in buildings is of great significance to a number of application domains for example: a handheld moving map for indoor navigation, "throwing" video content from a mobile phone to a wall-mounted display using gesture-based interaction, and emergency search and rescue and live coordination of personnel at disaster sites. In this paper, the term tracking is used to encompass both the initial position estimation of a person or object, followed by the continual monitoring of movement thereafter.

Accurate, fine-grained indoor tracking has been classically achieved by infrastructure-reliant systems, using ultrawideband (UWB) radio, optical or magnetic sensing [START_REF]See for instance commercial systems[END_REF]. Ultrasonic sensing has also been used for the provision of tracking services indoors [START_REF] Mccarthy | Accessible Ultrasonic Positioning[END_REF]. Compared to radio, the low frequency and propagation rate of airborne ultrasound yields advantages in accuracy, processing requirements, and ease of installation. However most systems to date have used piezoceramic transducers which have a usable bandwidth of less than five kilohertz. Such narrowband ultrasound poses restrictions on the real-world utility of these tracking systems, especially in scenarios where multiple users must be tracked concurrently, or environments with aggressive ultrasonic noise.

The use of broadband ultrasound in-air was first shown to afford a number of system-level enhancements to tracking [START_REF] Hazas | Broadband ultrasonic location systems for improved indoor positioning[END_REF]; this was demonstrated using a curved piezo film transducer design [START_REF] Toda | Theory of curved, clamped, piezoelectric film, air-borne transducers[END_REF] with a measured bandwidth of 40 kHz and code division multiple access (CDMA) signalling [START_REF] Hazas | Broadband ultrasonic location systems for improved indoor positioning[END_REF].

Others have since explored performance using MEMS transducer technology (bandwidth 15 kHz), and frequency-hopped signalling [START_REF] Gonzalez | High-Precision Robust Broadband Ultrasonic Location and Orientation Estimation[END_REF]. Another related work investigated the outdoor performance (bandwidth 14 kHz) also using direct-sequence (DS) signalling [START_REF] Álvarez | High reliability outdoor sonar prototype based on efficient signal coding[END_REF]. Close in spirit too, a communications work assessed airborne ultrasound's short-range ability to carry data using a practical bandwidth between 120 to 200 kHz with a bit rate of 200 kb/s [START_REF] Li | Short-range ultrasonic communications in air using quadrature modulation[END_REF]. However, the literature lacks a formal model for the airborne broadband ultrasonic channel, which is needed to better inform the design of future ABU tracking systems. Moreover, coded broadband ultrasonic signals are prone to significant Doppler distortion even in the presence of relatively slow rates of indoor movement, such as walking (1 m/s or faster). Consequently, elaborate tracking applications are best supported if ABU signals can continue to be detected in the presence of Doppler shifting.

This paper has two primary contributions. First, it provides (section II) a new description of the ABU channel. An information theoretic analysis to ABU is given in order to highlight the peculiarities of the channel, followed by statistical characterization using real ABU signals. Second, phase coherent direct sequence (DS) CDMA tracking for moving ABU tags is derived (section III), leading to two algorithmic configurations that belong to different classes of instantaneous frequency estimators. Empirical performance (section IV) of the two tracker configurations is then reported.

II. CHANNEL MODEL

In this section, the development of a channel model for ABU follows closely that of its underwater acoustics kin [START_REF] Stojanovic | On the relationship between capacity and distance in an underwater acoustic communication channel[END_REF].

A. Acoustic Propagation

The overall attenuation that characterizes the propagation of airborne acoustics indoors is given by

A(l, f ) = l β • α(f ) l (1)
A(l, f

) dB = β • 10 log l + l • 10 log α(f ) (2)
where l denotes the distance from the source and f is the signal frequency. As evident from equation [START_REF] Mccarthy | Accessible Ultrasonic Positioning[END_REF], the propagation of sound in free air is governed by two main phenomena: spreading and absorption. The parameter β models the spreading loss pattern, and assumes values between 1 and 2. The second term models the distance and frequency-dependent absorption of airborne acoustics, with the coefficient α(f ).

Formulae for calculating the absorption rate per meter at a single discrete frequency are previously developed [START_REF]International Organization for Standardization, Acoustics-Attenuation of sound during propagation outdoors-Part 1: Calculation of the absorption of sound by the atmosphere[END_REF], [START_REF]Committee S1, Acoustics-Method for Calculation of the Absorption of Sound by the Atmosphere, ser. ANSI S1.26-1995. 1430[END_REF]. Using these formulae, figure 1a illustrates the distance and frequency-dependent absorption coefficient of the airborne acoustics of interest. Always calculated at one standard atmospheric pressure, nine reference temperature and humidity permutations are shown in order to give a flavour for the behaviour of the absorption coefficient indoors. Typically, with a temperature of 25 • Celsius and a humidity of 40%, the absorption coefficient can be conveniently approximated at by

10 log α(f ) = 0.03592 × 10 -3 f (3)
as indicated by the dashed line in figure 1a.

B. Signal-to-Noise Ratio 1) Noise generated by receiver: The broadband receiver exhibits variability in terms of the ever present thermal noise across its frequency range. In order to account for this effect in the overall channel model, the receiver's ambient noise is next characterized using a dataset collected from a typical office environment.

Six seconds of data are collected at a sampling rate of 500 kHz. The power spectral density (PSD) of the receiver's ambient noise is estimated using the Welch averaged modified periodogram method. The PSD is referenced to the AC power density of receiver's bandwidth center (50 kHz), which is also the average AC power density of the recorded six-second sequence.

Figure 1b shows that a gradual decay of around 10 dB across the ultrasonic band characterizes in general the ambient noise PSD at the receiver. Outstanding ultrasonic noise spikes can be radiated by everyday indoor machinery, such as air conditioning fans. This is evident in figure 1b with the group of narrowband harmonics just below 50 kHz. While such ambient ultrasonic noise will vary subject to the specifics of a given indoor environment, the general noise behaviour of the broadband receiver can be approximately modelled as

N (f ) dB = 10 log N (f ) = 7.0560 -0.1278 × 10 -3 f (4)
where f is the frequency in Hz and the noise in decibels is referenced to 50 kHz as pointed to above. The noise model is illustrated by the dashed line in figure 1b.

2) SNR:

The expressions for the attenuation A(l, f ) and noise N (f ) can be combined to arrive at a frequency and distance-dependent signal-to-noise ratio (SNR) for the ABU channel, similar to underwater acoustics (UWA) [START_REF] Stojanovic | On the relationship between capacity and distance in an underwater acoustic communication channel[END_REF]. In this ABU formulation, the SNR does not factor in explicitly ambient noise-the accurate ambient noise PSD is dependent on a particular indoor environment. Nonetheless, it still provides a very informative account on the frequency and range variability of the SNR in the ABU channel. This variability translates into range-dependent bandwidth given a constant transmission power. Sophisticated power control in ABU should take advantage of this model in order to optimize bandwidth and capacity in a mobile medium with non-equidistant co-located users. Specifically, the most comprehensive power control could employ pre/de-emphasis techniques using this frequency and distance-dependent SNR model along with estimates of the range and speed of mobile users.

Defining a narrow band of frequencies ∆f around a single frequency f , the SNR is expressed as

SNR(l, f ) = P A(l, f )N (f )∆f ( 5
)
where P is the power of the signle frequency f .

From expression (5), we note that the variability in SNR is characterized by the factor 1/A(l, f )N (f ). Thus, this frequency and distance-dependent factor can be used in order to characterize SNR in the ABU channel. Figure 1c depicts such a factor evaluated at eight distances across the frequency span of interest, with a spherical spreading parameter of β = 2.

Near and far distances give rise to the most variability across the the wide frequency band of interest ([30 kHz, 70 kHz] in our system). At 1 meter, around 3.67 dB variability is in favour of higher frequencies, whereas at 8 meters, about 6.39 dB is in favour of lower frequencies. At 3 and 4 meters, SNR varies little across the frequency band to around 0.8 and 0.64 dB, respectively.

The A(l, f )N (f ) product is linear in dB with a slope that is dependent on distance. This fact gives rise to an optimal frequency f o (l) at which the narrowband SNR is maximized. Such optimal frequency is plotted versus distance in figure 2. Next we will elaborate on one optimization criterion by means of which the capacity of a mobile ABU channel with nonequidistant co-located users could be improved. Fig. 2: The optimal frequency that maximizes SNR.

C. Bandwidth and Capacity

In a mobile ABU tracking system, the capacity of the channel influences greatly the system's quality of service. That is, both the tracking update rate and number of participating users are increased with the enhancing of the channel's capacity.

Assuming a linear time-invariant channel during some observation interval, the overall capacity can be seen as the summation of a number of sub-narrowband capacities. That is, if the bandwidth is divided into I narrow sub-bands each of width ∆f and a noise PSD N (f i )-noise can be approximated as Gaussian-then the overall capacity is given as

C(l) = I i=1 ∆f log 2 1 + S l (f i )A -1 (l, f i ) N (f i ) (6)
where S l (f i ) is the i th frequency bin PSD of the transmitted signal designed to propagate over distance l. The optimal PSD S l (f ) which will maximize the capacity is given by [11, p. 745]

S l (f ) + A(l, f )N (f ) = K l (7)
where K l is a constant related to the total power constraint P (l), and S l (f ) ≥ 0. In turn, the power P (l) is related to SNR as

P (l) = B(l) S l (f )df (8) SNR l, B(l) = B(l) S l (f )A -1 (l, f )df B(l) N (f )df (9)
Varying K l according to the optimization criterion of equation [START_REF] Li | Short-range ultrasonic communications in air using quadrature modulation[END_REF] provides a means for achieving an SNR objective SNR 0 and a pre-specified transmission power. This can be seen by substituting equation [START_REF] Li | Short-range ultrasonic communications in air using quadrature modulation[END_REF] into equations ( 8) & ( 9) so as to obtain

P (l) = K l B(l) - B(l) A(l, f )N (f )df (10) SNR l, B(l) = K l B(l) A -1 (l, f )df B(l) N (f )df -1 (11)
The details of a numerical procedure for arriving at S l (f ) that maximizes capacity are supplied by Stojanovic [START_REF] Stojanovic | On the relationship between capacity and distance in an underwater acoustic communication channel[END_REF]. After such optimization, the distance-dependent transmission PSD is given by

S l (f ) = K l -A(l, f )N (f ), f ∈ B(l) 0, otherwise (12)
with the overall channel capacity being calculated from

C(l) = B(l) log 2 K l A(l, f )N (f ) df (13)

D. Numerical Model

The capacity-based model given previously is evaluated numerically next. As stated before, the bandwidth of interest is between 30 kHz and 70 kHz-outside this range, analysis is meaningless due to physical transduction limitations. A typical indoor temperature of 25 • Celsius and a humidity of 40% are chosen. A spherical spreading factor of β = 2 is used. The SNR objective SNR 0 is swept from -13 dB to 20 dB in 3 dB increments.

The effect of SNR on bandwidth allocation per distance is examined in figure 3a. The different slopes of the A(l, f )N (f ) product (figure 1c) result in a large bandwidth allocation variance across distances when operating in a low SNR regime. This is then gradually enhanced with the increasing of SNR until bandwidth becomes almost uniform starting at SNR 0 = -2 dB (not shown), which is in line with Shannon's limit on reliable communication at -1.6 dB [11, p. 207].

The capacity follows a similar pattern to that of the bandwidth as shown in figure 3b. The capacity peaks at the middle of the distance range in the low SNR regime owing to the almost constant A(l, f )N (f ) product. At SNR 0 = 4 dB, the capacity takes on a consistent shape whereby the channel changes little at farther distances. It is also observed that beyond SNR 0 = 1 dB, enhancements brought about with the increasing of SNR tend to diminish in comparison with those in the low SNR regime.

The power in micro Pascals is depicted in figure 3c. Affected by the bandwidth and capacity observations, around 40 dB transmission is required to clear the low SNR regime. Under the optimal capacity definition of bandwidth, this may or may not be possible to accommodate depending on a number of contending factors such as sound exposure safety, transducer excitation limitations, and operational volume. The 3 dB increments in SNR objective steadily offset the power curve in the high SNR regime.

A utilization measure of the available bandwidth is the bandwidth over capacity B/C ratio. To investigate the effect of SNR on the the utilization of the ABU channel, the B/C ratio at various distances is plotted in figure 4 against power and against SNR in 4a and 4b, respectively. Expectedly, the B/C ratio increases with power and SNR.

E. Time Variability

Away from the mobility-induced Doppler effect, the airborne acoustic channel indoors is influenced by factors such as temperature and humidity. Other contributors also affect the operation of the acoustic medium such as air conditioning drafts and noise emanating from everyday human activities. These factors are hard to model accurately in dynamic living environments and assessing their aggregate effect statistically remains the only viable approach at understanding the seemingly random changes in the airborne acoustic channel.

To this end, a transmitter-receiver pair was placed at an approximate distance of 2.25 m. The receiver was carefully oriented next to a smooth hard surface as to detect two path arrivals, direct and reflected. A Gold code sequence of length 511 was BPSK-modulated onto a 50 kHz carrier at a 20 kHz chip rate. A continuous transmission of about 18 seconds was recorded. The channel impulse response as obtained by correlation is shown in figure 5. The channel consists of a direct path and a reflected path which is always more attenuated as a result of the loss in acoustic energy. The variations over time in the magnitude and phase of the two path arrivals are also illustrated in figure 5, along with the respective histograms. These histograms suggest a Rician or even Gaussian distribution for the magnitude of two arrivals; while the phase seems to follow a uniform distribution. Nonetheless, further investigation is needed before formalizing the statistical nature of multipath in ABU under dynamic environmental conditions indoors.

III. TRACKING

This section focuses on direct-sequence spread spectrum (DSSS) modulation. DSSS allows for explicit self and multiuser-interference models that are computationally tractable, and as such highly amenable to embedded real-time realizations. Such real-time realizations were investigated by Alloulah using a high-level synthesis (HLS) methodology [START_REF] Alloulah | Real-Time Tracking for Airborne Broadband Ultrasound[END_REF].

A. DSSS System Model

The complex baseband signal at the receiver is modelled as

v(t) = m q(m)h(t -mT c ) × e jθ(t) + w(t) ( 14
)
where q is the spreading code, h(t) is the channel impulse response which accounts for transmitter and receiver filtering, T c is the chip duration, θ(t) is a phase function modeling the aggregate offset and/or distortion resulting from transmitterreceiver oscillator mismatch and/or mobility respectively, and w(t) is an additive noise which is assumed to be uncorrelated with the signal. At the receiver, the signal is oversampled by N s samples per chip and arranged in a column vector spanning the spreading code length L and whose topmost entry is the most recent input sample. The received signal vector expressed as the spreading code convolved with the channel can now be rewritten in vector format as

v(k) = m h[m]q[k -m] • e jθ(k) + w(k) (15)
where

v(k) =    v L-1 (k) . . . v 0 (k)    , h(k) =    h L-1 (k) . . . h 0 (k)    = h[0] (16)
are LN s × 1 vectors whose entries v l (k) and h l (k) are partial vectors containing the l th chip-worth of samples and the fractional channel coefficient, respectively. Further, noting the causality of the system, h[m] in the convolution kernel is defined as h(k) shifted down by m × N s samples.

h[m] =         0 . . . m h L-1 (k) . . . h m (k)         (17)
θ(k) =    θ(k + L -1) . . . θ(k)    (18)
Similarly, q[•] denotes modulo-L addressing of the spreading sequence. Also, θ(k) is a vector containing L-chip phase deviations that are multiplied element-wise with their respective L v l (k) as denoted with the dot operator • in equation [START_REF] Sharif | A computationally efficient doppler compensation system for underwater acoustic communications[END_REF]. Lastly, w(k) is an LN s × 1 vector of additive noise.

B. Receiver Algorithm

The proposed receiver structure is one that is user-oriented, decoupled, and extensible. A block abstraction of the receiver is presented in figure 6.

The derivation will proceed in a step-by-step manner to build hybrid Doppler-tolerant trackers. This approach is important owing to the absence of prior experience on phase tracking in the novel ABU band, hence the need to consider different methods in order to assess their utility in ABU and its applications.

1) Core DS CDMA Adaptive Engine: This subsection along with subsection III-B3 follow closely the derivation of Stajanovic et al. [START_REF] Stojanovic | Multiuser code acquisition in multipath channels[END_REF]. Later subsections will depart-building in places on the intuition from [START_REF] Stojanovic | Efficient processing of acoustic signals for high-rate information transmission over sparse underwater channels[END_REF] and [START_REF] Sharif | A computationally efficient doppler compensation system for underwater acoustic communications[END_REF]-to introduce the feedforward filter and linear interpolator stages developed for ABU.

a) Channel-based formulation: The acquisition algorithm is based on the fact that noise is uncorrelated with the spreading sequence

h(k) = E v(k) • e -jθ(k) q * (k) (19)
leading to a simple stochastic approximation of the form:

ĥ(k) = λ ch ĥ(k -1) + (1 -λ ch )v(k) • e -jθ(k) q * (k) (20)
where λ ch is an exponential forgetting factor. For the i.i.d. Gold code sequences, optimal tap selection through truncation in magnitude is readily achieved [START_REF] Stojanovic | Efficient processing of acoustic signals for high-rate information transmission over sparse underwater channels[END_REF]. b) Recursive interference-free signal reconstruction: Under inter-chip-interference (ICI)-be it caused by the timevarying Doppler distortion w.r.t chips within an entire code, or due to multipath-equation ( 15) can be expanded as

v(k) =   h[0]q[k] + m =0 h[m]q[k -m]   • e jθ(k) + w(k) (21)
The ICI signal is defined as

v ICI (k) ∆ = m =0 h[m]q[k -m] • e jθ(k) (22)
and the ICI-free signal as

v 0 (k) ∆ = h[0]q[k] • e jθ(k) + w(k) (23)
The composite signal, not counting the noise, can thus be expressed as

v(k) ∆ = L-1 m=0 h[m]q[k -m] • e jθ(k) = v ICI (k) + h[0]q[k] • e jθ(k) (24)
Adding and subtracting a h[0]q[k] • e jθ(k) term from the RHS of equation [START_REF] Sayed | Adaptive Filters[END_REF] gives:

v(k) = v 0 (k) + v ICI (k) = v 0 (k) + v(k) -h(k)q[k] • e jθ(k) (25)
In addition, v(k) obeys a shifting law such that where

v(k) = vL-1 (k) ↓ Ns v(k -1) (26)
vL-1 (k) = L-1 m=0 q[k + m]h m (k) (27)
and ↓ Ns v(k -1) symbolizes down-shifting by N s samples.

Replacing the channel vector by its estimate, we arrive at a channel adaptation procedure that can be performed in two coupled steps. Firstly, equation ( 25) is rearranged so as to allow the reconstruction of an estimate of the ICI-free signal:

v0 (k) = v(k) -v(k) + ĥ(k)q[k] • e j θ(k) (28)
Secondly, using the ICI-free signal in [START_REF] Gitlin | An algorithm for the stable operation of a digitally-implemented fractionally-spaced adaptive equalizer[END_REF] for better channel estimation yields

ĥ(k) = λ ch ĥ(k -1) + (1 -λ ch )v 0 (k)q * [k] (29)
Once the channel estimate and the ICI-free signal have been computed, the chip estimate is defined by the filtering operation

q(k) = 1 E ĥ(k) ĥ (k)v 0 (k) (30)
where E ĥ(k) = ĥ (k) ĥ(k)

This model allows us to monitor the state of selfinterference explicitly. It can be shown [START_REF] Stojanovic | Multiuser code acquisition in multipath channels[END_REF] that this formulation can be expanded in the multiuser case to account for multiple access interference (MAI), albeit at the expense of increased computation complexity and memory requirements. Specifically, with a total of U users present in the system, the recursive reconstruction of the average signal vector becomes

vL-1 (k) = U u=1 vu L-1 (k) (32)
where

vu L-1 (k) = L-1 m=0 q u [k + m]h u m (k) (33)
with the superscript u denoting the u th user's reconstructed signal, spreading code, and channel vector.

c) Acquisition test: Monitoring the MSE of a running despreader provides a measure of signal-to-noise-plusinterference ratio (SNIR) at the output:

d(k) = E{q(k)q * (k)} = 1 L k m=k-L+1 q(m)q * (m) (34)
e d(k) = 1 -d(k), SNIR out ∼ E |e d(k)| 2 (35) SNIR out (k) = -10 log 10 1 N d k m=k-N d +1 |e d(m)| 2 (36)
where N d is the length of despreader observation.

Once the MSE drops below a pre-specified threshold, it can be used as an indication for successful acquisition.

2) Implicit DFE Extension: Applying a decision feedback equalizer (DFE) to v(k) at the chip-rate, as if no spreading is taking place, a chip estimate can be expressed as

q(k) = a v(k) -b q(k -1) (37)
where a and b are the feedforward and feedback coefficient vectors respectively, and q(k -1) is the vector of the previously searched chips (during acquisition) whose length is in accordance with distortion span in chips M .

q(k -1) = q[k -1] • • • q[k -M ] T (38)
If we substitute our signal model ( 15) into (37) and neglect the noise for the moment, we obtain

q(k) = a ĥ[0]q[k] + a m>0 ĥ[m]q[k -m] • e j θ(k) -b q(k -1) (39)
Inspecting equation (39

) implies that b = a m>0 ĥ[m] • e j θ(k) (40)
is required so that the ICI contribution to the chip estimate be eliminated. In other words, the task of backward filtering is readily taking place through shifts of the channel estimate, where we assumed that the channel is zero outside the multipath span, i.e.

ĥn (k) = 0, n > M (41)
Equally, we can substitute (25) into (37) to obtain

q(k) = a v0 (k) + a vICI (k) -b q(k -1) (42) = a v0 (k) + a v(k) -a ĥ(k)q[k] • e jθ(k) -b q(k -1) (43)
Comparing equations ( 43), (39), and (28), one notices that it is possible to further filter the incoming signal v(k) with a linear adaptive equalizer a (k) independent of ĥ(k) in order to obtain an improved chip estimate at time kT c . The clean chip estimate is defined by the term containing the reconstructed ICI-free feedforward signal-which is dependent upon the goodness of the channel estimate

q(k) = a (k)v 0 (k) (44)
The combined operation effectively realizes an implicit DFE wherein the separation of the incoming signal into forward and backward components is denoted by the subscripts f and b respectively. Equating (44) and ( 30) is interesting:

a (k) = 1 E ĥ(k) ĥ (k) (45)
In the absence of multipath, the sole role of a (k) would be to combat phase distortion in the fractional chip fed into the spread spectrum algorithm. Equation (45) tells us that if the channel alone fails at equalizing the chip estimate perfectly, the feedforward filter of the implicit DFE would help in accomplishing the task by eliminating residual distortion. This is especially the case for a highly-Doppler-susceptible, nonstationary indoor acoustic environment.

To this end and with the intuition of equation (45) in mind, the modified fractional feedforward filter (FF) is made to operate on incoming chips at the oversampling rate with N s samples per chip, although its adaptation remains at the chip rate such that

v(t) = L ff m a(mT c /N s )s(t -mT c /N s ) ( 46
)
where s(t) is instead the raw complex baseband signal at the receiver now. In vector format, we define

c(k) ∆ = c Ns (k) • • • c 1 (k) T ( 47
)
to be an oversampled chip of N s entries at time kT c , and

v(k) ∆ = a (k)      c Ns (k) • • • c 1 (k) . . . . . . ↓ 1 s(k -1) c 1 (k) . . . ↓ Ns s(k -1)      T (48)
to be the modified fractional feedforward filtering at time kT c . The length of the filter is designed to span a certain number of oversampled chips L ff , and is devoted to combating phase distortion by means of matching the filter response to the incoming signal. Inevitably, it also introduces ICI which will be leaked away by the sparse channel-based backward filtering and futher combated within the interference cancellation in equation ( 43).

The LMS-adapted FF filter is given by

â(k + 1) = â(k) + µ LMS e * q (k) e -j θ(k) s(k) (49)
3) Integrated Second-order PLL: The signal baseband model ( 15) includes an explicit term to account for transmitterreceiver phase variations at the chip rate. While these variations are due to mismatch and/or node mobility, they also affect the order-comparable carrier and code frequencies. This phase term is estimated using a second-order, stochastic gradient phase-locked loop (PLL) according to the MMSE criterion of the chip estimate as follows.

The chip estimate error is first generated

e(k) = q[k] -q(k) ( 50
)
The gradient is then shown to be [START_REF] Stojanovic | Phase-coherent digital communications for underwater acoustic channels[END_REF] ∂|e

2 (k)| ∂ θ = 2 Re ∂e(k) ∂ θ e * (k) = -2 Im 1 E ĥ(k) ĥ (k) v(k) • e -j θ(k) e * (k)
(51) where v(k) is the chip-by-chip, adaptively forwarded vector. The PLL is then implemented as

ψ(k) = Im 1 E ĥ(k) ĥ (k) v(k) • e -j θ(k) e * (k) θ(k + 1) = θ(k) + K 1 ψ(k) + K 2 m≤k ψ(k) (52)
where K 1 and K 2 are the loop constants and often chosen such that K 2 = K 1 /10. 4) Linear Interpolation Stage: Under the narrowband assumption, the Doppler effect can be sufficiently modelled as a frequency shift. On the contrary, in wideband signals the ideal model is better represented by an accurate rate change operation [START_REF] Sharif | A computationally efficient doppler compensation system for underwater acoustic communications[END_REF]

r(nT s ) = s (1 + δ)nT s (53)
The removal of the non-uniform Doppler shift across the spectral content of the wideband signal is then performed as a sampling frequency rescaling operator

f s = (1 + δ)f D s (54)
where f D s is the original sampling frequency of the receiver at which the Doppler effect is manifest.

The framework afforded by the DS CDMA adaptive engine derived earlier facilitates a front-end linear interplator (LI) stage, also operating at the chip rate, that can seamlessly plug into the global cooperative MSE minimization criterion.

To this end, we define a time-varying stage that linearly interpolates the incoming complex baseband signal in accordance with a cost function driven by the core DS CDMA demodulator. In order to minimize the distortion introduced by linear interpolation, the complex baseband signal is well oversampled and decimation occurs tacitly, i.e. all oversamples participate in interpolation.

In a decision-directed mode, the maximum likelihood (ML) carrier phase estimator for the signal s(t; φ) with a known information sequence {I n } is given by [11, p. 304]

φML = -tan -1 Im K-1 n=0 I * n y n Re K-1 n=0 I * n y n ( 55
)
where y n is the output of the matched filter in the n th signal interval T , and K is a positive integer defining the observation interval T 0 = KT . Therefore, in a similar fashion [START_REF] Sharif | Adaptive doppler compensation for coherent acoustic communication[END_REF], the interpolation factor is adapted per step as

I(k + 1) = I(k) + K p φ(k) (56)
where K p is a phase tracking constant, and φ(k) is the ML phase estimate. Within our code acquisition DS CDMA formulation, φ(k) is updated at the chip rate using the chip phase as opposed to the symbol phase

φ(k) = ∠q(k) (57)
C. Motion inference 1) First-order moment-Velocity: Velocity inference is accomplished through the estimation of the the instantaneous frequency. The procedure for estimating IF varies depending on the tracking configuration.

a) DFE-PLL: In the joint DFE-PLL mode, there is an apparent phase tracking redundancy. In underwater acoustic (UWA) systems, PLL's tracking constants are made substantially higher than the FF's learning rate-which is characterized by its length and percentage of error feedback-such that the PLL dominates the FF equalizer [START_REF] Johnson | Improved doppler tracking and correction for underwater acoustic communications[END_REF]. In the ABU band, a special case arises when a competition between the feedforward complex filter and the second-order PLL is allowed. When both attempt to compensate for Doppler-induced phase distortion in near-equal strengths, it was observed empirically that the unwrapped phase of the chip estimates is no longer a monotonically increasing function. Rather, the unwrapped phase possesses some notable attributes that are correlated with motion:

1) It tends to inflect upon reversal of acceleration.

2) It drifts during motion pauses with a certain slope.

3) Provided that the responsiveness of the DFE-PLL configuration is tuned accordingly, the range covered by the unwrapped phase, in a time interval of continuous motion, is very close to the limit of the speed in that motion segment. This can be ascertained by converting the unwrapped phase into a unitless turns-per-second entity according to the direct weighing

Ω i (k) = c f carr 1 2π φ chip (k) (58)
where Ω i is the weighted instantaneous phase, measuring the turns per second that the phase makes throughout the course of the periodic motion stimulus. 4) The shape of the unwrapped phase is also a function of the responsiveness of this competition. These observations seem to suggest that the unusual competition between the PLL and DFE in tracking the phase results readily in a measure of the instantaneous frequency.

However, the low ABU speed of propagation is even less than one fourth that of the underwater environment. The timedilation or compression that airborne acoustics undergo is of the same order (up to a tenth of the center frequency) often encountered in RF inter-satellite communications. This is best visualized by examining the Doppler shift formula f dv ↑= f v↑ c↓ that translates increased velocity into accentuated Doppler shift. Moreover, the order-comparable chip rate further exasperate the Doppler distortion.

Thus the ultra high velocities in the ABU band would strain even a combined DFE-PLL equalizer being updated at the chip rate. Among the symptoms of such strain are excessive equalizer tap rotation and possible phase tracking loop instability [START_REF] Sharif | Closed loop doppler tracking and compensation for non-stationary underwater platforms[END_REF]. This strain calls for leaky adaptive techniques [START_REF] Gitlin | An algorithm for the stable operation of a digitally-implemented fractionally-spaced adaptive equalizer[END_REF] for FF as to aid stability, tracking, and combat drift.

Building on adaptive filtering concepts [21, p. 384-387] and techniques [START_REF] Nascimento | Unbiased and stable leakage-based adaptive filters[END_REF], we introduce the following ad hoc perturbation to the standard LMS adaptation criterion. We dub this arrangement Doubly Reinforced Code Epoch (DRCE) and define it as follows: Once every code period, the last chip is skipped and replaced by the first incoming chip which is consequently used to generate the adaptation error. The next iteration proceeds normally at the first chip again. This effectively injects data-dependent perturbation to weights. The modification introduced to the circular Gold sequence at time kT c is simply

q(k) = q[k] • • • q[k -L + 2] q[k -L] T (59) where q[k] = q[k -L].
From empirical grounding, it was found that this ad hoc technique has a critical effect on Doppler tracking for the DFE-PLL configurations.

A less formal inspiration for DRCE leaky adaptation can be found in the domain of pedestrian tracking using inertial sensors. According to the so-called zero-velocity updates (ZUPT) criterion, at the end of each estimation cycle, the velocity is set to zero upon detection of a stance phase in order to prevent excessive error growth in the navigation model [START_REF] Foxlin | Pedestrian tracking with shoe-mounted inertial sensors[END_REF]. This is vaguely related to DRCE-observing the periodicity of backto-back code transmissions, DRCE injects a cyclic perturbation to weights on a code-by-code basis. However, in DRCE, the perturbation is derived from the incoming data too. b) DFE-LI: Despite LI's agility in tracking Doppler, distortion introduced by linear interpolation needs to be dealt with. The matched filtering properties of the FF readily provides an equalizing effect for a variety of distortion sources. Therefore, the LI-DFE combination represents a very viable phase tracking mechanism whereby:

1) The interpolation factor's polarity determines the direction of movement: backward or and forward.

I sgn (k) = sgn I(k) (60)
I sgn is then used to delimit the motion segment and velocity estimation is performed as

v i (k) = c f carr 1 2π φ u chip -φ l chip t u -t l (61)
where the superscripts u and l denote respectively the upper and lower bounds of the motion segment as determined by I sgn . 2) Minimal tap translation and rotation per motion segment takes place in FF since LI will dominate phase tracking. 3) Acceleration and/or jerk can tip the balance of LI such that the interpolation factor will not return to rest. In this case, the FF will intervene briefly to adjust the phase resulting in ambiguity in the measured rate of change from the slope of the unwrapped phase of the chip estimate. 2) Second-order moment-Acceleration stress: The adaptive estimation of instantaneous acceleration (IA) is desirable. This desirability is because online mechanisms for the inference of additional physical properties of motion (e.g. acceleration) could add to the richness of the sensory model and allow for sophisticated workarounds in the presence of higher motion moments in the ABU band.

The IA estimator for coherent ABU is developed here in analogy with a single-tone, frequency-modulated (FM) signal. Specifically, we note that the running despreader d(k) in equation ( 34) can be thought of as a phasor at any given time kT c . This phasor represents a measure of how much, after demodulation, the receiver code is in-lock with the undistorted transmitted code. The rate of change in the demodulated despreader phase is now associated with the second motion moment, acceleration. That is, the complex depreader is already a first-order statistical moment of the chip estimates, which in turn can be seen as a complex-valued random variable (note the expectation operator in equation ( 34)). The despreader phasor can therefore be modelled in the continuoustime domain as

z d(t) = a d(t) e jφ d (t) = x d(t) + j y d(t) (62)
with

φ d(t) = 2π t -∞ f i (τ ) dτ (63)
and f i (t) denoting the time-varying, second-order instantaneous frequency.

Differentiating the despreader phase, using the classic continuous-time formula for FM discriminators [START_REF] Harris | Exact FM detection of complex time series[END_REF], yields the discrete-time estimator

f i (k) = 1 2π x d(k)y d(k) -x d(k)y d(k) x 2 d(k) + y 2 d(k) (64)
For a convergent d(k) i.e. after acquisition has been declared and tracking has begun, a number of approximations hold. First, if Re{ d(k)} is stable, then x d(k) ≈ 1, x d(k) ≈ 0, and x 2 d(k) ≈ 1. Second, if phase lock is maintained, y 2 d ≈ 0. These approximations lead to the following second-order IF estimator in our ABU system

f i (k) ∝ 1 2π y d(k) (65)
The second-order IF component is proportional to the rate of change (first derivative) of the imaginary part of the running despreader. In expression (65), equality is substituted for proportionality bearing in mind the limited code-length average which is provided by the running despreader. That is, additional processing (such as exponential statistical weighting) will be required if we were to convert this unitless entity into a meter-per-second-squared instantaneous acceleration estimator.

The final IA estimator is then obtained by scaling the second-order IF component according to

a i (k) = c f carr f i (k) (66)

IV. EXPERIMENTAL RESULTS

A. Signal Design and Experimental Setup

A Gold code of length 511 was BPSK-modulated onto a 50 kHz carrier at a 20 kHz chip rate. Chips were pulseshaped using a raised-cosine filter with roll-off factor of 1. The sampling rate utilized was 160 kHz.

One data set will be used for the characterization of the two algorithmic configurations. This data set has typical maximum speed (about 0.8 m/sec) and contains high-order motion moments. In order to facilitate characterization, we utilize a vision-based tracking system for the generation of real-time position and speed ground truth. A transmitter is mounted on a Lego Mindstorms robot moving forward and backward on a track while facing a stationary ABU receiver, at a distance of 2.9 m. The robot is tagged with a fiducial marker [25]. A PC-based application monitors the movement of the robot using the computer vision-based tracking. Upon motion, the application triggers an FPGA-based ABU transmitter and simultaneously commences data acquisition. The all-digital transmitter has an independent oscillator to the PC clock. The experimental setup used for data collection is shown in figure 8.

B. Performance

The acquired signal was down-converted to baseband where all processing takes place. Figure 7 illustrates the performance of the DFE-LI and DFE-PLL tracking configurations. The upper set of plots belong to DFE-LI. The first plot shows the channel which is maintained throughout tracking. The adaptation mean-square error (MSE) begins its decent initially until tacking reaches a steady-state around 1 sec. The second plot depicts of the evolution of the linear interpolation factor in time and the chip estimate scatter. The sensitivity of the DFE-LI tracker to acceleration and jerk was avoided initially by a priori knowledge of the Doppler stimulus. As such, the instance at which LI is activated was delayed. As shown, DFE-LI tracks Doppler very accurately during the first and most of the second motion legs. The velocity for the first motion segment was estimated to be 0.8047 m/s, which is virtually identical to ground truth. However, following the decelerating motion moments in the second leg, LI diverges away from rest condition owing to high second-order motion stresses. The third plot examines this stress resulting from second-order motion moments. The IA estimator shows distinct sustained acceleration spikes in three places: at the beginning of the Doppler stimulus around 0.5 sec; around the momentary pause of the robot at roughly 3.5 sec between leg 2 & leg 3; and at the end of the test case slightly before 7 sec. These instances are in line with the acceleration groundtruth of the Doppler stimulus.

Similarly, the performance of the DFE-PLL tracking configuration is illustrated in the second row of plots of figure 7.

Justifiably, the crude DFE-PLL tracker is noisier even though code timing stays in-lock. In the second plot the PLL phase evolution is shown alongside the chip estimate scatter, which is considerably more spread in comparison to that of the DFE-LI configuration. The third plot shows the scaled unwrapped phase of the chip estimate overlaid on the vision tracking ground truth. The scaled instantaneous phase is evidently correlated with motion in this competitive mode, which provides a crude means for estimating velocity. This method for inferring the motion pattern could be useful in certain applications not requiring precise velocity measurements.

V. CONCLUSION This paper has given a novel characterisation of the ABU channel, providing a summary of expectations for signalto-noise ratio, channel bandwidth and capacity, for typical indoor office conditions. Real-world ABU location systems require Doppler tracking in order to cope with typical indoor movements. This paper has developed the first Dopplertolerant ABU tracker (figure 6). It is based around an implicit decision feedback equaliser, incorporating a feedforward filter which has been modified to operate on incoming chips at the oversampling rate. The DFE works in conjunction with two other blocks: (1) an integrated second-order PLL to estimate phase variations (due to mismatch or node mobility); and (2) a linear interpolator which models Doppler shift as a rate change operation. Methods for computing the first-and second-order moments (velocity and acceleration) were also presented.

The ABU Doppler tracker is thus an hybrid one, operating as DFE-PLL or DFE-LI. Experimental analysis shows that the DFE-PLL tracker stays locked on, but produces noisier chip estimates and can only give a crude indication of velocity based on its unwrapped, scaled instantaneous phase. By contrast, the DFE-LI variant can produce highly accurate velocity estimates, but tends to lose track in the presence of significant acceleration and jerk.

  l, f )N (f )] -1 (dB) (a) Acoustic absorption in air (dB/m). (b) PSD of the receiver noise floor. (c) The A(l, f )N (f ) product.

Fig. 1 :

 1 Fig. 1: The dependency of SNR in the ABU channel on frequency and distance.

Fig. 3 :Fig. 4 :

 34 Fig. 3: The effect of SNR on (•) versus distance.

Fig. 5 :

 5 Fig. 5: Multipath intensity profile.

Fig. 6 :

 6 Fig. 6: The overall tracker for ABU.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Tracking performances: upper row DFE-LI and lower row DFE-PLL.