Francois-Xavier Socheleau
email: fx.socheleau@imt-atlantique.fr

Signal Processing in Acoustics: Underwater Communications and Networking Joint signal detection and channel estimation in multi-scale multi-lag underwater acoustic environments

We consider the problem of jointly detecting a known signal and estimating the channel in the frequency range used by underwater acoustic communication systems. A multi-scale multi-lag propagation environment is considered and no prior knowledge of the channel order is assumed. The proposed detection/estimation method is based on the framework of multifamily generalized likelihood ratio tests applied to a signal lying in a union of subspaces. The result is a "tuning-free" orthogonal matching pursuit algorithm with a stopping criterion that does not require the knowledge of the number of channel taps or the noise variance. The performance is illustrated with replay simulations using real shallow-water channels measured in the Mediterranean Sea. Numerical results show that the proposed method outperforms competing algorithms in terms of both detection probability and channel estimation error. In addition, channel estimation does not exhibit a performance floor as observed with fixed-order-based approaches.

INTRODUCTION

Acoustic signal detection and channel estimation is of interest for a wide range of applications, including underwater source localization and wireless communications. For instance, model-based localization techniques seek to estimate the channel structure in order to compute source position estimates by solving an inverse problem. [START_REF] Michalopoulou | Matched-Impulse-Response Processing for Shallow-Water Localization and Geoacoustic Inversion[END_REF][START_REF] Gall | Matched-Field Processing Performance Under the Stochastic and Deterministic Signal Models[END_REF][START_REF] Gomes | Building Location Awareness into Acoustic Communication Links and Networks through Channel Delay Estimation[END_REF] Similarly, channel estimation can be very useful to recover the transmitted information in a communication context. [START_REF] Stojanovic | Channel-Estimation-Based Adaptive Equalization of Underwater Acoustic Signals[END_REF][START_REF] Li | Estimation of Rapidly Time-Varying Sparse Channels[END_REF] In both cases, estimation is often made possible by detecting and processing a known signal referred to as beacon or preamble.

Detection of a known signal usually relies on a bank of matched-filters (BMF). Depending on the underlying assumptions on the channel and the noise, detection is performed by comparing the BMF outputs to some threshold after being possibly normalized, combined and/or ordered (maximization over Doppler branches, time delay averaging, accumulation of the strongest correlations, rank reduction, etc.). [START_REF] Li | Preamble Detection for Underwater Acoustic Communications Based on Sparse Channel Identification[END_REF][START_REF] Van Walree | The Watermark Benchmark for Underwater Acoustic Modulation Schemes[END_REF] Estimation methods rely on multipath channel models with various assumptions about their structure. For instance, the model may have a predefined degree of sparsity and/or it may explicitly consider the multi-scale structure due to path-dependent Doppler shifts. [START_REF] Gomes | Building Location Awareness into Acoustic Communication Links and Networks through Channel Delay Estimation[END_REF][START_REF] Li | Estimation of Rapidly Time-Varying Sparse Channels[END_REF][START_REF] Socheleau | Non Data-Aided Estimation of Time-Varying Multiscale Doppler in Underwater Acoustic Channels[END_REF][START_REF] Zhao | Parameter Estimation of Wideband Underwater Acoustic Multipath Channels Based on Fractional Fourier Transform[END_REF] Be it for the detection of a known signal or for channel estimation, the channel order is often assumed to be known or fixed a priori. This knowledge is either explicit by considering a fixed number of dominant taps or implicit by assuming some knowledge on the signal-to-noise ratio or by setting specific values for sparse regularization parameters. When the actual order departs from the presumed one, the performance degrades significantly. If it is overestimated, noise may be injected in the detection and estimation procedures. If underestimated, only part of the available signal energy is exploited.

The main contribution of this work is to propose a single method for joint signal detection and channel estimation with a multi-scale multi-lag (MSML) model and without prior knowledge of the channel order. This contribution is based on the framework of multifamily generalized likelihood ratio tests (MFGLRT) applied to a signal lying in a union of subspaces (UoS). [START_REF] Socheleau | A Multifamily GLRT for CFAR Detection of Signals in a Union of Subspaces[END_REF] In this framework, the channel output signal is assumed to lie in one out of a possible set of known subspaces. The candidate set of subspaces is here built from all the possible MSML channel structures. Although efficient, [START_REF] Socheleau | A Multifamily GLRT for CFAR Detection of Signals in a Union of Subspaces[END_REF] the resulting detection/estimation method is too complex to be applicable in our context. Therefore, we propose an approximate solution based on a modified implementation of orthogonal matching pursuit (OMP). This solution automatically estimates the signal subspace as well as the channel order from the data without the knowledge of the noise variance. The benefit of the proposed method is illustrated with replay simulations using real shallow-water channels measured in the Mediterranean Sea. Numerical results show that it outperforms competing algorithms in terms of both detection probability and channel estimation error.

The paper is organized as follows. Sec. 2 is devoted to the problem formulation. The derivation of the MFGLRT as well as its approximated formulation based on OMP is presented in Sec. 3. Numerical results are provided in Sec. 4, followed by conclusions in Sec. 5.

PROBLEM FORMULATION

We consider the possible transmission of a known signal x(t) over a MSML underwater acoustic communication channel. At reception, the baseband signal y(t) is either modeled as

H 0 : y(t) = w(t), (1)
when no signal is transmitted or as

H 1 : y(t) = L-1 � �=0 α � (t)x(t -τ � (t))e -i2πfcτ � (t) + w(t), (2)
when x(t) is transmitted. α � (t) denotes the complex attenuation of the �-th channel path, f c is the carrier frequency and τ � (t) is the delay assumed to be a linear function of time described by an initial delay τ 0 � and a drift � � (also called the Doppler scale or the time-compression factor), so that τ � (t) = τ 0 � + � � t. w(t) is the additive noise assumed to be Gaussian of unknown variance σ 2 . Based on the observation y(t) and the knowledge of x(t), the problem addressed in this paper is first to decide between the two hypotheses H 0 and H 1 , and then, when H 1 is decided, to estimate the unknown channel parameters

� L, � α � (t), τ 0 � , � � � {�=0,••• ,L-1}
� . The amplitudes of the channel paths are assumed constant over the duration of x(t), i.e α � (t) ≈ α � , so that the discrete-time baseband signal can be expressed, under H 1 , as

y = z + w with z = Xα, (3)
where y and

w ∈ C N , α ∈ C L and X ∈ C N ×L is a matrix defined as X = [x 0 , • • • , xL-1 ].
N denotes the number of signal samples and the k-th entry of each vector x� ∈ C N satisfies

x� (k) = x � kT s (1 -� � ) -τ 0 � � e -i2πfc(� � kTs+τ 0 � ) , (4)
where T s is the sampling period.

JOINT DETECTION AND ESTIMATION METHOD

Eq. ( 3) shows that signal z lies in a subspace spanned by the columns of X. In our context, the actual channel parameters are unknown so that X is also unknown. However, in practice, it is often reasonable to bound some of the channel parameters with known values. More specifically, 0 < L ≤ L max , -� max ≤ � � ≤ � max and 0 < τ 0 � ≤ τ max . In addition, since any digital estimator has a finite resolution, it is commonly assumed that each pair � τ 0 � , � � � approximately lies on a predefined delay-Doppler grid. Although

X is unknown, we can therefore consider that all possible values taken by X are known in advance such that X ∈ { X1 , • • • , XM }, where each Xi is known and where M denotes the number of all possible sets {L,

� α � (t), τ 0 � , � � � {�=0,••• ,L-1} }.
Note that M can be extremely large for practical values of L max , � max , τ max so that the prior knowledge of each Xi may be unusable because of computational limitations. However, for the sake of clarity, we first consider in Sec. 3.1 a fictitious scenario in which computational complexity is not an issue and present a joint detector/estimator that makes use of the knowledge of each Xi in an MFGLRT framework. An approximate solution with a much lower complexity is then proposed in Sec. 3.2.

A. MFGLRT

When all possible values taken by X are known in advance, z is said to lie in a union of M subspaces, i.e., z ∈ ∪ M i=1 S i . In this case, z belongs to one of the subspaces S i , spanned by the columns of Xi , but we do not know a priori to which one. The detection of an unknown signal that lies in a union of subspaces and that is observed in additive Gaussian noise with unknown variance has recently been addressed. [START_REF] Socheleau | A Multifamily GLRT for CFAR Detection of Signals in a Union of Subspaces[END_REF] It is shown that an efficient detector can be obtained by using the framework of multifamily generalized likelihood ratio tests. This detector is expressed as the following test

max 1≤i≤M g i (L i (y)) H 1 ≷ H 0 η, (5)
where

L i (y) = 1 2 ln � y H P S i y y H (I N -P S i ) y × N -n i n i � , (6)
and η is the detection threshold chosen so as to satisfy some predefined probability of false-alarm. P S i denotes the projection matrix in subspace S i defined as

P S i = Xi � XH i Xi � -1 XH
i . We assume each matrix Xi to have at most N -1 linearly independent columns so that n i = dim (S i ) < N . L i (y) is a GLR statistic and g i penalizes this statistic to counteract its tendency to increase with the model order n i . The analytic expression of g i is given in [10, Sec. III-B] and is shown to depend on the Legendre transform of the cumulant generating function (CGF) of L i (y) under hypothesis H 0 . g i is termed as the MFGLR statistic in the sequel. Not only Eq. ( 3) is a signal detector but it also provides an estimate of the subspace where z actually lies, that is

Ŝ = S i * with i * = argmax 1≤i≤M g i (L i (y)) . (7)
In other words, Eq. ( 7) provides an estimate of the channel parameters L and

� τ 0 � , � � � {�=0,••• ,L-1} and the estimate of X is � X = Xi * .
The amplitudes of the channel paths can then be obtained by the following least-square estimate:

α = � � XH � X� -1 � XH y. (8)

B. APPROXIMATED MFGLRT (AMFGLRT)

As mentioned previously, in most practical scenarios, M is extremely large so that the MFGLRT presented in Eqs. ( 5) and ( 7) is infeasible. For instance, if N grid designates the number of possible delay-Doppler pairs for each path, then

M = � Lmax �=1 � N grid �
� . However, the complexity can be drastically reduced by decoupling (at least in part) the estimation of subspace S from the computation of g i (L i (y)). Let Xgrid denote the matrix of N grid columns containing all possible time and Doppler shifted versions of the transmitted signal x. Using this definition, z can be seen as a sparse signal in the basis Xgrid . Therefore, an estimate of subspace S can be obtained by solving the following least-square problem

Ŝ = span{ Xgrid (U * )}, (9)
with

U * = argmin U : |U |≤Lmax min α � � �y -Xgrid (U ) α � � � 2 . (10)
U is any combination of the set {0, 1, • • • , L max -1} and |U | denotes the cardinality of U . Xgrid (U ) is the submatrix of Xgrid with columns indexed in U . U * is therefore the estimated sparse support set. Problem (10) can be solved using any sparse signal recovery algorithm such as matching or basis pursuit. [START_REF] Elad | Sparse and Redundant Representations: from Theory to Applications in Signal and Image processing[END_REF] However, for optimal performance, these algorithms require some prior knowledge on the sparsity degree, i.e., the actual channel order L, and/or on the noise power to efficiently tune their parameters. As illustrated in Sec. 4, estimation performance is highly sensitive to this tuning. The main contribution of this work is to combine the resolution of Eq. ( 10) using OMP with the computation of g i (L i (y)) as a decision statistic as well as a stopping criterion for OMP. More specifically, OMP is a greedy algorithm developed for sparse support set estimation. As recalled in Alg. 1, it iterates until some stopping condition is satisfied. In our context, since we do not know the actual value for |U |, we cannot set a priori the number of OMP iterations. Similarly, since the noise variance is unknown, we cannot set a predefined threshold on the energy of the so-called residual to stop the decomposition. Our approach relies on the property that MFGLR approaches are able to provide model order estimates by maximizing the test statistic g i . [START_REF] Kay | Exponentially Embedded Families-New Approaches to Model Order Estimation[END_REF][START_REF] Kay | The Multifamily Likelihood Ratio Test for Multiple Signal Model Detection[END_REF] In Problem (10), the model order is |U |. We use this property in OMP to stop the procedure. At each iteration i, the MFGLR statistic obtained with subspace Ŝi Algorithm 1: OMP OMP is stopped as soon as Ŝi leads to a MFGLR statistic smaller than the one obtained at the previous iteration with Ŝi-1 , meaning that the model order has been found. The great benefit of this approach is that no prior knowledge of the channel order (or sparsity degree) is required. OMP becomes "tuning-free".

Input: y, Xgrid , Output: Ŝ, � X, α 1 Initialization: i = 1, U 0 = ∅, r 0 = y 2 Support set estimation : U i = U i-1 ∪ s i with s i = argmax u=1,••• ,N grid � � �r H i-1 Xgrid (u) � � � 3 Subspace estimation : Ŝi = span{ Xgrid (U i )} 4 Residual update : r i = (I N -P Ŝi )y 5 if stopping condition is satisfied then 6 Ŝ = Ŝi , � X = Xgrid (U i ), α = � � XH � X� -
In addition, the same statistic is used for stopping the support set estimation procedure and for deciding whether a signal is present or absent. The full method is detailed in Alg. 2. As mentioned previously, the MFGLR statistic depends on the Legendre transform of the CGF of L i (y) under hypothesis H 0 . Although a closed-form expression of this statistic is available in the context of union of subspaces, this not the case anymore since the sparse support set is now estimated with OMP so that it becomes random and datadependent. However, an approximation can be found using a second-order expansion of the CGF of L i (y).

The resulting expression is given in step 5 of Alg. 2, where µ i = E(L i (y)) and

σ i = E(L 2 i (y)) -µ 2 i
under hypothesis H 0 . These values, as well as the detection threshold η, are obtained through Monte-Carlo simulations with noise only.

Algorithm 2: Approximated MFGLRT (AMFGLRT)

Input: y, Xgrid , µ i , σ 2 i , η Output: Ŝ, � X, α, Detection decision 1 Initialization: i = 1, U i = ∅, r 0 = y, g0 (y) = -∞ 2 Support set estimation : U i = U i-1 ∪ s i with s i = argmax u=1,••• ,N grid � � �r H i-1 Xgrid (u) � � � 3 Subspace estimation : Ŝi = span{ Xgrid (U i )} 4 Residual update : r i = (I N -P Ŝi )y 5 MFGLR statistic: gi (y) = 2 � λ i L i (y) -µ i λ i + σ 2 i 2 λ i � � {L i (y)-µ i } with L i (y) = 1 2 ln   y H P Ŝi y y H � I N -P Ŝi � y × N -i i   , λ i = L i (y)-µ i σ 2 i . � {•} denotes the indicator function. 6 if gi (y) < gi-1 (y) then 7 if gi-1 (y) > η then 8 Detection decision = true 9 Ŝ = Ŝi-1 , � X = Xgrid (U i-1 ), α = � � XH � X� -1 � XH y

NUMERICAL RESULTS

The performance of AMFGLRT is illustrated with replay simulations using real shallow-water channels measured in the Mediterranean Sea. Replay of time-varying impulse responses (TVIRs) distorts input signals by convolving them with measured channels. [START_REF] Van Walree | The Watermark Benchmark for Underwater Acoustic Modulation Schemes[END_REF][START_REF] Otnes | Validation of Replay-Based Underwater Acoustic Communication Channel Simulation[END_REF][START_REF] Socheleau | Stochastic Replay of non-WSSUS Underwater Acoustic Communication Channels Recorded at Sea[END_REF][START_REF] Socheleau | Parametric Replay-Based Simulation of Underwater Acoustic Communication Channels[END_REF][START_REF] Socheleau | Stochastic Replay of SIMO Underwater Acoustic Communication Channels[END_REF] Table 1 summarizes measurement conditions and test channel parameters. TVIRs are shown in Fig. 1. Drifts of multipath arrivals are clearly visible. This drifts result from the transmitter and receiver movements but also from the possible deviations of instrument clock frequencies from their nominal values (such deviations have the same effect as that of genuine Doppler shifts).

For both channels, the detection and estimation performance are measured with a BPSK-modulated signal x(t), whose simulation parameters are detailed in Table 2. The additive noise is white Gaussian, i.e., w ∼ CN (0, σ 2 w I N ), and the signal-to-noise ratio (SNR) is defined as the energy-per-bit over the noise power spectral density, that is: The delay-Doppler grid used to build matrix Xgrid has the following parameters: maximum Doppler � max = 5/1500 ≈ 3.3 10 -3 , Doppler step � step = 0.1/1500 ≈ 6.7 10 -5 , maximum delay spread τ max = 25 ms and delay step τ step = T s . Consequently, N grid = 10638 for the simulation with MED1 and 3650 for the one with MED2.

SNR = �z� 2 2 L seq σ 2 w . (11)

A. DETECTION PERFORMANCE

The performance of AMFGLRT is compared with two other detectors. The first one is the bank of constant-false-alarm-rate (CFAR) matched filters (BMF) made of Doppler-shifted signal replicas. This detector can be expressed as

max 1≤i≤N grid y H P Ŝi y y H � I N -P Ŝi � y H 1 ≷ H 0 η MFB with Ŝi = span{ Xgrid (i)}. (12)
Such a detector based its decision by only considering the main tap of the channel. The second detector is the one corresponding to the use of the GLRT statistic when the signal subspace is estimated with a fixed number n it of OMP iterations. After simplification, it satisfies

y H P Ŝn it y y H � I N -P Ŝn it � y H 1 ≷ H 0 η OMPn it ( 13
)
where Ŝn it is the output of Alg. 1 when the stopping condition in step 5 is set to i = n it . Setting the number of iterations to a fixed value n it is equivalent to assuming that the number of channel taps equals this value. We investigate the performance of the proposed detector in terms of probability of detection P D . The probability of false alarm P F A is set to 10 -3 . The detection thresholds and the probabilities of detection are obtained with 100/P F A and 2000 independent trials, respectively. For each trial, the signal x(t) is injected into the channel with a random delay. Fig. 2 shows the probability of detection as function of the SNR for AMFGLRT, BMF and OMP. Similar behaviors are observed in both channels: AMFGLRT outperforms the other methods and the performance of OMP degrades significantly as n it increases. In addition, performance in MED1 is better than in MED2, owing to a higher SNR gain provided by a longer m-sequence. The strength of AMFGLRT is that the number of channel taps (or equivalently, the dimension of subspace S) is estimated on the data and is not fixed a priori. For low SNR values and when the sparse support set is unknown, detectors better have to only consider the main taps and ignore taps with low amplitudes. This is because the actual delay-Doppler coordinates of taps with low amplitudes are very unlikely to be correctly estimated for low SNRs and may result in noise injection in the test statistic. However, as the SNR increases, it becomes beneficial to consider additional taps. The adaptation capacity of AMFGLRT is illustrated in Fig. 3 where the average number of estimated taps is shown as a function of the SNR for both channels. Note that the average number of taps is greater for MED1 than for MED2, which is expected since the delay resolution is inversely proportional to the signal bandwidth. As discussed in the next section, this adaptability feature is also very useful to optimize the channel estimation performance.

L, � α � (t), τ 0 � , � � � {�=0,••• ,L-1}
� are unknown. Therefore, we use the normalized mean-square error (NMSE) of the signal reconstruction as a performance metric. It is defined as

NMSE = E    � � �z - � X α� � � 2 2 �z� 2 2    , (14)
and is estimated with 2000 Monte-Carlo trials. The performance of AMFGLRT is compared with the same methods used in the previous section. For BMF, we set � X = Xgrid (U BMF ), where the sparse support set U BMF is estimated as the union of matched-filter branch indexes having a CFAR statistic greater than some threshold controlled by a false-alarm probability set to P FA = 10 -3 . Fig. 4 shows the NMSE as a function of SNR for all methods. Once again, AMFGLRT outperforms the other approaches in both channels. OMP offers a reasonable performance in the low SNR region with a small number of iterations but very quickly exhibits a performance floor due to the underestimation of the number of channel taps. This performance floor is shifted to large SNR values when n it increases but at the cost of degrading the performance for small SNRs. BMF also exhibits a performance floor whose position may change depending on P FA . In any case, it does not provide satisfactory performance. Channel-order adaptation, as provided by AMFGLRT, is clearly beneficial for all SNR ranges. Note that other estimation approaches, such as basis pursuit, have also been tested. The results are not shown here because their behavior is similar to OMP in that it does not adapt to the SNR and shows performance flooring depending on the value set for the regularization parameters.

CONCLUSION

We have proposed a joint signal detector and channel estimator in the framework of a sparse representation problem combined with a multifamily generalized likelihood ratio test. A key feature of this method, called AMFGLRT, is its capacity to estimate the channel order from the data. For low SNRs, AMFGLRT only considers the main channel taps and ignore taps with low amplitudes. This limits the amount of noise injected in the detection and estimation procedures. As the SNR increases, the additional signal energy provided by these weaker taps becomes beneficial and is taken into account by AMFGLRT. Using real channel data, this adaptability has been shown to provide a significant performance improvement, especially for channel estimation.

Although designed with the assumption of additive Gaussian noise, AMFGLRT is expected to be robust to impulsive noises with little adaptation. First, normalized CFAR statistics are known to be more robust to local non-stationarity than statistics assuming the knowledge of the noise power. Moreover, before being processed, observed data could be filtered by non-linear functions such as those used in robust estimation. [START_REF] Huber | Robust Statistics[END_REF] This analysis is left for future work.

15

 15 Set i = i + 1 and go to step 2

Figure 1 :

 1 Figure 1: Time-varying impulse responses: (a) MED1, (b) MED2

Figure 2 :Figure 3 :

 23 Figure 2: Probability of detection vs signal-to-noise ratio, P FA = 10 -3 . (a) MED1, (b) MED2

Figure 4 :

 4 Figure 4: Estimation performance vs signal-to-noise ratio. (a) MED1, (b) MED2

Table 1 :

 1 Channel parameters and sounding conditions

	Name	MED1	MED2
	Environment	Bay of Toulon	Bay of Toulon
	Time of year	July	July
	Range Water depth Transmitter depth	∼700 m ∼50 m 10 m (towed)	∼700 m ∼50 m 10 m (towed)
	Receiver depth	4 m (suspended) 12 m (suspended)
	Probe signal type	pseudo noise	pseudo noise
	-3dB freq. band	5.7-15.3 kHz	6.85-10.15 kHz
	Delay coverage	128 ms	42 ms
	Doppler coverage	7.8 Hz	23.8 Hz
	Play time	25 s	25 s

Table 2 :

 2 Simulation parameters for x(t)

	Name	MED1	MED2
	Modulation	BPSK m-sequence BPSK m-sequence
	Length of the sequence (L seq )	511	255
	Duration Symbol rate	∼60.1 ms 8.5 kbds	∼87.3 ms 2.92 kbds
	Roll-off (RRC filter)	0.1	0.1
	Sampling rate (T s )	0.235 ms	0.685 ms

ACKNOWLEDGMENTS

The author would like to thank the French Defense Procurement Agency (DGA) for providing the data used in this work.