A Framework for Multiple Object Tracking in Underwater Acoustic MIMO Communication Channels

This work presents a computational framework for the analysis and design of large-scale algorithms utilized in the estimation of acoustic, doubly-dispersive, randomly time-variant, underwater communication channels. Channel estimation results are used, in turn, in the proposed framework for the development of efficient high performance algorithms, based on fast Fourier transformations, for the search, detection, estimation and tracking (SDET) of underwater moving objects through acoustic wavefront signal analysis techniques associated with real-time electronic surveillance and acoustic monitoring (eSAM) operations. Particular importance is given in this work to the estimation of the range and speed of deep underwater moving objects modeled as point targets. The work demonstrates how to use Kronecker products signal algebra (KSA), a branch of finite-dimensional tensor signal algebra, as a mathematical language for the formulation of novel variants of parallel orthogonal matching pursuit (POMP) algorithms, as well as a programming aid for mapping these algorithms to large-scale computational structures, using a modified Kuck’s paradigm for parallel computation.


Introduction
This work formulates a computational framework for the development of efficient algorithms to effect computational signal processing operations to address the problem of electronic surveillance and acoustic monitoring (eSAM) of deep underwater moving objects.From the point of view of computational complexity theory, the combined problem of search, detection, estimation and tracking (SDET) of deep underwater moving objects through active signal sensing is considered an NP-hard problem [1][2][3].This is in the context of transmitter signal waveform design and signal-to-noise ratio (SNR) optimization.For this reason, large-scale modeling and simulation environments are needed in order to analyze, design and implement high performance computational signal processing algorithms that may provide optimal approximate solutions to this problem.In order to facilitate the exposition of the proposed computational framework, we proceed to simplify the formulation of the SDET problem in the following manner.
Problem formulation: There is a need to develop approximate algorithms for the estimation of the range and speed of deep underwater moving objects to facilitate tracking operations [4,5].These tracking operations are performed in a communication medium that manifests time-frequency dispersion effects and must be characterized as an acoustic, linear, stochastic (ALS), time-variant channel [6].
To obtain approximate solutions to this simplified SDET problem, we decomposed the problem itself into three main stages: (i) how to characterize the communication channel between the set of signal transmitters, assumed here to be M transmitters, and the set of scattering moving objects; (ii) how to characterize the scattering interaction among the transmitted signals and the scattering moving objects, assumed here to be K L scattering objects; and (iii) how to characterize the communication channel between the set of scattering moving objects and the set of signal receivers, assumed here to be N signal transducers?
Stages (i) and (iii) were integrated in a unified manner, and the communication channels were modeled as ALS time-variant systems resulting in a MIMO channel characterization formulation.Stage (ii) was addressed using multidimensional, multi-component, polynomial phase (MMPP) signals and operator theoretic properties associated with discrete Cohen distributions (DCDs) under a bidimensional harmonic analysis setting.This approach at solving the SDET problem allowed us to integrate the channel estimation problem of digital communications and the scattering object parameter estimation problem of sensor array signal processing.
This paper presents algorithmic formulations for the solution of the simplified SDET problem using Kronecker products signal algebra (KSA) as the mathematical language to convey our proposed computational signal processing methods in parallel format, which, to the best of our knowledge, are original formulations.Thus, this work extends the work of Li, W.and Preisig, J.C. on the estimation of rapidly time-varying sparse channels by providing a solution to the estimation problem for the multiple-input multiple-output (MIMO) case in parallel form, using Kronecker products signal algebra [7].The language of Kronecker products has been demonstrated to be a very useful tool for the mathematical formulation of fast unitary transforms.In 1990, Rodríguez, D. co-authored a 52-page tutorial article on a methodology for designing, modifying and implementing Fourier transform algorithms on various architectures using tensor or Kronecker products and an associated class of stride group permutations for numerical data flow management [8].
The results presented here are based on two additional works: first, on the work of Carrascosa, P.C. and Stojanovic, M. on underwater acoustic channel estimation, where they introduce spatial multiplexing to increase the data rate supported by the band-limited restriction of the channel [9,10]; second, on the work of Yatawatta, S. and Petropulu, A. on blind channel estimation, where they allow the users at the transmitter to transmit simultaneously, without any bandwidth restrictions, each user utilizing a single antenna [11].This work introduces innovation along three main venues: (i) using a mathematical language formulation to describe the processes of channel estimation and multiple object tracking in a unified manner; (ii) developing unique cyclic discrete time-frequency distributions for range-Doppler estimation; and (iii) exploiting the algebraic properties of block circulant structures in order to reduce the overall computational complexity of channel estimation algorithms.We provide in Table 1 a description of the salient attributes of the channel estimation algorithms of these selected works, as well as the salient features of multiple object tracking (MOT) MIMO estimation algorithms presented in this work.
In 1993, Van Loan, C. wrote a now well-known book on computational frameworks for fast Fourier transform algorithms.The book explains how "different FFT algorithms correspond to different factorizations of the discrete Fourier transform (DFT) matrix," and it shows how computer implementations of FFT algorithms are directly related to their mathematical formulations in Kronecker products notation [12].Van Loan's book extended, in a systematic manner, the approach presented in [8] for identifying Kronecker products expressions with computational constructs.In 2006, Franchetti, F. et al., presented a software program generator and optimizer for linear signal transforms, such as the DFT [13].This language generator was written using the Kronecker products notation formulated in [8].In 2007, Ali, A. et al. developed a portable framework for FFT algorithms to run on various parallel architectures.The computational framework was also formulated using the language of Kronecker products [14].In 2008, Rodríguez, D. co-authored an article where a methodology was presented for the high-level partitioning of signal transforms onto distributed hardware architectures using, again, the language of Kronecker products signal algebra [15].A direct mapping from Kronecker products expressions, which we termed functional primitives, to MATLAB (MATrix LABoratory) pseudocode, was clearly established by Franchetti, F. et al., in 2009, in an article that provides "the techniques needed to implement the discrete Fourier transform (DFT) efficiently on current multicore systems", and also discusses optimization techniques for good multicore performance [16].In 2009, Kepner, J. wrote a book for parallel programming in MATLAB, with specific examples of application algorithms written using the open source pMatlab software library [17].A good feature of this book is how it presents "a more sophisticate model [to] handle hierarchical collections of processors, memories, and networks".
The rest of this document is structured as follows.Section 2 introduces the concept of delay-Doppler MIMO estimation, describes our original MIMO channel model and presents three matching pursuit algorithm variants for its impulse response function estimation.Section 3 describes Kuck's paradigm for parallel computation used in our work as the computational model for the implementation of the Kronecker products-based algorithms.Section 4 presents a description of the implementation results obtained from our computational testbed for a particular solution of the simplified SDET problem when we simulate the detection and tracking of two deep underwater objects, at rest or moving at a constant speed.Two acoustic signal transmitters and two acoustic signal receivers are used to actively sense the communication medium.The section also describes an example of deep underwater tracking simulation of ten (10) objects, four (4) of them permanently fixed.Finally, Section 5 provides our conclusions.

Delay-Doppler MIMO Channel Characterization
In this section, we formulate the characterization of MIMO acoustic, linear, stochastic (ALS), time-variant systems and describe a set of matching pursuit greedy algorithms to estimate MIMO channel parameters.The characterizing channel function used in this work is the delay-Doppler spread function U(ξ, ν) [18].This function plays the role of surrogate function for the input-delay spread function; that is, for the impulse response function h(t, ξ) of acoustic linear stochastic (ALS) channels [19].
The parameters to be estimated are time-delays and the Doppler shifts associated with the delay-Doppler spread function.This estimation work must be carried out in a fraction of the time corresponding to the channel's coherence time, T C ; that is, the time period during which the impulse response function is assumed to be time-translation invariant.This approach for channel parameter estimation is done under a multiple-input multiple-output (MIMO) system assumption.Thus, we consider M transmitter transducers and N receiver transducers in a MIMO system scenario.It is possible in this scenario that a determined number of scatterers, L, may be present between the transmit and receive transducers.We assume that some scatterers may be in motion, moving at a constant speed.These movements introduce Doppler effects, which are manifested as frequency shifts acting over each copy of the transmitted signals z m (t) associated with each scattering point P l , l ∈ Z L , in motion.
An abstract depiction of the bistatic sonar nature of our proposed deep underwater multiple object tracking scenario is given in Figure 1.Under the MIMO assumption, we can express each received signal w n (t) as the sum of L copies of each transmitted signal z m (t), time delayed by ξ l , shifted by its respective Doppler frequency ν l and scaled by an attenuation factor α l .We define the following expression: Thus, we obtain the following result for each received signal: where t ∈ R and w, z ∈ L 2 (R).Here, α l,m,n ∈ C is the attenuation factor α l associated with each input signal z m (t) and received signal w n (t).The parameters ξ l,m,n , ν l,m,n ∈ R are the time delay and Doppler shift associated with each scatterer l ∈ Z L , transmitter m ∈ Z M and receiver n ∈ Z N .The parameters ξ 0,m,n and ν 0,m,n are set to zero, for line of sight (L.o.S.) consideration.
The signal η(t) is a real-valued, independent, wide-sense, stationary, Gaussian stochastic process that represents the channel noise.This approach captures the more important aspects of a physical realization of a deep underwater ALS MIMO channel.Under the MIMO assumption, we can obtain the impulse response function of the deep underwater ALS channel by substituting each input signal z m (t) with the impulse function δ m (t) in Equation (2).We proceed to define the following expression: ( Thus, we obtain: where t ∈ R and w, z ∈ L 2 (R).Here, h n,m (t, ξ) (kernel function) is the time-variant impulse response between the transmitter m ∈ Z M and the receiver n ∈ Z N .
As previously stated, the MIMO channel model formulated in this work has M transmitter transducers and N receiver transducers.Equation ( 5) presents the input-output relation of our MIMO channel model.
From Figure 2, we may obtain, after a series of algebraic manipulations, using concepts from our proposed computational framework, the following channel input/output relationship: where is the Kronecker products operator, is the vertical matrix concatenation operator, is the horizontal matrix concatenation operator, Z ∈ l 2 (Z KN × Z KMN ) is the input matrix, h ∈ l 2 (Z KMN ) is the impulse response channel vector, w ∈ l 2 (Z KN ) is the output vector and C is the circulant operator.We have taken advantage of the coherent condition of the channel and have invoked the properties of the convolution operation.This MIMO system model formulation allows the study of filter bank multi-carrier communication, partial FFT demodulation and waveform diversity design [20][21][22][23].While in the channel sounding mode, a series of predefined pulses z i , i ∈ Z M is transmitted [24].These pulses are compared with the received signals w j , j ∈ Z N for characterizing the MIMO channel using the inverse relation shown in Equation (6). where is the channel impulse response vector and w ∈ l 2 (Z KN ) is the signal output vector [25].It is important to point out that to arrive at the mathematical formulation for the MIMO channel input/output system configuration, although not quite difficult, is an intricate and laborious procedure.We illustrate this procedure by describing in some detail the mathematical formulations for the SISO, MISO and SIMO channel input/output configurations that we first derived in order to arrive at our proposed MIMO channel model system configuration [26].

SISO, MISO and SIMO Channel Configurations
To arrive at our MIMO channel model system configuration, we start by presenting the SISO, MISO and SIMO configurations.The idea behind this approach is to visualize the algebraic structures in each case.Each configuration has a particular matrix-vector formulation.These formulations allow us to develop data structures and procedures appropriate to address the efficient channel estimation problem.Another important objective in this section is to demonstrate the importance of Kronecker products signal algebra in matrix structure manipulation.Regularities in the matrix structures can be exploited for improving channel estimation algorithms.This work exploits the matrix-vector formulations and sparsity properties; however, many other regularities and symmetries can be identified and reviewed to develop new algorithms.

SISO Configuration
The SISO configuration has just one transmitter transducer and one receiver transducer.Equation (7) shows its input-output relationship.
where z ∈ l 2 (Z K ) is the transmitted signal, h ∈ l 2 (Z K ) is the impulse response function of the channel, w ∈ l 2 (Z K ) is the received signal, K is the cyclic convolution operator of order K and C is the circulant operator.Equation (7) shows a typical input-output relationship for a linear and time-invariant system, also known as a discrete-time filter.Figure 3 illustrates an SISO channel model system.The matrix-vector representation of the configuration is presented in Equation (8), which shows the first relevant matrix structure of our MIMO channel model configuration, known as a circulant matrix: . . .The computational complexity, under the direct computation approach, associated with this configuration is O(K 2 ), and it is given by the complexity of the inverse operator acting over the matrix Z of dimension K × K and the complexity of the matrix-vector product of the vector w, of length K. Therefore, the dominant complexity is O(K 2 ), where K is the number of samples in the impulse response h[k], k ∈ Z K .

MISO Configuration
The MISO configuration has M transmitter transducers and just one receiver transducer.Equation ( 9) shows its input-output relationship. where is the received signal, K is the cyclic convolution operator of order K, is the horizontal matrix concatenation operator, is the vertical matrix concatenation operator, Z ∈ l 2 (Z K × Z KM ) is the input channel matrix, h ∈ l 2 (Z KM ) is the impulse response channel vector and C is the circulant operator.Figure 4 illustrates the MISO channel model system.A matrix-vector representation for this model is shown in Equation (10).
The computational complexity, under direct computation approach, associated with the MISO channel estimation problem is O((KM) 3  3 ), and it is given by the complexity of the pseudo-inverse operator acting over the matrix Z of dimension K × KM and the complexity of the matrix-vector product of the vector h, of length KM.Therefore, the dominant complexity is O((KM) 3 ), where K is the number of samples in the impulse response h[k], k ∈ Z K , and M is the number of transmitter transducers.

SIMO Configuration
The SIMO configuration has just one transmitter transducer and N receiver transducers.Equation (11) shows its input-output relationship. where is the i-th received signal, K is the cyclic convolution operator of order K, is the vertical matrix concatenation operator, Z ∈ l 2 (Z KN × Z KN ) is the input channel matrix, h ∈ l 2 (Z KN ) is the impulse response channel vector, w ∈ l 2 (Z KN ) is the output vector and C is the circulant operator.Figure 5 illustrates the SIMO model system.A matrix-vector representation is shown in Equation (12).
. . .The computational complexity, under the direct computation approach, associated with the SIMO problem, is bounded by O((KN) 3 + 2(KN) 3 ) = O((KN) 3 + 2(KN) 3 ), and it is given by the complexity of the pseudo-inverse operator acting over the matrix Z, of order KN × KN, and the complexity of the matrix-vector product of the vector h, of length KN.Therefore, the dominant complexity is O((KN) 3 ), where K is the number of samples in the impulse response function h[k], k ∈ Z K , and N is the number of receiver transducers.
Three matching pursuit algorithm variants were implemented in this work: basic, orthogonal and order-recursive least-square matching pursuit.Basic matching pursuit (BMP) is the most simple implementation of these algorithms.BMP chooses the maximum projection (inner product) of the column set of Z on the output vector w and extracts the corresponding contribution of each maximum on w.The orthogonal matching pursuit variant adds a new stage to the estimation process.After each iteration, it recalculates all of the λ i coefficients associated with each chosen c i column.This is to force the orthogonal condition on the w residual vector on the subspace spanned by c i selected columns.
The process of selecting the columns is similar to the process used in the basic matching pursuit variant.In the order-recursive least-square matching pursuit algorithm variant, the method of selecting the columns is improved in the following manner.Each column is chosen such that it minimizes the residual value of the w vector.Therefore, each column is chosen considering the previous set of chosen columns for minimizing the w residual vector, in order to calculate the λ i coefficients, so that it is similar to the former algorithm variant.These coefficients provide a guide for the stopping criteria, as pertaining to residual errors, the channel impulse response formulation and the overall performance of the algorithms.

Kuck's Paradigm of Computation
David Kuck presented an abstract manner to represent parallel computational structures [27].This methodology of representation may be related to the parallel random machine model (PRAM) used to describe an abstract parallel machine, even though our formulation is a more sophisticated model.The goal of a parallelization effort is to reduce the overall computation time by distributing a computational task among available processors.
As its name indicates, the PRAM was intended as the parallel-computing analogy to the random-access machine (RAM).Like random access machine (RAM) models that are used by sequential algorithm designers to model algorithmic performance, the PRAM is used by parallel algorithm designers to model parallel algorithmic performance (such as time complexity, where the number of processors is a parameter).PRAM used the concept of shared memory to avoid the communication problem in message passing interface (MPI) models.The complexity is expressed as a function of the dimension of input set N and the number of processors Ncpus.
Kuck's approach creates a hierarchical structure to represent in a parallel manner the execution of an algorithm.This hierarchical structure reflects the connection between the computing units and shared memory stages.Each level is tagged with a specific label, like "1.5".This representation approach is sufficiently abstract to allow implementations in a wide spectrum of modern architectures.Figure 6 illustrates Kuck's diagram representation of the parallel MIMO ALS estimation problem [27].In this figure, Level 0 represents independent memory units, M 0 and uniprocessors P 0 , which contain registers and caches for fast access to frequently-used data.At level 0.5, these units are interconnected with a network N 0.5 , generally a bus, mesh, hypercube, shuffle or other mechanism of data interchange.This network provides communication between processors, but does not provide shared memory.Therefore, a protocol for message passing is necessary.The first level of shared-address memory space is SM 1 .The communication between the processors in Level 0 and SM 1 is managed by SMN 1 .There exists a great difference between direct access to memory via SMN 1 and indirect access to memory via N 0.5 .The difference is appreciated in the poor performance produced by the indirect access to memory.Now, we can apply a recursive principle to generate new levels.
The algorithmic formulations for solving the SDET problem were implemented in the pMatlab and Parallel Computing Toolbox (PCT) programming environments for MATLAB.These MATLAB parallel programming environments allowed us to conduct detailed computational complexity analysis on the mapping of Kronecker product-based algorithmic formulations to modeled distributed and shared memory structures (see Algorithm 1 and Table 2).when conducting the simulation on four or more processing units.We attribute this experimental observation to the following facts: (i) each MIMO channel link is difficult to parallelize by itself; and (ii) the communication time between processing units increases significantly beyond two units in a pMatlab parallel programming environment.The scattering interaction among the designed waveforms utilized at the transmitters was addressed using multidimensional, multi-component, polynomial phase (MMPP) signals and the processing of these signals using discrete Cohen distributions (DCDs) under a bidimensional harmonic analysis setting.In particular, the efficient computation of cross-ambiguity function operations was instrumental in solving the multiple objects scattering problem.Figures 10-23 show the result of our bidimensional harmonic analysis algorithm applied to two nearby scatterers when these scatterers get closer.The bidimensional harmonic analysis was conducted using a diverse set of designed waveforms, including, sinusoidal, square and linear modulation or "chirp" pulses.The results showed, as expected, that the bidimensional harmonic analysis using chirp waveforms outperforms the analysis conducted using square waveform pulses.The point target location accuracy of the SDET computational framework environment is addressed in this work on two fronts, i.e., the accuracy of the cross-ambiguity function of the transmitted and received pulses based on the design of the transmitted waveform and the accuracy of the relative position between the transmitter and receiver transducers in a MIMO configuration topology.To address the resolution issue through transmitter waveform generation for the cross ambiguity product, we followed the work of C. Y. Chen and P. P. Vaidyanathan on linear frequency modulation (LFM) chirp pulse time-bandwidth product design [28], as well as Brookner, E. work on MIMO array configuration design [29].The relative position between transmitter and receiver transducers must be known to a sub-wavelength precision in order to improve a given point target location accuracy, and we followed the work of Pailhas, Y. and Petillot, Y. to accomplish this task [30].

Deep Underwater Multiple Object Tracking Example
In this section, we describe a generalized testbed algorithm framework for simulating the tracking of multiple objects moving at a constant speed in a deep underwater environment.Figure 24 presents a tracking allocation chart when the number of tracked objects is equal to ten (10).However, the testbed algorithm framework was tested under multiple object tracking scenarios, each tracking scenario producing a different tracking allocation chart.Figures 25-34 show the performance of the ten-object tracking simulation example.For simplicity, the chirp pulse waveforms were normalized in the amplitude and were also time normalized by setting the sampling time equal to unity.

Conclusions
This paper presented a computational framework for the analysis and design of signal processing algorithms utilized in the estimation of acoustic, doubly-dispersive, randomly time-variant, deep underwater communication channels.It introduced innovation along three main venues, i.e., using a mathematical language formulation to describe channel estimation processes, developing unique cyclic discrete time-frequency distributions for range Doppler estimation and exploiting the algebraic properties of block circulant structures to reduce the overall computational complexity.The results presented here are based on the work of Li, W. and Preisig, J.C. on the estimation of rapidly time-varying sparse channels, the work of Carrascosa, P.C. and Stojanovic, M. on underwater acoustic channel estimation using spatial multiplexing and the work of Yatawatta, S. and Petropulu, A. on blind channel estimation.Multiple object tracking operations were successfully accomplished using a 2 × 2 MIMO channel estimation testbed structure.
The work presented in this paper demonstrated that the language of Kronecker products remains, after more than twenty (20) years, an invaluable tool for the formulation of signal processing algorithms.Kronecker products signal algebra, a branch of finite dimensional multilinear algebra, continues to produce fruitful results when dealing with the treatment of high dimensional signals.
For example, we took advantage of the fact that the discrete Fourier transform (DFT) matrix may be represented as a product of sparse matrices, both using decimation in time (DIT) and decimation in frequency (DIF) to formulate cyclic filtering operations, using the property of the DFT matrix for the diagonalization of circulant matrices.This procedure allowed us to use the commutation theorem of Kronecker products in order to reduce the communication time between processors when dealing with large-scale block circulants.To the best of our knowledge, the parallel formulations of the matching pursuit algorithms presented in Kronecker products signal algebra language in this paper are new.These algorithm formulations were readily implemented using the pMatlab and PCT parallel programming environments.Through this parallel algorithm implementations, we demonstrated that the ORLSMP formulation was the most accurate; however, it was also the slowest.

Table 1 .
Attributes of selected channel estimation algorithms.MOT, multiple object tracking.

Table 3 .
Time performance of MIMO matching pursuit algorithms.BMP, basic matching pursuit.