Particle Filter Design using Importance Sampling for Acoustic Source Localisation and Tracking in Reverberant Environments

—Sequential Monte Carlo methods have been recently proposed to deal with the problem of acoustic source localisation and tracking using an array of microphones. Previous implementations make use of the basic bootstrap particle ﬁlter, whereas a more general approach involves the concept of importance sampling. In this paper, we develop a new particle ﬁlter for acoustic source localisation using importance sampling, and compare its tracking ability with that of a bootstrap algo-rithm proposed previously in the literature. Experimental results obtained with simulated reverberant samples and real audio recordings demonstrate that the new algorithm is more suitable for practical applications due to its reinitialisation capabilities, despite showing a slightly lower average tracking accuracy. A real-time implementation of the algorithm also shows that the proposed particle ﬁlter can reliably track a person talking in real reverberant rooms.


I. INTRODUCTION
The concept of acoustic source localisation and tracking (ASLT) plays an important role in many practical speech acquisition systems.Domains of application include teleconferencing, multimedia information processing, and handsfree telephony, to name but a few.Other applications, such as automatic speech recognition and speaker identification systems, are also very sensitive to the quality of the audio input signals.In most cases, exact knowledge of the speaker position is the key to acquiring clean speech using such tools as beamforming or equalisation principles.
The multipath propagation of acoustic waves in practical environments however constitutes a major challenge to overcome for any tracking algorithm.Recently, methods based on a state-space approach (Bayesian filtering) have been developed to deal with this problem [1][2][3].Because Bayesian filtering algorithms deliver location estimates based on a series of past measurements rather than the current observation only, these methods are more efficient at dealing with the spurious effects of acoustic reverberation than traditional ASLT algorithms.Also, a tracker based on state-space filtering involves a model of the specific target dynamics, providing information regarding how the source is more likely to evolve from one time step to the next.This enables the tracker to effectively discriminate between observations originating from the true target and erroneous observations resulting from acoustic disturbances.
Among the different methods based on Bayesian filtering, the concept of particle filtering (PF) appears as a promising approach to tackle the ASLT problem [2][3][4].As a sequential Monte Carlo method, the PF technique can be used to deal with nonlinear and/or non-Gaussian problems, making it superior to algorithms such as the Kalman filter and its derivatives.This is of particular importance for ASLT, where the observations typically result from a nonlinear process due to the chosen localisation procedure (such as steered beamforming [5], cross-correlation [6] or eigenvalue decomposition [7]).Also, the observation noise in ASLT problems is usually non-Gaussian due to the effects of acoustic reverberation.Particle filtering can then be used to consider several observations per sensor, in order to represent multi-modal density functions reflecting the multiple hypotheses that each of the measurement modalities might originate from the target (see e.g.[2]).
Previous research works on particle filtering applied to ASLT, such as [3,4,8], make use of the basic bootstrap particle filter, introduced by Gordon, Salmond and Smith in [9].The conceptual simplicity of this algorithm leads to straightforward practical implementations and moderate computational requirements.The bootstrap PF however suffers from a major drawback: during each iteration, the particles are relocated in the state space without knowledge of the current observations.The PF might hence omit some important regions of the state space when searching for the target, which mainly precludes the PF from reinitialising after a target disappears or becomes occluded for a short period of time.Despite showing promising results, this algorithm con-sequently still lacks some important characteristics necessary for a smooth operation in practical scenarios, such as the automatic detection of new targets and the ability to recover from track loss.
In this research, we develop a particle filtering method based on the more general concept of importance sampling (IS), in which particles are generated during each iteration on the basis of both the particle set at the previous time step and the current measurement.This provides the resulting algorithm with the important property of reinitialisation.Importance sampling further allows the combination of different types of observations in a global statistical framework.
The development of a robust acoustic source tracking algorithm for reverberant environments is the main motivation behind the research described in this paper.In the next section, we review the generic approach to the problem of ASLT.The basic concepts of bootstrap filtering and importance sampling are briefly explained in Sections III and IV.We then develop a particle filter for ASLT using the IS approach in Section V, and finally present the results of experimental tests that demonstrate the performance of the newly proposed algorithm in Section VI.

II. SOURCE TRACKING AND BAYESIAN FILTERING
Consider an array of M acoustic sensors distributed at known locations in a reverberant environment with known acoustic wave propagation speed c.Assuming a single sound source, the problem is to estimate the location of this "target" for each time step k = 1, 2, . .., based on the signals s m (t), m ∈ {1, . . ., M }, provided by the array.Let X k represent the state variable at time k, corresponding to the position and velocity of the target in the state space:1 At each time step, each microphone in the array delivers a frame of audio signal which can be processed using some localisation technique such as, for instance, steered beamforming (SBF) or time-delay estimation (TDE).Let Y k denote the observation variable (or measurement) which, in the case of ASLT, typically corresponds to the localisation information resulting from this processing of the audio signals.
Using a Bayesian filtering approach and assuming Markovian dynamics, this system can be globally represented by means of the following two equations: where g( The solution to this Bayesian filtering problem consists in the following two steps of prediction and update [9].Assuming that the posterior density p(X k−1 |Y 1:k−1 ) is known at time k − 1, the posterior PDF p(X k |Y 1:k ) for the current time step k can be computed using the following equations: where is the transition density, and p(Y k |X k ) is the so-called likelihood function.

III. BOOTSTRAP PARTICLE FILTER
Particle filtering is an approximation technique that implements the recursion of (2) by representing the posterior density as a set of samples of the state space X (n) k (particles) with associated likelihood weights w A basic PF variant is the bootstrap filter [9] which can be described as follows.Assume that the set of particles and weights {(X n=1 is a discrete representation of the posterior density p(X k−1 |Y 1:k−1 ).The bootstrap PF then implements the following three iteration steps: 1. Resampling: draw N samples X Then normalise the weights so that they add up to unity: As a result, the set of particles and weights {(X is approximately distributed as the current posterior density p(X k |Y 1:k ).The sample set approximation of the posterior PDF can then be obtained via: where δ(•) is the Dirac delta function, and an estimate X k of the target state for the current time step k follows as: The disadvantage of this algorithm is that during the prediction step, the particles are relocated in the state space without knowledge of the current measurement Y k .Some regions of the state space with potentially high posterior likelihood might hence be omitted during the iteration, leading to a decreased tracking performance.This drawback can be addressed using the concept of importance sampling.

IV. IMPORTANCE SAMPLING
Assuming perfect Monte Carlo sampling, let {X be a set of N random samples drawn from the density p(X k |Y 1:k ), with uniform weights w This sample set allows the approximate computation of any statistical quantity of interest based on the PDF p(X k |Y 1:k ) such as its mean or mode, which can be used as an approximation of the current target state.In practice however, the posterior density is not usually available and it is hence impossible to sample directly from it.
An alternative solution is the use of importance sampling (IS), see e.g.[10].This method consists in choosing a so-called importance density q(X k |Y 1:k ) from which particles are easy to sample, X (n) k ∼ q(•).Then, for the approximation in (4) to remain a truthful representation of the desired posterior density p(X k |Y 1:k ), the computation of the weight must be updated to (see e.g.[11]): where the second line follows from (2b).The importance weights are hence defined as the product of the likelihood function and a correction term that compensates for a potentially uneven distribution of the particles that might result from the process of sampling the importance function.The generic IS algorithm can be summarised as follows: 1. Sample N particles according to the importance function, For each particle, compute the unnormalised importance weight as defined in (6): Then normalise the weights according to (3).The set of particles and weights {(X n=1 is then approximately distributed as the current posterior PDF p(X k |Y 1:k ), and an estimate of the current state can be computed using (5).To emphasize the fact that the particles are here sampled according to a specific PDF (rather than propagated from the previous time step as in the bootstrap implementation), the term importance particles will be used from now on to denote the samples X (n) k generated by drawing from the importance function q(•).
Note that, although described in this work as a separate algorithm, the bootstrap PF of Section III corresponds to a special case of the IS algorithm presented here.The bootstrap filter can indeed be derived from the IS procedure with the simplifying assumption q(•) p(X k |X k−1 ), emphasising the fact that particles are sampled without taking the current observations into account.Further information on existing PF algorithms and other Monte Carlo methods can be found in [10][11][12].
The importance sampling principle allows a decreased estimate variance by virtue of an improved sample-based representation.In terms of minimising the variance of the weights, which constitutes the so-called degeneracy problem in PF implementations, the optimal importance density q opt (•) has been shown to be [10]: It can be seen that this choice of importance density takes into account both the previous state X k−1 and the current observation Y k , making the IS algorithm more robust than the bootstrap method.
In theory however, any density (subject to some weak assumptions) could potentially be chosen as importance function, the main purpose of which is to redirect some of the particles in regions of the state space with potentially high posterior likelihood.In previous literature for instance, the importance function q(•) was implemented to take advantage of measurements from auxiliary sensors (see e.g.[13]), which provides an efficient way of fusing data obtained from different observations.Similarly, the algorithm presented in [14] implements the IS method to draw on information obtained from two different measurement processes derived from the same raw data.Contrary to the method consisting in combining the different observations in the representation of Y k , the IS technique hence offers a principled way of including these in a common framework, even when the statistical relationship between the different measurements is not completely known or hard to determine.This specific approach is applied here to the ASLT problem.

A. Algorithm Design
It can be seen that three design choices need to be made for a practical implementation of the IS principle, regarding the definition of the target dynamics, the likelihood function, and the importance function.These issues are discussed in detail below.
1) Target dynamics: in order to remain consistent with previous literature [2,3], a Langevin process is used to model the dynamics equation (1a).This model is typically used to characterise various types of stochastic motion, and it has proved to be a good choice for the current application.The source motion in each of the Cartesian coordinates is assumed to be an independent first-order process, which can be described by the following equation: with the noise variable: where N (µ, Σ) denotes the density of a multi-dimensional Gaussian random variable with mean vector µ and covariance matrix Σ.The parameter T U corresponds to the time interval separating two consecutive updates of the particle filter.The model parameters in (8) are defined as: with v the steady-state velocity parameter and β the rate constant.The transition PDF p(X k |X k−1 ) then simply follows from the noise characteristics defined in this model: with N (α; µ, Σ) the density of a Gaussian variable with mean µ and covariance matrix Σ evaluated at α.
2) Likelihood function: experimental results from previous research carried out on particle filtering for ASLT have shown that steered beamforming (SBF) delivers an improved tracking performance compared to TDE-based methods [3,15].The SBF principle is hence used here to implement a pseudolikelihood (PL) function, as introduced in [3]. 2 With S m (ω) = F{s m (t)} the Fourier transform of the m-th signal data, the likelihood function is defined as the output P Ω (ℓ) of a delayand-sum beamformer (DSB) steered to the location ℓ = [x y] T , and computed over the frequency domain Ω: where ℓ m = [x m y m ] T is the known position of the m-th microphone.In the sequel, the likelihood function is hence computed according to p(Y k |X k ) P ΩL (ℓ), with the location vector ℓ reflecting the current state of the variable X k and with the integration in (9) carried out over the frequency range 3) Importance function: the purpose of q(•) is to relocate some of the particles in the state space taking the current observation into account, and potentially also taking advantage 2 The pseudo-likelihood is defined as a pseudo-density, which differs from a true PDF in that it is not necessarily suitably normalised.The reader is referred to [3,8] for a description of the pseudo-likelihood approach.of a different measurement process.Rather than a fine-scale and accurate representation of the particle sampling areas, the importance function is typically meant to give a coarse indication of where the particles should be sampled in the state space.Based on the signals received at the sensors, several principles could be used to implement this function.The SBF output computed for low frequencies is however known to possess these desired properties.The SBF beampattern at high frequencies generally exhibits a narrow main lobe and suffers from aliasing effects which typically generate spurious peaks in the observations. 3For low frequencies however, the aliasing effects are reduced and the width of the main lobe in the beampattern becomes more important, leading to less accurate but also less ambiguous localisation results.Hence, this approach is of particular interest in the context of importance sampling, and the importance function is defined here as q(•) ∝ P ΩS (ℓ), which is computed according to (9) with the integration carried out over the frequency band Ω S : ω ∈ 2π • [100Hz, 400Hz].Note that because the importance function is typically evaluated on a grid defined across the entire state space (see Section VI-A), this function can be easily normalised and it is hence not defined as a pseudodensity.

B. Proposed IS Algorithm for ASLT
The proposed IS algorithm for ASLT, which will be denoted SBF-IS from now on, is given in Algorithm 1.It must be noted that the previously defined importance function is only a coarse approximation of the optimal density q opt (•) defined in (7), since it only relies on the current SBF measurements.In order to generate some of the state samples on the basis of the previous particle set {X , a standard bootstrap option is included in the algorithm (iteration step C).Also, in a manner similar to [14], the reinitialisation step (iteration option A) has been added to allow the PF to deal efficiently with speech pauses or detect a new target entering the scene.This procedure can be seen as a mixed-state bootstrap step, with particles distributed according to a combination of the original bootstrap density and the reinitialisation density.To this purpose, the reinitialisation density has been simply defined to be the same PDF as the importance function, implicitly defining iteration option A of Algorithm 1 as an importance sampling step without compensation of the corresponding importance weights.
The resampling process involved in iteration step C of the IS algorithm can be easily implemented using a scheme based on a cumulative weight function [9].Alternatively, several other resampling methods are also available from the particle filtering literature, see e.g.[11].Any of these methods may also be used to efficiently implement the process of sampling particles from the (discrete) importance function q(•), in steps A and B of Algorithm 1.
Assumption: at time k − 1, the set of particles and weights {(X n=1 is a discrete representation of the posterior distribution p(X k−1 |Y 1:k−1 ).Iteration: for each particle, i.e. for n = 1, . . ., N, choose randomly one of the following sampling methods according to their respective probabilities: A. Reinitialisation (probability P R ): sample the particle Importance sampling (probability P S ): sample the particle X (n) k ∼ q(X k |Y 1:k ), and compute the unnormalised importance weight according to (6): Finally, normalise the weights according to (3).
k )} N n=1 of particles and weights is approximately distributed as the posterior density p(X k |Y 1:k ), and the current target state can be estimated according to (5).

C. Discussion of Practical Implementation Aspects
The respective probabilities of each sampling method are free parameters in the IS algorithm.They can be determined in various ways, including setting them to constant values, as done in [14].Here, these probabilities are determined at every time step on the basis of whether the current importance function is suitable for sampling or not.Ideally, the importance function is expected to present one peak only, explicitly defining one single region where particles are to be generated.If this function presents several local maxima, it is obviously not appropriate for single-target tracking.Hence, during each PF iteration, the importance function is first computed across the state space, and the number N P of peaks above a certain threshold (defined here as 90% of the largest measured value) is then determined.The reinitialisation and bootstrap probabilities are then computed as: P R = P R /N P and P S = P S /N P , where P R and P S are the prior probabilities of each method, respectively, and have been optimised on the basis of practical tests as P R = 0.01 and P S = 0.25.
In practice, the density p(X in the computation of the importance sampling weights (iteration step B) can be approximated as follows, using (2a) and (4): However, because the importance particles are sampled in the state space in a manner that usually violates the propagation model described by (8), the transition PDF p(X k |X k−1 ) in ( 10) must be updated in order to allow these sampled particles to be given non-zero weights.In the sequel, the following transition PDF will be used in the implementation of (10): where U(•) denotes the uniform distribution (defined over the considered state space), and the background probability ψ is set to a small constant to account for the fact that importance particles are not governed by the same dynamics model as particles used in a standard bootstrap step.More information about tracking models with switching parameters is provided in [17].
Finally, it can be seen that the importance function q(•) defined in Section V-A only contains spatial information about the state vector X k .As a result, the velocity component of the importance particles is here set to some random value upon sampling from the importance density:

VI. PRACTICAL EXPERIMENTS A. Experimental Setup
The setup defined for the following experiments was based on a medium-sized room measuring roughly 2.9m × 3.8m × 2.7m, and fitted with an array of M = 8 omnidirectional microphones positioned at a constant height and organised as one pair on each wall.In each pair, the distance between the sensors was 0.6m.
The microphone signals used in the experiments were samples of audio data sampled at 8kHz, either recorded in a real office room or generated using the image method [18].For the practical recordings, the sound source was simulated with a loudspeaker moving along a predefined path across the enclosure.The signals were split into frames of 512 samples (processed using a Hamming window), and subsequently used as observation to compute both the importance and likelihood functions.The data processing was carried out using a 50% overlapping factor, yielding the update interval T U = 0.032s.The numerical values defined for the transition model parameters were set to v = 0.7m/s and β = 10Hz.
For the SBF-IS algorithm, the importance function was computed over a horizontal grid of points uniformly distributed across the state space with a spacing of 0.1m.
In the following results, the performance of the IS algorithm is compared to that of the SBF-PL method, a bootstrap-only algorithm described in [3].For both methods, the number of particles was set to N = 30.Other algorithm-specific parameters were optimised empirically to achieve a satisfactory tracking performance, using a reference sample of real audio data recorded in the environment described above.

B. Tracking Examples
A typical example of the tracking results achieved with the SBF-IS algorithm is depicted in Figure 1.It contains the plots of the estimated source position versus time resulting from the two PF methods.The grey lines above and below the estimated source position represent plus/minus one standard deviation of the particle set for both the x and y-coordinates.The audio data used in this example was recorded in a real office room with reverberation time T 60 = 0.39s and average SNR 9.4dB.The acoustic source was moving at a constant speed along a straight line over a distance of about 1.6m.The signal recorded with one of the array sensors is given as an example in the top plot of Figure 1.This practical result also demonstrates the reinitialisation capabilities of the IS method, with the set of particles purposely initialised in a random room location at the start of the simulation, about 2m away from the true start position of the target.As soon as the source starts emitting an acoustic signal, the IS method is able to relocate its particles towards the true source position and subsequently tracks the target as it moves across the state space.The non-IS filter is unable to detect the source due to the current measurement data not being taken into account when propagating the particles.The situation described in Figure 1 typically constitutes an example of target detection (track acquisition), for which the IS method clearly shows its superiority over a pure bootstrap implementation.More results on the tracking performance of algorithm SBF-PL can be found in [3].
The results depicted in Figure 2 were obtained with a scenario where two speakers take part in an alternating conversation.The simulation was carried out using the image method to generate signals originating from two different locations in the above mentioned setting, with a reverberation time T 60 = 0.35s.White noise was added to the microphone signals with an SNR level of about 20dB.The plot at the top of Figure 2 shows an example of signal resulting for one of the sensors.The vertical dotted lines represent time instants at which a speaker change occurs in the original source signal.The plots at the bottom of Figure 2 show the tracking results obtained with the SBF-IS algorithm.This demonstrates once again the efficiency of this method which automatically switches between talkers as soon as a speech signal is detected at a different location in the state space.

C. Image Method Results
Results presented in the previous section specifically demonstrate the performance of the IS algorithm during the phase of target detection, i.e. in localisation mode.This section deals with a more specific assessment of the PF operating in tracking mode only.To this purpose, the particles were initialised at the true source location at the beginning of each simulation in the following results.
For this experiment, the microphone signals were generated with the image method [18] for varying values of reverberation time T 60 .White noise was added to the resulting signals with an approximate SNR level of 20dB.A single example of target trajectory and source signal was considered, with a path corresponding to a 1.6m straight line across the room.The source signal was a sentence uttered by a male speaker, defining a 7.3s audio sample.
The results presented in Figure 3 were obtained by simulating each PF algorithm 100 times for the considered audio data.For each run, an estimate of the tracking accuracy was computed as the average deviation (root mean squared error, RMSE) of the PF location estimate from the true source trajectory.The statistical distribution of this assessment parameter (for each value of T 60 ) is plotted in Figure 3 using a boxplot representation, which contains information about interquartile range and median of the RMSE data set.
For low to medium reverberation times, i.e. up to T 60 ≈ 0.6s, these results show that the median tracking accuracy of both IS-based and non-IS methods are similar.Simulation runs for which the PF does not recover after losing track of the target result in the appearance of a second mode in the distribution of the RMSE parameter.This effect can be seen easily in the SBF-PL results for reverberation times greater than about 0.4s, whereas the reinitialisation capabilities of the SBF-IS method allows such cases to be mostly avoided.On the other hand, algorithm SBF-IS exhibits distributions of the RMSE results that are more spread out: the outliers appear here as the tail of the distribution rather than a separate mode.This results from the SBF-IS algorithm occasionally reinitialising off-track (i.e.erroneously) and then recovering, rather than due to a complete and definitive loss of the target as with SBF-PL.

D. Further Discussion
When designing any tracking algorithm, a compromise must be found between its localisation ability and its tracking accuracy.With the proposed IS algorithm, this can be achieved very efficiently by tuning the prior probabilities of the reinitialisation and importance sampling options, P R and P S respectively.A bootstrap implementation constitutes an extreme limit in this tradeoff with P R = P S = 0.
On the basis of a (non-optimised) Matlab implementation, it can be seen that the SBF-IS algorithm requires roughly twice more computational power than SBF-PL to process the same amount of input data.This is of course due to the additional task of computing the importance function over a fixed grid of points across the state space.However, a realtime implementation of the SBF-IS algorithm, running on a 1.7GHz computer in conjunction with a 16 sensor array, shows that this additional processing power requirement does not represent any difficulties for modern desktop computers.Given this hardware setup, the number of particles in the IS algorithm can be increased up to 120 before reaching the limits of the system resources, which proves to be more than sufficient for the considered application.This practical implementation demonstrates the robustness of the IS algorithm when localising sources and tracking fast target motions in the setting of a 3.5m×4.5m×2.7moffice room with a practically measured reverberation time T 60 = 0.5s.Demonstration movies (originally recorded in real-time) showing some typical examples of the IS algorithm output delivered by this implementation can be found online at http://www.watri.org.au/¬ericl/tracking.html.
Finally, it must be kept in mind that the tracking performance of the IS method developed in this paper can be potentially largely improved by using some additional information (such as e.g.voice activity detection) to adjust the reinitialisation probability P R .The use of a more elaborate beamforming principle providing improved localisation estimates would also lead to a better tracking performance.

VII. CONCLUSION
Speaker localisation and tracking are complicated array processing applications, made especially challenging by complex reverberation effects and the discontinued nature of speech signals.Adopting a Bayesian filtering approach to this problem leads to superior tracking performance compared to traditional acoustic localisation methods.In this paper, we have developed a particle filtering technique using the principle of importance sampling.The resulting algorithm is able to automatically recover from track loss, detect a new source entering the acoustic scene, and switch between speakers taking turns, thus making it more suitable than bootstrap methods in practice.In a practical tracking system, a bootstrap-only algorithm would typically necessitate additional processing units to deal with such scenarios, whereas the IS method already integrates these functionalities at a low level in the algorithm.

Fig. 1 :
Fig. 1: Tracking results obtained with an IS-based and a non-IS method.Top plot: example of signal recorded with one array sensor for this simulation.Bottom plots: true source position (dotted line), source location estimate (solid lines) and lines representing + − one standard deviation of the particle set (grey lines).

Fig. 2 :
Fig. 2: tracking results with alternating conversation scenario.Top plot: example of audio signal generated for one of the array sensors.Vertical dotted lines denote a change of speaker.Bottom plots: tracking results in x and y-coordinates.Dotted lines represent the position of the active source.

Fig. 3 :
Fig.3: Statistical tracking performance results obtained with simulated reverberant data (image method) for various levels of reverberation.In each boxplot, the dots represent RMSE data points, the lines at the top and bottom of the box correspond to the 75th and 25th percentile of the data set, respectively, and the horizontal line in the middle of the box is the median of the data set.
{Y 1 , . . ., Y k } represents the concatenation of all measurements up to time k.The posterior PDF p(X k |Y 1:k ) contains all the statistical information available regarding the current condition of the state variable X k .An estimate X k of the state then follows, for instance, as the mean or the mode of this PDF.
•) and h(•) are possibly nonlinear functions, and u k and v k are possibly non-Gaussian noise variables.