A Novel Approach for the Estimation of Doubly Spread Acoustic Channels Based on Wavelet Transform

: In this paper, the estimation of doubly spread underwater acoustic (UWA) channels is investigated. The UWA channels are characterized by severe delay spread and signiﬁcant Doppler effects, and can be well modeled as a multi-scale multi-lag (MSML) channel. Furthermore, exploiting the sparsity of UWA channels, MSML channel estimation can be transformed into the estimation of parameter sets (amplitude, Doppler scale factor, time delay). Based on this, the orthogonal matching pursuit (OMP) algorithm has been widely used. However, the estimation accuracy of OMP depends on the size of the dictionary, which is related with both delay spread and Doppler spread. Thus it requires high computational complexity. This paper proposes a new method, called wavelet transform (WT) based algorithm, for the UWA channel estimation. Different from OMP algorithm which needs to search in both time domain and Doppler domain, WT-based algorithm only needs to search in time domain by using the Doppler invariant characteristic of hyperbolic frequency modulation (HFM) signal. The performance of the proposed algorithm is evaluated by computer simulations based on BELLHOP. The simulation results show that WT-based algorithm performs slightly better than OMP algorithm in low signal to noise ratio (SNR) while can greatly reduce computational complexity


Introduction
Underwater acoustic (UWA) channels pose grand challenges for reliable high data-rate communications, due to significant Doppler effects and severe delay spread [1][2][3].In underwater communications, acoustic waves propagate at 1500 m/s, much lower than the speed of electromagnetic waves in terrestrial wireless systems (which is 3 × 10 8 m/s).Thus, the motion of platform will cause more significant Doppler effects, which is expressed as signal compressing or dilating in time domain and can be treated as Doppler scale [4].In addition, the delay spread results in severe inter-symbol interference (ISI) [5].To fully understand the channel characteristics and overcome challenges it poses, accurate channel models and estimation methods are essential to investigate.
As observed in many experiments [6,7], signals from different paths will experience different Doppler scales, arrive at different time and have different energies.Therefore, the multi-scale multi-lag (MSML) channel, denoted in [8], can well model acoustic channel and has been adopted in many researches [7,9].In MSML channel model, each path can be parameterized by Doppler scale factor, time delay and amplitude.However, for severe delay spread, this model will be too complex to estimate.
To overcome such difficulties, many researchers investigated the sparsity of UWA channel, that is, most of the energy is concentrated in some small regions.Hence only a few channel taps in MSML channel model are nonzero and need to be tracked.As a result, the computational complexity has been reduced and many sparse channel estimation algorithms based on compressed-sensing (CS) have been proposed [7,[10][11][12][13][14][15].
These algorithms can be generally grouped in two categories: dynamic programming like matching pursuit (MP), and linear programming like basis pursuit (BP) [7].BP aims at the l 1 minimization which needs high computational complexity, thus is less attractive for practical large-scale applications.Therefore, we will mainly focus on MP algorithm and its successors.
In [10], matching pursuit (MP) algorithm is applied to estimate Doppler scale factors of different paths.It iteratively selects one column from the dictionary that is most relevant with the residual signal, and subtracts the estimated path component to update the residual signal.Compared with MP, its orthogonal version, orthogonal matching pursuit (OMP) algorithm [11] makes the residual signal be orthogonal with all the selected columns, thus has better convergence speed and accuracy.The traditional subspace methods and CS-based methods are compared in [7], and the author concludes that CS-based methods have better performances in channel estimation.Meanwhile, some improved algorithms, which focus on adaptively estimating path numbers, like sparsity adaptive matching pursuit (SaMP) [12] and adaptive step size SaMP (AS-SaMP) [13] have been applied to sparse channel estimation.Furthermore, some references propose methods to reduce computation: a two-stage OMP algorithm is proposed by [14], which estimates the Doppler scale factor and time delay respectively, rather than simultaneously as OMP does.However, it requires some preprocessings before channel estimation.The author of [15] proposes that fast Fourier transform (FFT) can be utilized in OMP to simplify calculation.However, the reduction is limited as it only focuses on the computing process rather than on the reducing of column dimensions in the dictionary.Therefore, the main limitation of MP algorithm and its successors is that the estimation accuracy depends on the size of the dictionary, which is determined by the product of grid point numbers on the time delay dimension and the Doppler scale dimension.To guarantee a fine resolution, the number of columns in the dictionary could be extremely large.Thus calculating inner products of the received signal and each column lead to extensive calculation, especially for UWA channel, where the delay-scale spread is large.To overcome this difficulty, this paper proposes a novel wavelet transform (WT) based algorithm by employing the hyperbolic frequency modulation (HFM) signal as training signal.
In UWA environment, Doppler effects cause frequency shifts of the transmitted signal.HFM signal is a preferable training signal in UWA communication as the frequency shifts caused by Doppler effects can be treated as a time shift [16,17].Thus the Doppler effects and time delay can be simultaneously estimated by analyzing the time-frequency characteristics of the signal.In addition, we choose WT in this paper as it is an effective measurable tool for time-frequency analysis and has been used for the signal analysis in engineering, biomedical and researcher application [18,19].WT methods consist of two categories: continuous wavelet transform (CWT) and discrete wavelet transform (DWT), DWT is nothing but discretization of CWT.
The main contributions of this paper are as follows: • We propose to use a superimposed HFM signal as training signal and prove that the HFM up sweep and down sweep signals will not interfere each other for WT.

•
We propose an iterative algorithm to estimate MSML channel based on WT.

•
We use computer simulations based on BELLHOP to investigate the performance of the proposed algorithm, and make comparisons with OMP algorithm.

•
We present the complexity analysis of the proposed algorithm as well as OMP algorithm.
The rest of this paper is organized as follows: Section 2 gives the MSML channel model.In Section 3, brief introductions of HFM signal and WT are included and we present the WT-based algorithm.Section 4 focuses on the simulation results and complexity analysis.Section 5 concludes the paper.
Notation: we will use the following notations in this paper: Upper (lower) bold-face letters stand for matrices (vectors); Superscript * denotes conjugate.

Channel Model
UWA channel features large multipath spread for the interaction with ocean surface, bottom medium, as well as inhomogeneous particles of water column.MSML channel model can well describe the UWA multipath channel, that is: where L is the number of channel taps.A l (t) is the time-varying amplitude of the lth path and can be assumed to be constant during a short period of time, for example, the transmission duration of one data frame [7].τ l and a l are the time delay and Doppler scale factor of the lth path, respectively.In addition, δ(t) is a delta function defined as follows: Let s(t) be the transmitted signal and the corresponding received signal r(t) can be written as where n(t) is the additive noise.
Given the sparsity of UWA channel, only some channel taps are nonzero in (1), which means that r(t) is a superposition of only a few scale-delayed versions of s(t).Therefore, the calculation complexity for channel estimation is significantly reduced.

HFM Signal
HFM signal is a kind of non-linear frequency modulation signal whose instantaneous frequency is a hyperbolic function of time, which can be expressed as [16]: where In addition, T is the duration of the HFM signal, f 1 is the starting frequency and f 2 is the end frequency.If f 1 < f 2 then s(t) is HFM up sweep signal, otherwise represents HFM down sweep signal.The signal instantaneous frequency is the derivation of the phase term of the cosine function: We consider to apply HFM signal as the channel training signal for its good Doppler-invariability.Introducing a Doppler scale factor a to (4), the phase term ϕ(t) becomes Correspondingly, the instantaneous frequency is In addition, from [17], there is a constant time shift ∆t = (a−1)t 0 a to satisfy the following equation.
Thus, the HFM signal is Doppler invariant.

CWT
As mentioned before, WT is an effective measurable tool to analyze non-stationary signal, and the instantaneous frequency of HFM varies with time.Therefore, WT is a preferable choice for the time-frequency analysis of HFM signal.The CWT of the received signal r(t) can be expressed below: where ψ(t) is the basis function that called mother wavelet [20], and a W and b are the time scale and time shift, respectively.Selection of proper mother wavelet is very important for effectively time-frequency analysis.There are different types of mother wavelet, such as Haar, Daubachies, Sine and Gaussian.In addition, they all satisfy the following admissible condition [21]: where ψ(w) is the Fourier transform of ψ(t).
It is easy to verify that a signal with finite length in both time domain and frequency domain will satisfy (10), Agrawal, S. et al. [22] points out that the WT efficiency will increase if characteristic of mother wavelet is closer to the input signal.So we consider to employ the transmitted HFM signal as mother wavelet, and since the transmitted signal has limited samples, thus it's length is limited and can satisfy (10).

WT-Based Algorithm
Firstly, we consider the transmitted HFM signal as s(t), then the received signal is as in (3), plugging (3) in ( 9) and let ψ(t) = s(t), we get CWT displays the distribution of signal amplitude and phase into time scale a W and time shift b.Thus, with proper a W and b, which match to the signal from the lth path, WT coefficients will appear a peak.Therefore, we can get the frequency information of r l (t) at time shift b.With this, we can set the values of a W and b to be a list of possible values and get the WT coefficients matrix.Thus get the time-frequency information of r l (t).To get better time-frequency analysis, the values of a W and b should be as many as possible.The two-dimensional searching requires large computation cost.
Considering the good Doppler-invariability of HFM signal as we mentioned before, we can simplify the two-dimensional searching to one-dimension, that is, by using the transmitted HFM signal as the mother wave, the time scale a W can be set as 1 and the WT coefficients will meet a peak at a time shift for each path as the following: The detail information can be seen in Appendix A. From (12), the time shift is determined by both time delay and Doppler scale factor.So it is difficult to estimate a l and τ l simultaneously when only sending one HFM signal.We consider to send a superimposed signal, which is defined as The variables in ( 13) are defined as in ( 4) and f 1 < f 2 .Therefore, the first part in ( 13) is a HFM up sweep signal and can be defined as s + (t); the second part is a HFM down sweep signal and can be denoted as s − (t).
The received signal will be where n(t) is the additive noise.In addition, at the receiver, WT will be performed twice and the mother wave will be s + (t) and s − (t), respectively.We have verified that when use s + (t) as mother wave, a constant time shift does not exist which will make s + (t) match to s − (a l (t − τ l )) in frequency in Appendix B. Thus the superposition of s + (t) and s − (t) will not interfere each other, and corresponding time shifts for the lth path are: The Doppler scale factor and time delay can be calculated by In addition, the path amplitude can be calculated by normalizing (11), that is For MSML acoustic channels, we propose an iterative algorithm to estimate parameters for each path.Specifically, at each iteration, path with the highest energy can be identified and will be separated from the received signal.In addition, parameters {A l , a l , τ l } can be estimated at the same time.The whole estimation scheme is given as follows.
WT-based iterative channel estimation algorithm: Input: Transmitted HFM up sweep signal vector S + , and HFM down sweep signal vector S − ; received signal vector r; a small positive ε.
where l represents the lth path, and the path amplitude Âl can be calculated as (19).( 4) Update residual signal: where N is the length of r e .( 5) If E − ||r e || 2 2 /||r|| 2 2 < ε, stop the iteration, else, set i = i + 1 and go to step 1).Output: The estimated channel parameters { Âi , âi , τi } L i=1 .Note: since the number of paths L is unknown, we use ε as a threshold to determine the termination of the algorithm.ε represents the normalized decrease of residual signal, and when all paths have been detected, the decrease will be small and ε can be set as a small positive number, such as 0.01 in this paper.

Performances
In this part, we use computer simulations to evaluate the performance of the proposed WT-based algorithm, and comparisons with OMP algorithm will also be included.
The parameters of the transmitted superimposed HFM signal are set as follows: the starting frequency and end frequency are f 1 = 6 kHz, and f 2 = 10 kHz.The time duration is set to be T = 200 ms, and the sampling rate is f s = 40 kHz.
The computer simulation is based on BELLHOP and is implemented in Matlab software.BELLHOP is a beam tracing model for predicting acoustic pressure fields in ocean environments [23].The underwater environment is set as: the water is 100 m deep, the initial horizontal range between the transmitter and receiver is 1000 m. the transmitter is fixed at the depth of 80 m, and the receiver is at 30 m depth with a horizontal speed of 4.5 m/s toward the transmitter.The speed of the acoustic wave is set to be 1500 m/s.In addition, the reflection coefficients of the bottom and the surface are set to be 0.7 and −0.9, respectively.As shown in Figure 1.At the receiver, both WT-based algorithm and OMP algorithm will be applied to channel estimation.For WT-based algorithm, WT will be applied twice to search for b l + and b l − at each iteration, while for OMP algorithm, we build a dictionary with a resolution of 1 × 10 −3 in the Doppler scale factor and 1/ f s in the tap delay.Figures 2 and 3 display the estimating errors of the Doppler scale factor and time delay versus signal to noise ratio (SNR), as defined in the following: Error  From Figures 2 and 3, we can see that WT-based algorithm is slightly better than OMP algorithm when SNR is lower than 12 dB.In addition, at high SNR, both algorithm can achieve low Doppler scale estimation errors, and the delay estimation errors are nearly 0, which implies that WT-based algorithm is an effective algorithm for channel estimation as the classical OMP algorithm while only needs one-dimension searching.
Figure 4 illustrates the residual signal energy rate, ||r e || 2 /||r|| 2 , versus SNR of both estimation algorithms.In addition, a reference which uses the true channel information is also included.It shows that the performance of WT-based algorithm approaches to the lower bound while OMP algorithm has some performance losses at low SNR.Residual signal energy rate is a comprehensive evaluation of the channel estimation, that is, the estimation accuracies of Doppler scale, delay and amplitude all contribute to it.Thus the performance in Figure 4 is also coincident with which we analyzed in Figures 2 and 3   In addition, according to the reviewer's suggestions, we also include residual signal energy rate versus SNR at different sampling rates: f s = 20 kHz in Figure 5 and f s = 60 kHz in Figure 6.It shows that WT-based algorithm can maintain good performances at different sampling rates, while OMP algorithm shows performance degradation at low sampling rate.

Complexity Analysis
Both WT-based and OMP algorithms estimate path parameters iteratively, and for proper selection of ε, the iteration number is nearly the true channel path number for both of them.Therefore, the mainly differences of computation lie in the process of one iteration.Hence we will focus on the computation of one iteration.
For OMP algorithm, parameters in both time domain and Doppler domain need to be searched.Let K L be the length of the training sequence, N τ be the grid points on the time dimension, N a be the grid points on the Doppler scale dimension, thus the column numbers of the dictionary is N = N τ N a .Thus, the inner products of the received signal and columns in the dictionary require ρ = N τ N a K L complex multiplications.
For WT-based algorithm, WT will be performed twice at each iteration.As mentioned before, HFM signal is Doppler invariant, thus it only needs to search in the time domain.let K L be the length of the training sequence, and N τ be the grid points on the time dimension, and then the computation Therefore, WT-based algorithm only needs one-dimension searching and the computational complexity is much lower than OMP algorithm, especially when the grid points on the Doppler scale dimension is large.

Conclusions
We have investigated the problem of doubly spread UWA channels estimation in this paper and proposed the WT-based algorithm.By using the superimposed HFM signal as training signal, the proposed algorithm only needs to search in time domain, and the Doppler scale factor and time delay can be simultaneously estimated.Therefore, it greatly reduce the computational complexity than OMP algorithm, while has slightly better performance.
where a l is the Doppler scale factor, a l > 1 when the transmitter and the receiver are approaching each other; a l ∈ (0, 1) when they are away from each other; and a l = 1 when they are static.Therefore a l = 0 and b varies with t.The assumption that b is constant is incorrect.In addition, it is the same for s + (t).

Initialization:
Set the residual signal vector r e = r; iteration index i = 1; Iterate: (1) Calculate and save the energy rate of r e as: E = ||r e || 2 2 /||r|| 2 2 .(2) Use S + , S − as mother waves and do WT with r e , respectively.Select the maximum peaks and get corresponding time shifts b + max , b − max .(3) Calculate the path parameters: âl

Figure 1 .
Figure 1.Acoustic ray paths of UWA channel.

Figure 2 .
Figure 2. Errors of the estimating Doppler scale factor versus SNR.

Figure 3 .
Figure 3. Errors of the estimating delay versus SNR. .

Figure 5 .
Figure 5. Residual signal energy rate versus SNR with sampling rate at 20 kHz.

Figure 6 .
Figure 6.Residual signal energy rate versus SNR with sampling rate at 60 kHz.