# Fast raypath separation based on low-rank matrix approximation in a shallow-water waveguide

Longyu Jiang, Wenbo Song, Zhe Zhang, Chunfeng Yang, Shijie Wang,

Philippe Roux

► To cite this version:

Longyu Jiang, Wenbo Song, Zhe Zhang, Chunfeng Yang, Shijie Wang, et al. Fast raypath separation based on low-rank matrix approximation in a shallow-water waveguide. Journal of the Acoustical Society of America, 2018, 143 (4), pp.EL271-EL277. 10.1121/1.5030916. hal-01807117

HAL Id: hal-01807117

https://univ-rennes.hal.science/hal-01807117v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fast raypath separation based on low-rank matrix approximation in a shallow-water waveguide

Longyu Jiang; Wenbo Song; Zhe Zhang; Chunfeng Yang; Shijie Wang; Philippe Roux

<figure><img src="image_2.png" /></figure>

Check for updates

J. Acoust. Soc. Am. 143, EL271–EL277 (2018)

https://doi.org/10.1121/1.5030916

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

## Articles You May Be Interested In

Active wideband higher-order raypath separation in multipath environment

J. Acoust. Soc. Am. (January 2017)

Two-dimensional active raypath separation using examination of the roots of the spectrum polynomial

J. Acoust. Soc. Am. (October 2017)

Three-dimensional higher-order raypath separation in a shallow-water waveguide

JASA Express Lett. (July 2022)

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>

# Fast raypath separation based on low-rank matrix approximation in a shallow-water waveguide

Longyu Jiang,<fnref n="1" />,ᵃ) Wenbo Song,<fnref n="1" /> Zhe Zhang,<fnref n="1" /> Chunfeng Yang,<fnref n="1" />,ᵇ) Shijie Wang,<fnref n="1" /> and Philippe Roux<fnref n="2" />

**Abstract:** Subspace algorithms based on higher-order cumulants were developed to achieve high-resolution separation in non-Gaussian processes. However, singular value decomposition (SVD) of a huge matrix is an unavoidable step of these algorithms. The memory space and running time required by the decomposition are super-linear with respect to the size of the matrix, which is prohibitive in terms of practical applications. Thus, in this paper, a fast raypath separation algorithm based on low-rank matrix approximation is proposed in a shallow-water waveguide. The experimental results illustrate that the proposed algorithm dramatically reduces the consumption of time and space, with arbitrarily small error, compared to conventional higher-order cumulant-based algorithms.

© 2018 Acoustical Society of America
[DHC]

Date Received: November 7, 2017

Date Accepted: March 21, 2018

## 1. Introduction

In shallow water, acoustic waves propagate in multi-paths due to refraction and reflection at the ocean surface and the continental shelf. Multi-path propagation has been widely used in many applications, such as active sonar, ocean acoustic tomography,<fnref n="1" /> and source localization, because more raypaths provide more information. However, the successful use of the information is based on the premise of first eliminating the raypaths' interferences with each other produced by the multi-path propagation. Among the existing algorithms to separate raypaths, subspace-based algorithms are noticeable for their large resolution improvement. In particular, the multiple signal classification (MUSIC) algorithm,<fnref n="2" /> as a representative subspace-based method, was applied to separate raypaths under the assumption of an uncorrelated narrow-band signal. However, MUSICA fails to separate raypaths when they are fully correlated or coherent. To overcome this drawback, Jiang et al.³,⁴ developed a high-resolution method called the smoothing-multiple signal classification active large band (MUSICAL) algorithm, which greatly improves the separation performance, especially when two fully correlated raypaths arrive close together. Although the smoothing-MUSICAL algorithm was developed under the assumption of white Gaussian noise, colored noise exists in real ocean environments. Recently, Jiang et al.⁵ presented a higher-order wideband active algorithm to separate raypaths interrupted by colored noise. Naturally, it takes advantage of fourth-order cumulants containing non-Gaussian information in non-Gaussian processes. The experimental results obtained from tank data at small scale⁵ illustrate that the higher-order algorithm can obtain more accurate results and suppress noise; however, it has large computational cost (both time and memory space) for computing the singular value decomposition of a huge cumulants matrix. This cost is prohibitive from a practical perspective. To overcome this drawback, we propose in this study a fast higher-order separation algorithm, which has the advantages of low computational cost and high resolution. In addition, its performance will be tested using synthetic data, experimental data obtained in a tank and ocean data. The most important contribution is that the proposed algorithm provides a strategy to find the signal subspace efficiently, which can be extended to all

<footnote n="1">¹The Laboratory of Image Science and Technology, Southeast University, Nanjing 210096, China</footnote>

<footnote n="2">²Institut des Sciences de la Terre, Universite Joseph Fourier, Centre National de la Recherche Scientifique, 1381 Rue de la Piscine, Saint-Martin d'Heres, France
JLY@seu.edu.cn, 220161523@seu.edu.cn, 220151575@seu.edu.cn,
chunfeng.yang@seu.edu.cn, mike@seu.edu.cn, philippe.roux@ujf-grenoble.fr</footnote>

<footnote>ᵃ)Also at: Key Laboratory of Computer Network and Information Integration, Ministry of Education, Southeast University, Nanjing 210096, China. Author to whom correspondence should be addressed.</footnote>

<footnote>ᵇ)Also at: Centre de Recherche en Information Biomedicale Sino-Francais (LIA CRIBs), Rennes, France.</footnote>
the subspace-based algorithms. Based on the strategy, their computational cost can be largely reduced especially when the cumulants matrix has large dimensions.

This paper is organized as follows. Our fast raypath separation algorithm is introduced in Sec. 2. In Secs. 3 and 4 we discuss the performance and the error of the proposed method. In Sec. 5, its performance test using synthetic data, experimental data obtained in a tank, and ocean data is conducted, respectively. We draw conclusions in Sec. 6.

## 2. A fast higher-order separation algorithm

We consider the case where P raypaths arrive on a vertical array composed of M sensors. The received signal of the mth ($m = 1, 2, ..., M$) sensor in the frequency domain is modeled by

$$x_m(\nu) = \sum_{p=1}^{P} a_p e(\nu) \exp(-j2\pi\nu(T_p + (m-1)t_m(\theta_p))) + n_m(\nu), \quad (1)$$

where the term $e(\nu)$ is the deterministic amplitude of the emitted signal at frequency $\nu$. The amplitude of each raypath, $a_p$, is considered random, where $T_p$ represents the arrival time of the pth raypath on the reference sensor and $t_m(\theta_p)$ is the delay between the reference sensor and the mth sensor. Additionally, $t_m(\theta_p) = d \sin(\theta_p)/c$, where $\theta_p$ is the arrival direction of raypath p on the array, d is the distance between two adjacent sensors, c is the propagation velocity of the acoustic signal, and $n_m(\nu)$ is the additive noise received at the mth sensor.

We use F frequency bins of the wideband signal. Hence, Eq. (1) is written in matrix form as follows:

$$\mathbf{x}_g = \mathbf{H} \cdot \mathbf{A} + \mathbf{n}_g, \quad (2)$$

where $\mathbf{x}_g = [\mathbf{x}^+(\nu_1), \mathbf{x}^+(\nu_2), ..., \mathbf{x}^+(\nu_F)]^+$ with $\mathbf{x}(\nu_i) = [x_1(\nu_i), x_2(\nu_i), ..., x_M(\nu_i)]^+$ ($i = 1, ..., F$); that is, $\mathbf{x}_g$ is a long-vector of size $M \times F$, which is arranged according to all the frequency bins on all the sensors. $\mathbf{n}_g$ is similarly the concatenation of the observed noise vectors at each frequency. Specifically, $\mathbf{n}_g = [\mathbf{n}^+(\nu_1), \mathbf{n}^+(\nu_2), ..., \mathbf{n}^+(\nu_F)]^+$ with $\mathbf{n}(\nu_i) = [n_1(\nu_i), n_2(\nu_i), ..., n_M(\nu_i)]^+$. $\mathbf{A} = [a_1, a_2, ..., a_P]^+$ is a vector of dimension P. $\mathbf{H} = [\mathbf{h}_1, \mathbf{h}_2, ..., \mathbf{h}_P]^+$ with $\mathbf{h}_p = [e(\nu_i)e^{-2j\pi\nu_i\tau_{lp}}, ..., e(\nu_F)e^{-2j\pi\nu_F\tau_{Mp}}]^+$, where $e^{-2j\pi\nu_i\tau_{lp}}$ ($i = 1, ..., F; m = 1, ..., M$) is the transfer function between the source and the mth sensor, and + indicates transposition.

We first apply spatial-frequency smoothing⁴,⁵ as a preprocessing step to eliminate the interferences among raypaths. That is, for a spatial smoothing factor $K_s$, an antenna composed of M sensors is divided into $2K_s + 1$ sub-antennas. The $k_s$th sub-antenna is composed of the sensors in the range of [$k_s, k_s + 1, ..., M - 2K_s$]. At the same time, a similar frequency-smoothing process is also applied, where $K_f$ denotes the frequency-smoothing factor. Finally, using the recurrences produced by the spatial-frequency smoothing, we can estimate the wideband trispectral matrix with the following formula:

$$\hat{\mathbf{C}} = E\{(\mathbf{x}_{g,k_s,k_f} \otimes \mathbf{x}_{g,k_s,k_f}^*)(\mathbf{x}_{g,k_s,k_f} \otimes \mathbf{x}_{g,k_s,k_f}^*)^H\} \\ - E\{(\mathbf{x}_{g,k_s,k_f} \otimes \mathbf{x}_{g,k_s,k_f}^*)\}E\{(\mathbf{x}_{g,k_s,k_f} \otimes \mathbf{x}_{g,k_s,k_f}^*)^H\} \\ - E\{(\mathbf{x}_{g,k_s,k_f} \otimes \mathbf{x}_{g,k_s,k_f}^H)\} \otimes E\{(\mathbf{x}_{g,k_s,k_f} \otimes \mathbf{x}_{g,k_s,k_f}^H)^*\}, \quad (3)$$

where $\otimes$ represents the Kronecker product, $*$ indicates conjugation, and $H$ indicates the conjugate transpose. Based on the preprocessing, the rank of $\hat{\mathbf{C}}$ is equal to $K = (2K_s + 1)(2K_f + 1)$ (it is necessary to select $K$ greater than $P$).

The traditional subspace-based algorithms use methods such as the singular value decomposition (SVD) to span a signal subspace, however, they require memory and time which are super-linear in the size of the trispectral matrix (or spectral matrix). It is prohibitive for a practical application, especially when the sizes are very large. Inspired by Frieze et al.,⁶⁻⁸ we introduce a randomized strategy to span the signal subspace at low computation cost.

The strategy is based on the following three principles.

**Principle 1:** We can approximate the product of two matrices by an algorithm called the basic matrix multiplication.⁷ That is, when we input two matrices, $\hat{\mathbf{C}} \in \mathbb{R}^{(MF)^2 \times (MF)^2}$ and $\hat{\mathbf{C}}^T$, a probability distribution $\{q_i\}_{i=1}^{(MF)^2}$ and a number $c \le (MF)^2$, two output matrices $\hat{\mathbf{C}}_1$ and $\hat{\mathbf{C}}_1^T$ are returned, such that $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1^T \approx \hat{\mathbf{C}} \hat{\mathbf{C}}^T$, where $\hat{\mathbf{C}}_1 \in \mathbb{R}^{(MF)^2 \times c}$ is a matrix whose columns are c randomly chosen columns of $\hat{\mathbf{C}}$ with suitably rescale.
Principle 2: An appropriate choice of the probability distribution is crucial steps of the proposed algorithm. Although a uniform distribution could be always used, the optimal probabilities {$q_i$}$_{i=1}^{MF^2}$ are proved and defined as the following equation:⁷ Suppose $\hat{\mathbf{C}} \in \mathbb{R}^{(MF)^2 \times (MF)^2}$, $c \in \mathbb{Z}^+$, such that $1 \le c \le n$, and {$q_i$}$_{i=1}^{(MF)^2}$ are such that $q_i \ge 0$, $\sum_{i=1}^n q_i = 1$, and such that for some positive constant $\beta \le 1$: $q_i \ge \beta |\hat{\mathbf{C}}|^2 / ||\hat{\mathbf{C}}||_F^2$.

Principle 3: According to linear algebra, for symmetric matrices the left and right singular vector are the same. The singular values of $\hat{\mathbf{C}}$ are the nonnegative square roots of the eigenvalues of $\hat{\mathbf{C}}^T \hat{\mathbf{C}}$ and of $\hat{\mathbf{C}} \hat{\mathbf{C}}^T$; in addition, the left singular vectors of $\hat{\mathbf{C}}$ are eigenvectors of $\hat{\mathbf{C}} \hat{\mathbf{C}}^T$ and the right singular vectors are eigenvectors of $\hat{\mathbf{C}}^T \hat{\mathbf{C}}$.

The general process of the proposed algorithm is first described: (1) based on principle 3, finding the eigenvectors of the matrix $\hat{\mathbf{C}}$ is equivalent to find the eigenvector of the matrix $\hat{\mathbf{C}} \hat{\mathbf{C}}^T$; (2) based on principles 1 and 2, the matrix $\hat{\mathbf{C}} \hat{\mathbf{C}}^T$ can be approximated by $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1^T$ with the optimal probabilities, where $\hat{\mathbf{C}}_1$ consists of those c columns of $\hat{\mathbf{C}}$ (after appropriate rescaling); it has been proven that by choosing appropriate sampling probabilities and enough columns the approximation error can be made arbitrarily small; (3) based on principle 3, finding the eigenvectors of the matrix $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1^T$ can be transformed to find the eigenvectors of the matrix $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1$, then from $\hat{\mathbf{C}}_1^T \hat{\mathbf{C}}_1$ approximations to the top singular values and corresponding singular vectors of $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1^T$ may be computed, which are the approximations of the eigenvectors of the matrix $\hat{\mathbf{C}} \hat{\mathbf{C}}^T$. Thus, the fast algorithm are finally operated in four specific steps.

(1) To sample c columns of $\hat{\mathbf{C}}$ ($\hat{\mathbf{C}} \in \mathbb{R}^{(MF)^2 \times (MF)^2}$) randomly according to a judicious probability q, we first compute $q_i = |\hat{\mathbf{C}}^{(i)}|^2 / ||\hat{\mathbf{C}}||_F^2$, $i = 1, 2, ..., (MF)^2$ for all $(MF)^2$ columns (where $\hat{\mathbf{C}}^{(i)}$ is the ith column of matrix $\hat{\mathbf{C}}$ and $|| \cdot ||_F$ indicates Frobenius norm). Subsequently, we get a random number r, where $0 < r < 1$. If $\sum_{j=1}^{j=i-1} q_j < r < \sum_{j=1}^{j=i} q_j$, we retain the ith column. The sampling step will be repeated another c - 1 times and finally we will get c columns of $\hat{\mathbf{C}}$.

(2) According to the rule provided by the first step, we sample c columns where $1 < P < c < (MF)^2$, rescale each column by the factor $1/\sqrt{cq_i}$, $t = 1, 2, ..., c$, and arrange all c columns to construct the approximation matrix of $\hat{\mathbf{C}}$, which is denoted as $\hat{\mathbf{C}}_1 \in \mathbb{R}^{(MF)^2 \times c}$. We would like to approximate the eigenvectors of square matrix $\hat{\mathbf{C}}$. Normally, we can obtain these vectors directly through performing a SVD on matrix $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1^T \in \mathbb{R}^{(MF)^2 \times (MF)^2}$, but its size is still the same as that of $\hat{\mathbf{C}}$. Thus, we first perform a SVD on $\hat{\mathbf{C}}_1^T \hat{\mathbf{C}}_1 \in \mathbb{R}^{c \times c}$, whose size is quite smaller, to find its eigenvectors, which are also the right singular vectors of $\hat{\mathbf{C}}_1$; then, they can be transformed into the left singular vectors of $\hat{\mathbf{C}}_1$, which are also the eigenvectors of the matrix $\hat{\mathbf{C}}_1 \hat{\mathbf{C}}_1^T$. Therefore in the next steps first we carry out the SVD on matrix $\hat{\mathbf{C}}_1^T \hat{\mathbf{C}}_1$ to obtain right singular vectors of $\hat{\mathbf{C}}_1$; Second we compute $P^2$ left singular vectors of $\hat{\mathbf{C}}_1$, which successively approximate $P^2$ eigenvectors of square matrix $\hat{\mathbf{C}}$.

(3) SVD is performed on matrix $\hat{\mathbf{C}}_1^T \hat{\mathbf{C}}_1$, and the $P^2$ eigenvector $\mathbf{y}^t(t = 1, 2, ..., P^2)$, corresponding to the $P^2$ largest eigenvalues $\sigma_t(t = 1, 2, ..., P^2)$, is identified.

(4) Eigenvectors $\mathbf{h}^t(t = 1, 2, ..., P^2)$ are calculated using the equation $\mathbf{h}^t = \hat{\mathbf{C}}_1 \mathbf{y}^t / \sigma_t(\hat{\mathbf{C}}_1)$.

Based on the eigenvectors computed previously, the signal subspace $\hat{\mathbf{C}}_s$ is projected as follows:

$$ \hat{\mathbf{C}}_s = \mathbf{H}_{P^2} \times \mathbf{H}_{P^2}^*, \quad (4) $$

where $\mathbf{H}_{P^2}$ is composed of $\mathbf{h}^t(t = 1, 2, ..., P^2)$.

Finally, the direction of arrival (DOA) and the arrival time of each raypath are obtained by maximizing the following estimator:

$$ P(\theta, T) = \mathbf{a}(\theta, T)^* \hat{\mathbf{C}}_s \mathbf{a}(\theta, T), \quad (5) $$

where the wideband steering vector $\mathbf{a}(\theta, T)$ is the Kronecker product of $\mathbf{b}(\theta, T)$,

$$ \mathbf{a}(\theta, T) = \mathbf{b}(\theta, T) \otimes \mathbf{b}(\theta, T), \quad (6) $$
where $\mathbf{b}(\theta, T)$ is the concatenation of the classical steering vectors $\mathbf{d}(\nu_i, \theta)$, that is,

$$ 
\begin{aligned}
\mathbf{b}(\theta, T) &= \begin{bmatrix} e(\nu_1)e^{-2j\pi\nu_1 T} \mathbf{d}^{+}(\nu_1, \theta) \\ \vdots \\ e(\nu_F)e^{-2j\pi\nu_F T} \mathbf{d}^{+}(\nu_F, \theta) \end{bmatrix} \\
\mathbf{d}(\nu_i, \theta) &= [1, e^{-2j\pi\nu_i\tau_{1,2}(\theta)}, \dots, e^{-2j\pi\nu_i\tau_{1,M-1}(\theta)}]^{+}, \quad i = 1, 2, \dots, F.
\end{aligned} 
\quad (7) $$

### 3. Error analysis of the low-rank matrix approximation

The error analysis of the low-rank matrix approximation is performed in two steps: (1) define the bound of the approximation error and (2) choose $c$ and define the sampling probabilities.

The low-rank approximation incurs two error components: first, the error between $\hat{\mathbf{C}}$ and $\hat{\mathbf{C}}_{1P2}$ (assuming $\hat{\mathbf{C}}_{1P2}$ is the most optimal $P^2$-rank approximation to matrix $\hat{\mathbf{C}}$); second, $\mathbf{H}_{P2}\mathbf{H}_{P2}^T\hat{\mathbf{C}}$ incurs an error depending on $\|\hat{\mathbf{C}}\hat{\mathbf{C}}^T - \hat{\mathbf{C}}_1\hat{\mathbf{C}}_1^T\|_F$ or, equivalently, $\|\hat{\mathbf{C}}\|_F^2$. It has been proven⁷ that the error appears entirely within the following two bounds with respect to the Frobenius norm and spectral norm, respectively,

$$ \|\hat{\mathbf{C}} - \mathbf{H}_{P2}\mathbf{H}_{P2}^T\hat{\mathbf{C}}\|_F^2 \le \|\hat{\mathbf{C}} - \hat{\mathbf{C}}_{1P2}\|_F^2 + 2P\|\hat{\mathbf{C}}\hat{\mathbf{C}}^T - \hat{\mathbf{C}}_1\hat{\mathbf{C}}_1^T\|_F, \quad (8) $$

$$ \|\hat{\mathbf{C}} - \mathbf{H}_{P2}\mathbf{H}_{P2}^T\hat{\mathbf{C}}\|_2^2 \le \|\hat{\mathbf{C}} - \hat{\mathbf{C}}_{1P2}\|_2^2 + 2\|\hat{\mathbf{C}}\hat{\mathbf{C}}^T - \hat{\mathbf{C}}_1\hat{\mathbf{C}}_1^T\|_2, \quad (9) $$

where $\|\|_F$ indicates the Frobenius norm and $\|\|_2$ indicates the spectral norm. The optimal $P^2$-rank approximation matrix is $\hat{\mathbf{C}}_{1P2} = \mathbf{U}_{P2}\mathbf{U}_{P2}^T\hat{\mathbf{C}}_1$, where $\mathbf{U}_{P2}$ is composed of the top $P^2$ left singular vectors of $\hat{\mathbf{C}}_1$.

To control the approximation error, the method used to select $c$ and to define the sampling probabilities is important. According to the theory of Drineas et al.,⁷ we can determine the sampling probabilities as follows. We assume $H_p^2$ is constructed by sampling $c$ columns of $\hat{\mathbf{C}}$ with probabilities $q_i$ ($i = 1, \dots, (MF)^2$) to make $q_i \ge \beta|\hat{\mathbf{C}}^i|^2 / \|\hat{\mathbf{C}}\|_F^2$ for some positive $\beta \le 1$. Set $\eta = 1 + \sqrt{(8/\beta) \log(1/\delta)}$. $\epsilon$ represents an additional error coefficient and $\delta$ denotes the failure probability.

* If a Frobenius norm bound is desired: set $q_i = |\hat{\mathbf{C}}^{(i)}|^2 / \|\hat{\mathbf{C}}\|_F^2$, $c \ge 4P^2\eta^2/\epsilon^2$; then, with probability of at least $1 - \delta$,

$$ \|\hat{\mathbf{C}} - \mathbf{H}_{P2}\mathbf{H}_{P2}^T\hat{\mathbf{C}}\|_F^2 \le \|\hat{\mathbf{C}} - \hat{\mathbf{C}}_{1P2}\|_F^2 + \epsilon\|\hat{\mathbf{C}}\|_F^2. \quad (10) $$

* If a spectral norm bound is desired: set $q_i = |\hat{\mathbf{C}}^{(i)}|^2 / \|\hat{\mathbf{C}}\|_F^2$, $c \ge 4\eta^2/\epsilon^2$; then, with probability of at least $1 - \delta$,

$$ \|\hat{\mathbf{C}} - \mathbf{H}_{P2}\mathbf{H}_{P2}^T\hat{\mathbf{C}}\|_2^2 \le \|\hat{\mathbf{C}} - \hat{\mathbf{C}}_{1P2}\|_2^2 + \epsilon\|\hat{\mathbf{C}}\|_F^2. \quad (11) $$

It is necessary to note that the approximation error can be made arbitrarily small by choosing enough columns.

### 4. Analysis of the implementation and running time of the low-rank approximation

In the proposed algorithm, we can sample one column in one pass over the matrix without previously knowing its size. With the algorithm running in parallel, we can obtain more samples in a single pass in $O(c)$ additional space and time. Given the index to be sampled, we need additional $O((MF)^2c)$ space and time to construct matrix $\hat{\mathbf{C}}_1$. Computing $\hat{\mathbf{C}}_1^T\hat{\mathbf{C}}_1$ requires $O((MF)^2c)$ additional space and $O(((MF)^2c)^2)$ additional time, and we consider the SVD of $\hat{\mathbf{C}}_1^T\hat{\mathbf{C}}_1$ to require $O(c^2)$ additional space and $O(c^3)$ additional time. Then, computing $\mathbf{H}_{P2}$ requires $P^2$ matrix-vector

<table><caption>Table 1. The configuration parameters of simulation, small-scale experiment and the at-sea experiment.</caption><thead><tr><th></th><th><i>M</i></th><th><i>z</i><sub>s</sub> (m)</th><th><i>z</i><sub>r<sub>1</sub></sub> ~ <i>z</i><sub>r<sub>M</sub></sub> (m)</th><th><i>d</i> (m)</th><th><i>D</i> (m)</th></tr></thead><tbody><tr><td>Simulation</td><td>6</td><td>25</td><td>22.75 - 26.5</td><td>0.75</td><td>750</td></tr><tr><td>Small-scale experiment</td><td>11</td><td>27.125 × 10<sup>-3</sup></td><td>12.125 × 10<sup>-3</sup> - 42.125 × 10<sup>-3</sup></td><td>0.3 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>Ocean data</td><td>32</td><td>96</td><td>80-112</td><td>0.75</td><td>4.701 × 10<sup>3</sup></td></tr></tbody></table>
<table><caption>Table 2. The test parameters of simulation, small-scale experiment and the at-sea experiment.</caption><thead><tr><th></th><th>ν<sub>c</sub> (Hz)</th><th>ν<sub>w</sub> (Hz)</th><th>ν<sub>n</sub>(ν<sub>b</sub> - ν<sub>e</sub>) (Hz)</th><th>N</th><th>δ</th><th>ε</th></tr></thead><tbody><tr><td>Simulation</td><td>3 × 10<sup>3</sup></td><td>3 × 10<sup>3</sup></td><td>39(2 × 10<sup>3</sup> - 5 × 10<sup>3</sup>)</td><td>128</td><td>0.01</td><td>2</td></tr><tr><td>Small-scale experiment</td><td>1.2 × 10<sup>6</sup></td><td>1 × 10<sup>6</sup></td><td>17(0.48276 × 10<sup>6</sup> - 1.5862 × 10<sup>6</sup>)</td><td>145</td><td>0.01</td><td>4</td></tr><tr><td>Ocean data</td><td>3.2 × 10<sup>3</sup></td><td>1 × 10<sup>3</sup></td><td>7(2.8 × 10<sup>3</sup> - 4 × 10<sup>3</sup>)</td><td>60</td><td>0.01</td><td>2</td></tr></tbody></table>

multiplications for $O((MFP)^2c)$ additional space and time. Overall, the proposed fast algorithm requires $O((MF)^2c + c^2)$ additional space and $O((MF)^2c^2 + c^3)$ additional time, which have linear relationships with the number of rows in the high-order cumulants square matrix.

## 5. Performance test

In this section, we test the performance of the proposed algorithm with many sets of simulation data, a set of real data obtained in a tank⁹ (the experiment is based on the principle that if we amplify the frequency of the signals by a factor while reducing the spatial distance by the same factor, the actual physical phenomena occurring in nature can be reproduced on a smaller scale) and finally a set of ocean data.¹⁰ (The experiment was performed in July 2005 north of Elba Island, Italy.) The conventional higher-order algorithm in Ref. 5 is provided for comparison. These tests are performed on a server with 6 Intel Xeon CPUs(X7542) with 48 cores @2.67 GHz and 256.0 GB RAM on the Redhat Linux 5.6 operating system. All the tests are performed between a point source and a vertical array, which are composed of *M* equally spaced sensors. We show the parameters including the number of sensors (*M*), the source location (*z*<sub>s</sub>), the range of the array covering the ocean from the position at *z*<sub>r1</sub> m to the position at *z*<sub>rM</sub>, the interval between two adjacent sensors *d* and the distance between the source and the reference sensor *D* in each configuration in Table 1. We give the central frequency (*ν*<sub>c</sub>), the bandwidth of the emitted signals (*ν*<sub>w</sub>), the number of the frequencies (*ν*<sub>n</sub>: from *ν*<sub>b</sub> to *ν*<sub>e</sub>) and the taps of the broadband time-domain signal. (*N*) used in the tests in Table 2.

<figure><img src="image_7.png" /><figcaption>Fig. 1. (Color online) Separation results with the proposed fast algorithm. Black crosses indicate the expected positions of the raypaths. (A) Separation results with the simulation data. (B) The time reduction ratio with the simulation data (the proposed algorithm vs 4-smoothing-MUSICAL). (C) Separation results with the small-scale data.</figcaption></figure>
The proposed fast algorithm is applied to the simulation data, a set of real data obtained in a tank and finally a set of ocean data and the separation results are shown in Fig. 1 and Fig. 2, respectively. The expected positions are marked with black crosses. Figure 1(A) illustrates the separation result with a set of simulation data. We can see that the five raypaths are all successfully separated with few errors compared with their expected values [(0.5, 0), (0.9, -0.09), (1.2, 0.15), (1.9, -0.254), and (2.5, 0.235)]. Moreover, the efficiency of the proposed algorithm is tested with many sets of simulation data for various numbers of snapshots. Ten realizations are used for each number of snapshots. Figure 1(B) shows that the time reduction ratio increases sharply with increasing number of snapshots. Figure 1(C) shows the separation result with a set of real data obtained in a tank. The proposed algorithm successfully separates the seven raypaths with few errors. Their expected values are (0.6789, 0), (0.6799, 3.1), (0.68, -3.2), (0.683, 6.3), (0.683, -6.3), (0.688, 9.3), and (0.6882, -9.4). Figures 2(A) and 2(B) show the separation results for the at-sea experimental data with the beamforming algorithm and our proposed algorithm, respectively. The proposed fast algorithm are more efficient for noise suppression and obtain the smaller spots representing the positions of the raypaths compared to the beamforming algorithm. Figure 2(C) displays the three raypaths in Figs. 2(A) and 2(B), which propagate between the centers of source and receive arrays. The centers of source and receive arrays both located at 96 m in the water. Due to the sound speed variation described in the Ref. 10, the three raypaths arrived at close time. There are a surface-reflected ray and two refracted rays. Each refracted ray has a turning point near the depth of maximum sound-speed variability. Figures 2(a) and 2(b) illustrate that the proposed fast algorithm is robust to the fluctuation of environmental factors, for instance the water column sound speed variation we have described. The performance of the 4-smoothing-MUSICAL algorithm and the proposed fast algorithm are almost the same, which are not included in the

<figure><img src="image_8.png" /><figcaption>Fig. 2. (Color online) Results comparison for the data obtained from at-sea experiment. Black crosses indicate the expected positions of the raypaths. (A) Separation results with the beamforming algorithm. (B) Separation results with the proposed fast algorithm. (C) Raypaths between the source and the reference sensor for the three rays.</figcaption></figure>
<table><caption>Table 3. The comparison of time and space consumption.</caption><thead><tr><th></th><th></th><th>Conventional EVD</th><th>Linear time EVD</th><th>Reduction ratio</th></tr></thead><tbody><tr><td rowspan="2">Simulation</td><td>CPU time (s)</td><td>5195.208</td><td>19.642</td><td>99.62%</td></tr><tr><td>Peak memory (Kb)</td><td>15 167 030</td><td>682 820</td><td>95.5%</td></tr><tr><td rowspan="2">Small-scale experiment</td><td>CPU time (s)</td><td>19 042</td><td>58</td><td>99.7%</td></tr><tr><td>Peak memory (Kb)</td><td>38 797 300</td><td>291 500</td><td>99.2%</td></tr><tr><td rowspan="2">Ocean data</td><td>CPU time (s)</td><td>55 080</td><td>55</td><td>99.9%</td></tr><tr><td>Peak memory (Gb)</td><td>95.3</td><td>0.235</td><td>99.75%</td></tr></tbody></table>

paper due to space limitations. Table 3 compares the time and space consumed by the 4-smoothing-MUSICAL algorithm and the proposed fast algorithm, illustrating that with the proposed algorithm, the cost of the time in decomposition step has been reduced by approximately 99.62% (simulation), 99.7% (small-scale experiment), and 99.9% (ocean data), respectively. The peak memory has been reduced by approximately 95.50% (simulation), 99.2% (small-scale experiment), and 99.75% (ocean data), respectively.

## 6. Conclusion

In this paper we proposed a fast raypath separation method that achieves large reductions in time and space consumption compared with the conventional high-order cumulant-based algorithm. In the future, the method will be extended to improve the efficiency of the multi-dimensional subspace-based algorithms.

## Acknowledgments

This paper was supported by the National Natural Science Foundation of China (Grant Nos. 61401085 and 31400842), the State Key Laboratory of Acoustics, Chinese Academy of Sciences (Grant No. SKLA201604), and the Scientific Research Foundation for the Returned Overseas Chinese Scholars. The Focused Acoustic Forecasting experiment (FAF05) experiment was performed in collaborative experiments with the NATO Underwater Research Centre (NURC), La Spezia, Italy, with Mark Stevenson as Chief Scientist. Scientists who contributed to these experiments include Tuncay Akal, W. A. Kuperman, W. H. Hodgkiss, H. C. Song, B. D. Cornuelle, Piero Boni, Piero Guerrini, other NURC staff, and the officers and crew of the RV Alliance.

## References and links

¹W. Munk, P. Worcester, and C. Wunsch, “Ocean acoustic tomography,” in *Ocean Acoustic Tomography* (Cambridge University Press, San Diego, 1995), Chap. 2, pp. 30–114.

²R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” *IEEE Trans. Antennas Propag.* 34(3), 276–280 (1986).

³L. Jiang, F. Aulanier, G. Le Touzé, B. Nicolas, and J. Mars, “Raypath separation with high resolution processing,” in *OCEANS’11 IEEE Santander-Oceans of Energy for a Sustainable Future* (2011), Vol. 1, pp. 1–4.

⁴L. Jiang, P. Roux, and J. I. Mars, “Raypath separation with a high-resolution algorithm in a shallow-water waveguide,” *IEEE J. Ocean. Eng.* 43(1), 119–130 (2018).

⁵L. Jiang, Y. Hong, P. Roux, J. Wu, and H. Shu, “Active wideband higher-order raypath separation in multipath environment,” *J. Acoust. Soc. Am.* 141(1), EL38–EL44 (2017).

⁶A. Frieze, R. Kannan, and S. Vempala, “Fast Monte-Carlo algorithms for finding low-rank approximations,” *J. ACM* 51(6), 1025–1041 (2004).

⁷P. Drineas, R. Kannan, and M. W. Mahoney, “Fast Monte Carlo algorithms for matrices. II: Computing a low-rank approximation to a matrix,” *SIAM J. Comput.* 36(1), 158–183 (2006).

⁸P. Drineas, R. Kannan, and M. W. Mahoney, “Fast Monte Carlo algorithms for matrices. I: Approximating matrix multiplication,” *SIAM J. Comput.* 36(1), 132–157 (2006).

⁹P. Roux, W. A. Kuperman, W. S. Hodgkiss, H. C. Song, T. Akal, and M. Stevenson, “A nonreciprocal implementation of time reversal in the ocean,” *J. Acoust. Soc. Am.* 116, 1009–1015 (2004).

¹⁰P. Roux, B. D. Cornuelle, W. A. Kuperman, and W. S. Hodgkiss, “The structure of raylike arrivals in a shallow-water waveguide,” *J. Acoust. Soc. Am.* 124(6), 3430–3439 (2008).