# Multivariate Decomposition of Acoustic Signals in Dispersive Channels

Miloš Brajović, Isidora Stanković, Jonatan Lerga, Cornel Ioana, Eftim Zdravevski, Miloš Daković

► To cite this version:

Miloš Brajović, Isidora Stanković, Jonatan Lerga, Cornel Ioana, Eftim Zdravevski, et al. Multi-variate Decomposition of Acoustic Signals in Dispersive Channels. Mathematics, 2021, Special Issue New Trends in Graph and Complexity Based Data Analysis and Processing, 9 (21), pp.2796. 10.3390/math9212796. [https://arxiv.org/abs/2003.09174](https://arxiv.org/abs/2003.09174)

HAL Id: hal-03589944

https://hal.science/hal-03589944v1

Submitted on 10 Nov 2024

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>
Article

# Multivariate Decomposition of Acoustic Signals in Dispersive Channels

Miloš Brajović<fnref n="1" />,*, Isidora Stanković <fnref n="1" />,*, Jonatan Lerga <fnref n="2" />,*, Cornel Ioana <fnref n="3" />,*, Eftim Zdravevski <fnref n="4" />
and Miloš Daković <fnref n="1" />

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

check for
updates

**Citation:** Brajović, M.; Stanković, I.; Lerga, J.; Ioana, C.; Zdravevski, E.; Daković, M. Multivariate Decomposition of Acoustic Signals in Dispersive Channels. Mathematics **2021**, *9*, 2796. https://doi.org/10.3390/math9212796

Academic Editor: João Nuno Prata

Received: 23 September 2021

Accepted: 28 October 2021

Published: 4 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** We present a signal decomposition procedure, which separates modes into individual components while preserving their integrity, in effort to tackle the challenges related to the characterization of modes in an acoustic dispersive environment. With this approach, each mode can be analyzed and processed individually, which carries opportunities for new insights into their characterization possibilities. The proposed methodology is based on the eigenanalysis of the autocorrelation matrix of the analyzed signal. When eigenvectors of this matrix are properly linearly combined, each signal component can be separately reconstructed. A proper linear combination is determined based on the minimization of concentration measures calculated exploiting time-frequency representations. In this paper, we engage a steepest-descent-like algorithm for the minimization process. Numerical results support the theory and indicate the applicability of the proposed methodology in the decomposition of acoustic signals in dispersive channels.

**Keywords:** concentration measures; dispersive channels; multivariate signals; non-stationary signals; multicomponent signal decomposition

## 1. Introduction

Signals with time-varying spectral content, known as non-stationary signals, are analyzed using time-frequency signal (TF) signal analysis [1–17]. Some commonly used TF representations include short-time Fourier transform (STFT) [1,3], pseudo-Wigner distribution (WD) [1,9,12], and S-method (SM) [3]. Time-scale, multi-resolution analysis using the wavelet transform is an additional approach to characterize non-stationary signal behavior [4]. Various representations are primarily applied in the instantaneous frequency (IF) estimation and related applications [8–15], since they concentrate the energy of a signal component at and around the respective instantaneous frequency. Concentration measures provide a quantitative description of the signal concentration in the given representation domain [18], and can be used to assess the area of the time-frequency plane covered by a signal component.

In order to characterize multicomponent signals, it is quite common to perform signal decomposition, which assumes that each individual component is extracted for separate analysis, such as for the IF estimation. Decomposition techniques for multicomponent signals are quite efficient if components do not overlap in the time-frequency plane [19–26]. The method originally presented in [26] can be used to completely extract each component by using an intrinsic relation between the PWD and the SM. In the analysis of multicomponent signals, it is, however, common that various components partially overlap in the time-frequency plane, making the decomposition process particularly challenging [19–26].

<footnote n="1">¹ Faculty of Electrical Engineering, University of Montenegro, 81000 Podgorica, Montenegro; isidoras@ucg.ac.me (I.S.); milos@ucg.ac.me (M.D.)</footnote>

<footnote n="2">² Faculty of Engineering, University of Rijeka, 51000 Rijeka, Croatia</footnote>

<footnote n="3">³ Gipsa-Lab, Université Grenoble Alpes, 38400 Grenoble, France; cornel.ioana@gipsa-lab.grenoble-inp.fr</footnote>

<footnote n="4">⁴ Faculty of Computer Science and Engineering, University Ss. Cyril and Methodius, 1000 Skopje, North Macedonia; eftim.zdravevski@finki.ukim.mk</footnote>

<footnote>* Correspondence: milosb@ucg.ac.me (M.B.); jlerga@riteh.hr (J.L.)</footnote>
In this rather unfavorable scenario, overlapped components partially share the same domains of supports, and existing decomposition techniques provide only partial results in the univariate case, limited to very narrow signal classes. For example, linear frequency modulated signals are decomposed using the chirplet transform, Radon transform, or similar techniques [20,25], whereas sinusoidally modulated signals are separated using the inverse Radon transform [27]. However, these techniques cannot perform the decomposition when components have a general, non-stationary form.

In the multivariate (multichannel) framework, it is assumed that the signals are acquired using multiple sensors, [28–44]. The sensors modify component amplitudes and phases. However, the interdependence of values from various channels can be utilized in the signal decomposition. This concept has also been exploited in the empirical mode decomposition (EMD) [39–43]. It was previously shown that WD-based decomposition is possible if signals are available in the multivariate form [28–30]. Moreover, the decomposition can be performed by directly engaging the eigenanalysis of the auto-correlation matrix, calculated for signals in the multivariate form [31–34]. It should also be noted that the problem of multicomponent signal decomposition has some similarities with the blind source separation [45–48]. However, the basic difference is in the aim to extract each signal component in the decomposition framework, whereas in the blind source separation, the aim is to separate signal sources (although one source may generate several components). The mixing scheme from the blind source separation framework is used in a recently proposed mode decomposition approach [49]. Another line of the decomposition-related research includes mode decomposition techniques, which could be used for separation of modal responses and identification of progressive changes in modal parameters [50].

Overlapped components pose a challenge in various applications, such as in biomedical signal processing [44,51,52], radar signal processing [53], and processing of lamb waves [54]. Popular approaches, such as the EMD and multivariate EMD (MEMD), [39–43] cannot respond to the challenges posed by components overlapped in the time-frequency plane and do not provide acceptable decomposition results in this particular case [28]. Additionally, the applicability of these methods is highly influenced by amplitude variations of the signal components. In this paper, we present a framework for the decomposition of acoustic dispersive environment signals into individual modes based on the multivariate decomposition of multicomponent non-stationary signals. Even when simple signal forms are transmitted, acoustic signals in dispersive channels appear in the multicomponent form, with either very close or partially overlapped components. Being reflected from the underwater surfaces and objects, each individual component carries information about the underwater environment. That information is inaccessible while the signal is in its multicomponent form. This makes analyzing acoustic signals (mainly their localization and characterization) a challenging problem for research [55–60]. The presented decomposition approach enables complete separation of components and their individual characterization (e.g., IF estimation, based on which knowledge regarding the underwater environment can be acquired).

We aim at solving this notoriously difficult practical problem by exploiting the interdependencies of multiply acquired signals: such signals can be considered as multivariate and are subject to slight phase changes across various channels, occurring due to different sensing positions and due to various physical phenomena, such as water ripples, uneven seabed, and changes in the seabed substrate. As each eigenvector of the autocorrelation matrix of the input signal represents a linear combination of the signal components [31,33], slight phase changes across the various channels are actually favorable for forming an undetermined set of linearly independent equations relating the eigenvectors and the components. Moreover, we have previously shown that each component is a linear combination of several eigenvectors corresponding to the largest eigenvalues, with unknown weights [31] (the number of these eigenvalues is equal to the number of signal components). Among infinitely many possible combinations of eigenvectors, the aim is to find the weights producing the most concentrated combination, as each individual signal compo
nent (mode) is more concentrated than any linear combination of components, as discussed in detail in [31]. Therefore, we engage concentration measures [18] to set the optimization criterion and perform the minimization in the space of the weights of linear combinations of eigenvectors.

We revisit our previous research from [28,31,33], and the main contributions are two-fold. The decomposition principles of the auto-correlation matrix [31,33] are reconsidered. Instead of exploiting direct search [31] or a genetic algorithm [33], we show that the minimization of concentration measure in the space of complex-valued coefficients acting as weights of eigenvectors, which are linearly combined to form the components, can be performed using a steepest-descent-based methodology, originally used in the decomposition from [28]. The second contribution is the consideration of a practical application of the decomposition methodology.

The paper is organized as follows. After the Introduction, we present the basic theory behind the considered acoustic dispersive environment in Section 2. Section 3 presents the principles of multivariate signal decomposition of dispersive acoustic signals. The decomposition algorithm is summarized in Section 4. The theory is verified on numerical examples and additionally discussed in Section 5. Whereas the paper ends with concluding remarks.

## 2. Dispersive Channels and Shallow Water Theory

Our primary goal is the decomposition of signals transmitted through dispersive channels. Decomposition assumes the separation of signal components while preserving the integrity of each component. Signals transmitted through dispersive channels are multicomponent and non-stationary, even in cases when emitted signals have a simple form. This makes the challenging problem of decomposition, localization, and characterization of such signals a fairly studied topic [55-67]. The decomposition can be performed using the time-frequency phase-continuity of the signals [55], or using the mode characteristics of the signal [56]. After being transmitted through a dispersive environment, measured signals consist of several components called modes. The non-stationarity of these modes is a consequence of frequency dependent properties of the signal propagation media.

The dispersive acoustic environment is commonly studied within the context of shallow waters, defined by the depth of the sea/ocean, which are less than $D = 200 \text{ m}$ [55,57-67]. The speed of signals traveling through water is affected by many factors, such as the salinity, the temperature, or the pressure of the water, but it is usually approximated around 1480-1500 m/s. Note that this speed is larger than the speed of signals traveling through the air, which is estimated at approximately 340-360 m/s. Such setups typically have extremely complex analyses. Moreover, bottom properties and water volume add up to this complexity, as well as noise caused by activities on the water surface and on the coastlines (commonly related to cavitation). Dispersivity of shallow waters occurs due to many reasons, among which are the roughness of the bottom, strength of the waves, and cavity level. Dispersive channels have varying frequency characteristics (phase and spectral content) during the transmission of the signal.

### 2.1. Normal Mode Solution

The propagation of sound in a shallow water environment is mathematically represented by the wave equations. Among several methods of deriving the solution of the wave equation, the most commonly used is the normal mode solution, based on solving depth-dependent equations using the method of variable separation. Further analysis will be developed based on the isovelocity waveguide model presented in Figure 1, which characterizes a rigid boundary of the seabed. This further yields to an ideally spread velocity $c$. Furthermore, channel models assume that the structure of a channel is a waveguide, where multiple normal-modes are received as delayed and scaled versions of the transmitted signal [56,58,59,65]. Our aim is to decompose the received signal by extracting each mode
separately. Such extracted modes can be used in further processing, such as IF estimation, characterization, and classification.

More general models assume a more complicated environment, where the boundary of the bottom depends on the nature of the ocean, such as the roughness, depending on the weather conditions and different environments in the ocean itself. These models take into account the scattering of the transmitted signal as well. Our future work will be oriented towards these models as well.

<figure><img src="image_4.png" /><figcaption>Figure 1. The considered underwater isovelocity setup. Water depth is D, transmitter depth is zt, the receiver is positioned at depth zr, and the transmitter-receiver range is r.</figcaption></figure>

## 2.2. Problem Formulation—Signal Processing Approach

The practical setup, shown in Figure 1 is further considered. In this setup, it is assumed that the transmitter is located in the water at the depth of zt, whereas the receiver is located at the depth of zr meters. It is assumed that the wave is transmitted through an isovelocity channel as in [55–57,61–63,67]. The distance between the transmitter and the receiver is r.

Taking into account the spectrum of the received signal, in the normal-mode case, the transfer function reads

$$H(\omega) = \sum_{p=1}^{+\infty} G_p(z_t) G_p(z_r) \frac{\exp(jk_r(p, \omega)r)}{\sqrt{k_r(p, \omega)r}} = \sum_{p=1}^{+\infty} A_t(p, \omega) \exp(jk_r(p, \omega)r), \quad (1)$$

with $G_p(z_t)$ and $G_p(z_r)$ being the modal functions of the p-th mode corresponding to the transmitter and the receiver [55,56,65], with the attenuation rate is $A_t(p, \omega) = A(p, \omega) / \sqrt{r}$. Angular frequency is denoted by $\omega$. The modes are dependent on wavenumbers $k_r(p, \omega)$ [55]

$$k_r^2(p, \omega) = \left(\frac{\omega}{c}\right)^2 - \left((p - 0.5)\frac{\pi}{D}\right)^2. \quad (2)$$

The multicomponent structure of the transfer function is dependent on the number of modes. The speed of sound propagation underwater is c = 1500 m/s.

The response to a monochromatic signal

$$s(n) = \exp(j\omega_0 n) \quad (3)$$

at the p-th mode can be written as

$$s_p(n) \approx A_t(p, \omega_0) \exp(j\omega_0 n - jk_r(p, \omega_0)r). \quad (4)$$

The phase velocity of this signal is

$$v_p(\omega) = \frac{\omega}{k_r(p, \omega)} = \frac{\omega}{\sqrt{\left(\frac{\omega}{c}\right)^2 - \left((p - 0.5)\frac{\pi}{D}\right)^2}}. \quad (5)$$
This is the horizontal velocity of the corresponding phase for the p-th mode. The energy propagation of the signal component is represented by the group velocity

$$u_p(\omega) = \frac{dr(t)}{dt} = \frac{d\omega}{dk_r(p, \omega)} = \frac{1}{\frac{dk_r(p, \omega)}{d\omega}} = \frac{1}{\frac{d}{d\omega}\sqrt{\left(\frac{\omega}{c}\right)^2 - \left((p - 0.5)\frac{\pi}{D}\right)^2}} \quad (6)$$

The received signal can be represented in the Fourier transform domain as a product of the Fourier transform of the transmitted signal, $S(\omega)$ and the transfer function $H(\omega)$ of the channel in the normal-mode form; that is

$$X(\omega) = S(\omega)H(\omega). \quad (7)$$

In time domain, the received signal, $x(n)$, is the convolution of the transmitted signal, $s(n)$ and the impulse response, $h(n)$, from (1), i.e.,

$$x(n) = s(n) * h(n). \quad (8)$$

In the following sections, we present an efficient methodology for the decomposition of mode functions, which will make the problem of detecting and estimating the received signal parameters straightforward.

### 3. Multivariate Decomposition

#### 3.1. Multivariate (Multichannel) Signals

Multivariate or multichannel signals are acquired using multiple sensors. It is further assumed that C sensors at the receiving position are used for the acquisition of signal $x_R(n)$. Here, subscript R is used to denote the fact that the acquired signal is real-valued. All C sensors placed at the depth $z_r$ are part of the receiver. In the range direction, sensor distances from the transmitter are $r + \delta_c$, $c = 1, 2, \dots, C$. Deviations $\delta_c$, $c = 1, 2, \dots, C$, are small as compared to the distance, $r$, between the transmitter and receiver locations in Figure 1 in range direction.

Since the measured signal, $x_R(n)$, is real-valued, its analytic extension

$$x(n) = x_R(n) + j\mathcal{H}\{x_R(n)\} \quad (9)$$

is assumed in the further multivariate decomposition setup, where $\mathcal{H}\{x_R(n)\}$ is the Hilbert transform of this signal. This analytic form assumes only non-negative frequencies. Each sensor modifies the amplitude and the phase of the acquired signal. Therefore, the channel signals take the form $a_c(n) \exp(j\phi_c(n)) = \alpha_c \exp(j\varphi_c) x(n)$, for each sensor $c = 1, 2, \dots, C$. When a monocomponent signal $x(n) = A(n) \exp(j\psi(n))$, is measured at sensor c, this yields

$$a_c(n) \exp(j\phi_c(n)) = \alpha_c \exp(j\varphi_c) x(n),$$

or $a_c(n) \cos(\phi_c(n))$ in the case of a real-valued signal. The corresponding analytic signal, $a_i(n) \exp(j\phi_i(n)) = a_i(n) \cos(\phi_i(n)) + j\mathcal{H}\{a_i(n) \cos(\phi_i(n))\}$ is a valid representation of the real amplitude-phase signal $a_c(n) \cos(\phi_c(n))$ if the spectrum of $a_i(n)$ is nonzero only within the frequency range $|\omega| < B$ and the spectrum of $\cos(\phi_i(n))$ occupies a non-overlapping (much) higher frequency range [5]. If variations of the amplitude, $a_c(n)$, are much slower than the phase $\phi_c(n)$ variations, then this signal is monocomponent [31]. A unified representation of a multichannel (multivariate) signal, acquired using C sensors, assumes the following vector form

$$\mathbf{x}(n) = \begin{bmatrix} x^{(1)}(n) \\ x^{(2)}(n) \\ \vdots \\ x^{(C)}(n) \end{bmatrix} = \begin{bmatrix} a_1(n)e^{j\phi_1(n)} \\ a_2(n)e^{j\phi_2(n)} \\ \vdots \\ a_C(n)e^{j\phi_C(n)} \end{bmatrix}, \quad n = 1, 2, \dots, N. \quad (10)$$
## 3.2. Multivariate Multicomponent Signals

When the measured signal consists of a linear combination of P linearly independent components $s_p(n) = A_p(n)e^{j\psi_p(n)}$, $p = 1, 2, \dots, P$, then it is commonly referred to as a multicomponent signal

$$x(n) = \sum_{p=1}^{P} s_p(n) = \sum_{p=1}^{P} A_p(n) e^{j\psi_p(n)}. \quad (11)$$

Component amplitudes, $A_p(n)$, are characterized by slow-varying dynamics as compared to the variations of the component phases, $\psi_p(n)$. Linear independence of the components assumes that neither component can be represented as a linear combination of other components for any considered time instant $n$.

Incorporation of multicomponent signal definition (11) into the multichannel form (10), yields

$$\mathbf{x}(n) = \begin{bmatrix} x^{(1)}(n) \\ x^{(2)}(n) \\ \vdots \\ x^{(C)}(n) \end{bmatrix} = \begin{bmatrix} a_1(n) e^{j\phi_1(n)} \\ a_2(n) e^{j\phi_2(n)} \\ \vdots \\ a_C(n) e^{j\phi_C(n)} \end{bmatrix} = \begin{bmatrix} \sum_{p=1}^{P} \alpha_{1p} s_p(n) e^{j\phi_{1p}} \\ \sum_{p=1}^{P} \alpha_{2p} s_p(n) e^{j\phi_{2p}} \\ \vdots \\ \sum_{p=1}^{P} \alpha_{Cp} s_p(n) e^{j\phi_{Cp}} \end{bmatrix}, \quad n = 1, 2, \dots, N, \quad (12)$$

or, more briefly,

$$\begin{bmatrix} x^{(1)}(n) \\ x^{(2)}(n) \\ \vdots \\ x^{(C)}(n) \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1P} \\ a_{21} & a_{22} & \cdots & a_{2P} \\ \vdots & \vdots & \ddots & \vdots \\ a_{C1} & a_{C2} & \cdots & a_{CP} \end{bmatrix} \begin{bmatrix} s_1(n) \\ s_2(n) \\ \vdots \\ s_P(n) \end{bmatrix}, \quad (13)$$

that is

$$\mathbf{x}(n) = \mathbf{A}\mathbf{s}(n), \quad (14)$$

where the vector of signal components, $\mathbf{s}(n)$ is, for instant $n$, given by

$$\mathbf{s}(n) = \begin{bmatrix} s_1(n) \\ s_2(n) \\ \vdots \\ s_P(n) \end{bmatrix}, \quad n = 1, 2, \dots, N. \quad (15)$$

Matrix $\mathbf{A}$ of size $C \times P$, which relates the signal in the c-th channel, $x^{(c)}(n)$ with signal components, $s_p(n)$, in form of a linear combination

$$x^{(c)}(n) = \sum_{p=1}^{P} a_{cp} s_p(n) = \sum_{p=1}^{P} \alpha_{cp} s_p(n) e^{j\phi_{cp}} \quad (16)$$

has the following form

$$\mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1P} \\ a_{21} & a_{22} & \cdots & a_{2P} \\ \vdots & \vdots & \ddots & \vdots \\ a_{C1} & a_{C2} & \cdots & a_{CP} \end{bmatrix},$$
with elements being complex constants $a_{cp} = \alpha_{cpe} e^{j q_{cp}}$, $c = 1, 2, \dots, C$, $p = 1, 2, \dots, P$. These constants linearly relate the channel signals with signal components. Clearly, the maximum number of independent channels $x^{(1)}(n), x^{(2)}(n), \dots, x^{(C)}(n)$ in $\mathbf{x}(n)$ is

$$M = \min\{C, P\}, \qquad (17)$$

since $\min\{\mathbf{A}\} \le \min\{C, P\}$.

The relation between the C measured channel signals, $x^{(c)}(n)$, and P components, $x_p(n)$, can be, taking into consideration all time instants, formed by introducing $C \times N$ matrix $\mathbf{X}_{sen}$ with elements being the sensed signal values, and $\mathbf{X}_{com}$ comprising the samples of signal components $s_p(n)$. In that case, the relation is

$$\begin{bmatrix} x^{(1)}(1) & \dots & x^{(1)}(N) \\ x^{(2)}(1) & \dots & x^{(2)}(N) \\ \vdots & \ddots & \vdots \\ x^{(C)}(1) & \dots & x^{(C)}(N) \end{bmatrix} = \mathbf{A} \begin{bmatrix} s_1(1) & \dots & s_1(N) \\ s_2(1) & \dots & s_2(N) \\ \vdots & \ddots & \vdots \\ s_p(1) & \dots & s_p(N) \end{bmatrix}. \qquad (18)$$

or

$$\mathbf{X}_{sen} = \mathbf{A}\mathbf{X}_{com}. \qquad (19)$$

Now we can introduce the autocorrelation matrix $\mathbf{R}$ of the sensed signal, whose eigenvectors will be used in the multivariate decomposition framework:

$$\mathbf{R} = \mathbf{X}_{sen}^H \mathbf{X}_{sen}, \qquad (20)$$

where $(\cdot)^H$ denotes the Hermitian transpose. Individually, elements of this matrix are products of $\mathbf{x}(n_1)$ and $\mathbf{x}^H(n_1)$ at given instants $n_1$ and $n_2$:

$$R(n_1, n_2) = \mathbf{x}^H(n_2)\mathbf{x}(n_1) = \sum_{i=1}^{C} x^{(i)*}(n_2)x^{(i)}(n_1), \qquad (21)$$

where $\mathbf{x}(n_1) = [x^{(1)}(n_1) \ x^{(2)}(n_1) \ \dots \ x^{(C)}(n_1)]^T$ is the column vector of sensed values at a given instant $n_1$. As it will be shown next, the eigenvectors of the autocorrelation matrix, $\mathbf{R}$, corresponding to the largest eigenvalues, consist of linear combinations of signal components. This fact will be used to develop the algorithm for the extraction of those components.

### 3.3. Eigendecomposition of the Autocorrelation Matrix

It is well-known that any square matrix $\mathbf{R}$, of dimensions $K \times K$, can be subject of eigenvalue decomposition

$$\mathbf{R} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^H = \sum_{p=1}^{K} \lambda_p \mathbf{q}_p \mathbf{q}_p^H, \qquad (22)$$

with $\lambda_p$ being the eigenvalues and $\mathbf{q}_p$ being the corresponding eigenvectors of matrix $\mathbf{R}$. Matrix $\mathbf{\Lambda}$ contains eigenvalues $\lambda_p, p = 1, 2, \dots, K$ on the main diagonal and zeros at other positions. Matrix $\mathbf{Q} = [\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_K]$ contains the eigenvectors $\mathbf{q}_p$ as its columns. Since $\mathbf{R}$ is symmetric, eigenvectors are orthogonal.

From definition (20) and based on relation $\mathbf{X}_{sen} = \mathbf{A}\mathbf{X}_{com}$, autocorrelation matrix $\mathbf{R}$ can be rewritten as

$$\mathbf{R} = \mathbf{X}_{sen}^H \mathbf{X}_{sen} = \mathbf{X}_{com}^H \mathbf{A}^H \mathbf{A} \mathbf{X}_{com} = \sum_{i=1}^{p} \sum_{j=1}^{p} \tilde{a}_{ij} \mathbf{s}_i \mathbf{s}_j^H, \qquad (23)$$
where $\bar{a}_{ij}$ is used to denote elements of matrix $\mathbf{A}^H\mathbf{A}$ and $\mathbf{s}_i = [s_i(1), s_i(2), \dots, s_i(N)]^H$. Elements of matrix $\mathbf{R}$ are

$$R(n_1, n_2) = \sum_{i=1}^{P} \sum_{j=1}^{P} \bar{a}_{ij} s_i(n_1) s_j^*(n_2) = [s_1^*(n_2), s_2^*(n_2), \dots, s_P^*(n_2)] \mathbf{A}^H \mathbf{A} \begin{bmatrix} s_1(n_1) \\ s_2(n_1) \\ \vdots \\ s_P(n_1) \end{bmatrix}. \quad (24)$$

Based on the decomposition of matrix $\mathbf{R}$ on its eigenvalues and eigenvectors, we further have

$$\mathbf{R} = \sum_{p=1}^{M} \lambda_p \mathbf{q}_p \mathbf{q}_p^H = \sum_{i=1}^{P} \sum_{j=1}^{P} \bar{a}_{ij} \mathbf{s}_i \mathbf{s}_j^H, \quad (25)$$

with $M = \min\{C, P\}$. It will be further assumed that the number of sensors, $C$ is such that $C \ge P$. In that case, there are $M = P$ eigenvectors in (25). Two general cases can be further discussed:

*   **Non-overlapped components.** Note that the case when no components $\mathbf{s}_i$ and $\mathbf{s}_j$ overlap in the time-frequency plane implies that these components are orthogonal. In that case, the right side of (25) becomes:

$$\mathbf{R} = \sum_{i=1}^{P} \mathbf{s}_i \mathbf{s}_i^H \sum_{j=1}^{P} \bar{a}_{ij} = \kappa_p \sum_{p=1}^{P} \mathbf{s}_p \mathbf{s}_p^H = \sum_{p=1}^{P} \lambda_p \mathbf{q}_p \mathbf{q}_p^H \quad (26)$$

where $\kappa_p = \sum_{j=1}^{P} \bar{a}_{ij}$. The considered case of non-overlapped (orthogonal) components further implies that

$$\kappa_p \mathbf{s}_p = \lambda_p \mathbf{q}_p, \quad p = 1, 2, \dots, P. \quad (27)$$

*   **Partially overlapped components.** Based on (25), since the partially overlapped components are non-orthogonal; that is, such components are linearly dependent, eigenvectors can be expressed as linear combinations of such components

$$\begin{align*}
\mathbf{q}_1 &= \zeta_{11} \mathbf{s}_1 + \zeta_{21} \mathbf{s}_2 + \dots + \zeta_{P1} \mathbf{s}_P \\
\mathbf{q}_2 &= \zeta_{12} \mathbf{s}_1 + \zeta_{22} \mathbf{s}_2 + \dots + \zeta_{P2} \mathbf{s}_P \\
&\vdots \\
\mathbf{q}_M &= \zeta_{1M} \mathbf{s}_1 + \zeta_{2M} \mathbf{s}_2 + \dots + \zeta_{PM} \mathbf{s}_P,
\end{align*} \quad (28)$$

with $M = \min\{C, P\}$, i.e., for assumed $C \ge P$, $M = P$.

### 3.4. Components as the Most Concentrated Linear Combinations of Eigenvectors

Based on (28) and for assumed $M = P$, each signal component, $\mathbf{s}_p$ can be expressed as a linear combination of eigenvectors $\mathbf{q}_p$ of matrix $\mathbf{R}$, $p = 1, 2, \dots, P$; that is

$$\mathbf{s}_p = \gamma_{ip} \mathbf{q}_1 + \gamma_{2p} \mathbf{q}_2 + \dots + \gamma_{Pp} \mathbf{q}_P, \quad (29)$$

where $\gamma_{ip}$, $i = 1, 2, \dots, P$, $p = 1, 2, \dots, P$ are unknown coefficients. Obviously, there are $M = P$ linear equations for $P$ components, with $P^2$ unknown weights. Among infinitely many solutions of this undetermined system of equations, we aim at finding those combinations that produce signal components. Moreover, since components are partially overlapped, in the case when one component is detected, its contribution should be removed from all equations (linear combinations of eigenvectors) in order to avoid that it is detected again.

Obviously, for the detection of linear combinations of eigenvectors, which represent signal components, a proper detection criterion shall be established. Since non-stationary signals can be suitably represented using time-frequency representations, and signal com-
ponents tend to be concentrated along their instantaneous frequencies, our criterion will be based on time-frequency representations.

Time-frequency signal analysis provides a mathematical framework for a joint representation of signals in time and frequency domains. If $w(m)$ denotes a real-valued, symmetric window function of length $N_w$, then signal $s_p(n)$ can be represented using the STFT

$$STFT_p(n,k) = \sum_{m=0}^{N_w-1} w(m)s_p(n+m)e^{-j2\pi mk/N_w}, \quad (30)$$

which renders the frequency content of the portion of signal around the each considered instant $n$, localized by the window function $w(n)$.

To determine the level of the signal concentration in the time-frequency domain, we can exploit concentration measures. Among various approaches, inspired by the recent compressed sensing paradigm, measures based on the $\ell_p$ norm of the STFT have been used lately [18]

$$\begin{aligned} \mathcal{M}\{STFT_p(n,k)\} &= \|STFT(n,k)\|_{\rho}^{\rho} \\ &= \sum_{n} \sum_{k} |STFT(n,k)|^{\rho} = \sum_{n} \sum_{k} \mathrm{SPEC}^{\rho/2}(n,k), \end{aligned} \quad (31)$$

where $\mathrm{SPEC}(n,k) = |STFT(n,k)|^2$ represents the commonly used spectrogram, whereas $0 \le \rho \le 1$. For $\rho = 1$, the $\ell_1$-norm is obtained.

We consider $P$ components, $s_p(n), p = 1, 2, \dots, P$. Each of these components has finite support in the time-frequency domain, $\mathbb{P}_p$, with areas of support $\Pi_p, p = 1, 2, \dots, P$. Supports of partially overlapped components are also partially overlapped. Furthermore, we will make a realistic assumption that there are no components that overlap completely. Assume that $\Pi_1 \le \Pi_1 \le \dots \le \Pi_p$.

Consider further the concentration measure $\mathcal{M}\{STFT_p(n,k)\}$ of

$$\mathbf{y} = \eta_1 \mathbf{q}_1 + \eta_2 \mathbf{q}_2 + \dots + \eta_P \mathbf{q}_P, \quad (32)$$

for $p = 0$. If all components are present in this linear combination, then the concentration measure $\|STFT(n,k)\|_0$, obtained for $p = 0$ in (31), will be equal to the area of $\mathbb{P}_1 \cup \mathbb{P}_2 \cup \dots \mathbb{P}_p$.

If the coefficients $\eta_p, p = 1, 2, \dots, P$ are varied, then the minimum value of the $\ell_0$-norm based concentration measure is achieved for coefficients $\eta_1 = \gamma_{11}, \eta_2 = \gamma_{21}, \dots, \eta_P = \gamma_{P1}$ corresponding to the most concentrated signal component $s_1(n)$, with the smallest area of support, $\Pi_1$, since we have assumed, without the loss of generality, that $\Pi_1 \le \Pi_1 \le \dots \le \Pi_P$ holds. Note that, due to the calculation and sensitivity issues related with the $\ell_0$-norm, within the compressive sensing area, $\ell_1$-norm is widely used as its alternative, since under reasonable and realistic conditions, it produces the same results [31]. Therefore, it can be considered that the areas of the domains of support in this context can be measured using the $\ell_1$-norm.

The problem of extracting the first component, based on eigenvectors of the autocorrelation matrix of the input signal, can be formulated as follows

$$[\beta_{11}, \beta_{21}, \dots, \beta_{P1}] = \arg \min_{\eta_1, \dots, \eta_P} \|STFT(n,k)\|_1. \quad (33)$$

The resulting coefficients produce the first component (candidate)

$$\bar{\mathbf{s}}_1 = \beta_{11} \mathbf{q}_1 + \beta_{21} \mathbf{q}_2 + \dots + \beta_{P1} \mathbf{q}_{P1}. \quad (34)$$

Note that if $\beta_{11} = \gamma_{11}, \beta_{21} = \gamma_{21}, \dots, \beta_{P1} = \gamma_{P1}$ holds, then the component is exact; that is, $\bar{\mathbf{s}}_1 = \mathbf{s}_1$ holds. In the case when the number of signal components is larger than two, the concentration measure in (33) can have several local minima in the space of unknown coefficients $\eta_1, \eta_2, \dots, \eta_P$, corresponding not only to individual components but also to
linear combinations of two, three or more components. Depending on the minimization procedure, it can happen that the algorithm finds this local minimum; that is, a set of coefficients producing a combination of components instead of an individual component. In that case, we have not extracted successfully a component since $\bar{s}_1 \neq s_1$ in (34), but as it will be discussed next, this issue does not affect the final result, as the decomposition procedure will continue with this local minimum eliminated.

### 3.5. Extraction of Detected Component and Further Decomposition

Upon detection of the first local minimum, being a signal component or a linear combination of several components, $\bar{s}_1$, first eigenvector, $q_1$ should be replaced by $\bar{s}_1$ in the linear combination

$$ \mathbf{y} = \eta_1 \mathbf{q}_1 + \eta_2 \mathbf{q}_2 + \dots + \eta_P \mathbf{q}_P, \quad (35) $$

i.e., $q_1 = \bar{s}_1$ is further used as the first eigenvector. However, since (28) holds, the contribution of this detected component (or linear combination of components) is still present in remaining eigenvectors $q_p, p = 2, 3, \dots, P$ and shall be removed from these eigenvectors as well. To this aim, we utilize the signal deflation theory [31], and remove the projections of this component from remaining eigenvectors using

$$ \mathbf{q}_p = \frac{\mathbf{q}_p - \mathbf{q}_1^H \mathbf{q}_1}{\sqrt{1 - |\mathbf{q}_1^H \mathbf{q}_p|^2}}. \quad (36) $$

This ensures that $\bar{s}_1$ is not repeatedly detected afterward. If $\bar{s}_1 = s_1$, then the first component is found and extracted, whereas its projection on other eigenvectors is removed.

The described procedure is then repeated iteratively, with linear combination $\mathbf{y} = \eta_1 \mathbf{q}_1 + \eta_2 \mathbf{q}_2 + \dots + \eta_P \mathbf{q}_P$ with first eigenvector $\mathbf{q}_1 = \bar{s}_1$ and eigenvectors $\mathbf{q}_p, p = 1, 2, \dots, P$, modified according to (36). Upon detecting the second component (or linear combination of a small number of components), $\bar{s}_2$, the second eigenvector is replaced, $\mathbf{q}_1 = \bar{s}_2$, whereas its projections from remaining eigenvectors is removed using

$$ \mathbf{q}_p = \frac{\mathbf{q}_p - \mathbf{q}_2^H \mathbf{q}_2}{\sqrt{1 - |\mathbf{q}_2^H \mathbf{q}_p|^2}}. \quad (37) $$

The process repeats until all components are detected and extracted. These principles are incorporated into the decomposition algorithm presented in the next section.

## 4. The Decomposition Algorithm and Concentration Measure Minimization

### 4.1. Decomposition Algorithm

The decomposition procedure can be summarized with the following steps:

1. For given multivariate signal

$$ \mathbf{x}(n) = \begin{bmatrix} x^{(1)}(n) \\ x^{(2)}(n) \\ \vdots \\ x^{(C)}(n) \end{bmatrix} $$

calculate the input autocorrelation matrix

$$ \mathbf{R} = \mathbf{X}_{sen}^H \mathbf{X}_{sen} \quad (38) $$

where
$$ \mathbf{X}_{sen} = \begin{bmatrix} x^{(1)}(1) & \dots & x^{(1)}(N) \\ x^{(2)}(1) & \dots & x^{(2)}(N) \\ \vdots & \ddots & \vdots \\ x^{(C)}(1) & \dots & x^{(C)}(N) \end{bmatrix}. \qquad (39) $$

2. Find eigenvectors $\mathbf{q}_p$ and eigenvalues $\lambda_p, p = 1, 2, \dots, P$ of matrix $\mathbf{R}$.

It should be noted that the number of components, $P$, can be estimated based on the eigenvalues of matrix $\mathbf{R}$. Namely, as discussed in [31], $P$ largest eigenvalues of matrix $\mathbf{R}$ correspond to signal components. These eigenvalues are larger than the remaining $N-P$ eigenvalues. This property holds even in the presence of a high-level noise: a threshold for separation of eigenvalues corresponding to signal components can be easily determined based on the input noise variance [28].

3. Initialize variable $N_u = 0$ and variable $k = 0$. Variable $N_u$ will store the number of updates of eigenvectors $\mathbf{q}_p, p \neq i$ when projection of detected component (candidate) is removed from eigenvectors $\mathbf{q}_p, p \neq i$. Variable $k$ represents the index of the detected components.

4. For $i = 1, 2, \dots, P$, repeat the following steps:

(a) Solve minimization problem

$$ \min_{\beta_{1k}, \dots, \beta_{Pk}} \left\| \text{STFT} \left\{ \frac{1}{C} \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p \right\} \right\|_1 \quad \text{subject to } \beta_{ik} = 1 $$

where STFT{} is the STFT operator. Signal $\mathbf{y} = \frac{1}{C} \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p$ is scaled with

$$ C = \sqrt{\left\| \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p \right\|_2} $$

in order to normalize energy of the combined signal to 1. Coefficients $\beta_{1k}, \beta_{2k}, \dots, \beta_{Pk}$ are obtained as a result of the minimization.

(b) Increment component index $k \leftarrow k + 1$

(c) If for any $p \neq i, \beta_{pk} \neq 0$ holds, then

* Increment variable $N_u \leftarrow N_u + 1$

* Upon replacing the $i$-th eigenvector by the detected component,

$$ \mathbf{q}_i = \frac{1}{C} \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p \qquad (40) $$

remove projections of detected component (candidate) from remaining eigenvectors. For $l = i + 1, i + 2, \dots, P$ repeat:

i. $b = \mathbf{q}_i^H \mathbf{q}_l$

ii. $\mathbf{q}_l = \frac{1}{\sqrt{1-|b|^2}} (\mathbf{q}_l - b\mathbf{q}_i)$

5. If $N_u > 0$, return to Step 3.

Finally, as a result, we obtain the number of components, $P$, and the set of extracted components, $\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_P$.

It should be noted that checking whether $N_u > 0$ holds in Step 5 is crucial for removing possibly detected local minima of concentration measure not corresponding to individual components, but to a linear combination of several components. Namely, if this situation happens, upon applying signal deflation by removing projection of the linear combination of components from other eigenvectors, a linear dependence among eigenvectors will still remain, and it will not allow $N_u$ to fall to zero. This returns the algorithm to Step 3, and the procedure for the component detection repeats, but with the local minimum removed from
the concentration measure since all the eigenvectors are already updated in the previous cycle. Note that the component index k is reset to zero in this case.

Moreover, it should be emphasized that, while the presented procedure produces P eigenvectors, which is exactly equal to the given number of components, this number is not always a priori known. In practical applications, it can be determined based on eigenvalues of matrix R. As it will be illustrated in numerical examples, the largest eigenvalues correspond to signal components [28,31]. The remaining eigenvalues correspond to the noise. Therefore, a simple threshold can be used to calculate the exact number of signal components. Namely, we simply count the number of eigenvalues larger than a predefined threshold T, being a small positive constant. In the presence of the noise, threshold T should be at least equal to the noise variance. The larger the noise, the larger are the eigenvalues corresponding to the noise (thus, larger should the threshold T be). Of course, the procedure works without the exact information about the number of components: for p > P, eigenvectors q_p contain only the noise after the decomposition is finished.

## 4.2. Concentration Measure Minimization

The concentration measure minimization is performed in the steepest descent manner, as presented in Algorithm 1. The coefficient $\beta_{pk} = 1$ is fixed for $p = i$, whereas the values of other coefficients are varied for $\pm\Delta$. Note that real and imaginary parts are varied separately.

For the real part, and for each $p = 1, 2, \dots, P$, $p \neq i$, the $l_1$-norm based concentration measure is calculated in both cases, for auxiliary signal formed when given coefficient is increased by $\Delta$, and for the other auxiliary signal formed when $\Delta$ is subtracted from the given coefficient.

For illustration, observe linear combination $\mathbf{y} = \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p$. When $\Delta$ is added to given $\beta_{pk}$, $p \neq i$, $p = p_0$, signal

$$ \mathbf{y}_r^+ = \sum_{\substack{p=1 \\ p \neq p_0}}^{P} \beta_{pk} \mathbf{q}_p + (\beta_{p_0 k} + \Delta) \mathbf{q}_{p_0} = \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p + \Delta \mathbf{q}_{p_0} = \mathbf{y} + \Delta \mathbf{q}_{p_0} \quad (41) $$

is formed. For this signal, with energy normalized using the $l_2$-norm; that is

$$ \frac{\mathbf{y}_r^+}{\|\mathbf{y}_r^+\|_2} = \frac{\mathbf{y} + \Delta\mathbf{q}_{p_0}}{\|\mathbf{y} + \Delta\mathbf{q}_{p_0}\|_2}, \quad (42) $$

concentration measure $M_r^+$ is calculated, as the $l_1$-norm of the corresponding STFT coefficients

$$ M_r^+ = \|\text{STFT}\{\mathbf{y}_r^+\}\|_1. \quad (43) $$

Similarly, for the coefficient $\beta_{p_0 k}$ changed in the opposite direction; that is, for $-\Delta$, measure

$$ M_r^- = \|\text{STFT}\{\mathbf{y}_r^-\}\|_1. \quad (44) $$

is calculated for signal

$$ \frac{\mathbf{y}_r^-}{\|\mathbf{y}_r^-\|_2} = \frac{\mathbf{y}_r - \Delta\mathbf{q}_{p_0}}{\|\mathbf{y}_r - \Delta\mathbf{q}_{p_0}\|_2}. \quad (45) $$
**Algorithm 1 Minimization procedure**

**Input:**

* Vectors $\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_p$

* Index $i$ where corresponding vector $\mathbf{q}_i$ should be kept with unity coefficient $\beta_{pk} = 1$

* Required precision $\varepsilon$

$$1: \quad \beta_{pk} = \begin{cases} 1 & \text{for } p = i \\ 0 & \text{for } p \neq i \end{cases} \quad \text{for } p = 1, 2, \dots, P$$

$$2: \quad \mathcal{M}_{old} \leftarrow \infty$$

$$3: \quad \Delta = 0.1$$

**4: repeat**

$$5: \quad \mathbf{y} \leftarrow \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p$$

$$6: \quad \mathcal{M}_{new} \leftarrow \left\| \text{STFT} \left\{ \frac{\mathbf{y}}{\|\mathbf{y}\|_2} \right\} \right\|_1$$

**7: if** $\mathcal{M}_{new} > \mathcal{M}_{old}$ **then**

$$8: \quad \Delta \leftarrow \Delta/2$$

**9:** $\beta_{pk} \leftarrow \beta_{pk} + \nabla_p$, for $p = 1, 2, \dots, P \triangleright$ Cancel the last coefficients update

$$10: \quad \mathbf{y} \leftarrow \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p$$

**11: else**

$$12: \quad \mathcal{M}_{old} \leftarrow \mathcal{M}_{new}$$

**13: end if**

**14: for** $p = 1, 2, \dots, P$ **do**

**15: if** $p \neq i$ **then**

$$16: \quad \mathcal{M}_r^+ \leftarrow \left\| \text{STFT} \left\{ \frac{\mathbf{y} + \Delta \mathbf{q}_p}{\|\mathbf{y} + \Delta \mathbf{q}_p\|_2} \right\} \right\|_1$$

$$17: \quad \mathcal{M}_r^- \leftarrow \left\| \text{STFT} \left\{ \frac{\mathbf{y} - \Delta \mathbf{q}_p}{\|\mathbf{y} - \Delta \mathbf{q}_p\|_2} \right\} \right\|_1$$

$$18: \quad \mathcal{M}_i^+ \leftarrow \left\| \text{STFT} \left\{ \frac{\mathbf{y} + j \Delta \mathbf{q}_p}{\|\mathbf{y} + j \Delta \mathbf{q}_p\|_2} \right\} \right\|_1$$

$$19: \quad \mathcal{M}_i^- \leftarrow \left\| \text{STFT} \left\{ \frac{\mathbf{y} - j \Delta \mathbf{q}_p}{\|\mathbf{y} - j \Delta \mathbf{q}_p\|_2} \right\} \right\|_1$$

$$20: \quad \nabla_p \leftarrow 8\Delta \frac{\mathcal{M}_r^+ - \mathcal{M}_r^-}{\mathcal{M}_{new}} + j8\Delta \frac{\mathcal{M}_i^+ - \mathcal{M}_i^-}{\mathcal{M}_{new}}$$

**21: else**

$$22: \quad \nabla_p \leftarrow 0$$

**23: end if**

**24: end for**

**25:** $\beta_{pk} \leftarrow \beta_{pk} - \nabla_p$, for $p = 1, 2, \dots, P \triangleright$ Coefficients update

**26:** **until** $\sum_{p=1}^{P} |\nabla_p|^2$ is below required precision $\varepsilon$

**Output:**

* Coefficients $\beta_{1k}, \beta_{2k}, \dots, \beta_{Pk}$
Since each considered coefficient $\beta_{p0k}$ is complex-valued in general, the same procedure is repeated for the imaginary parts of coefficients. Therefore, signals

$$ \frac{\mathbf{y}_i^+}{\|\mathbf{y}_i^+\|_2} = \frac{\mathbf{y}_i + j\Delta\mathbf{q}_{p0}}{\|\mathbf{y}_i + j\Delta\mathbf{q}_{p0}\|_2}. \qquad (46) $$

and

$$ \frac{\mathbf{y}_i^-}{\|\mathbf{y}_i^-\|_2} = \frac{\mathbf{y}_i - j\Delta\mathbf{q}_{p0}}{\|\mathbf{y}_i - j\Delta\mathbf{q}_{p0}\|_2}. \qquad (47) $$

are formed, serving as a basis to calculate the corresponding concentration measures

$$ \mathcal{M}_i^+ = \|\text{STFT}\{\mathbf{y}_i^-\}\|_1 \qquad (48) $$

and

$$ \mathcal{M}_i^- = \|\text{STFT}\{\mathbf{y}_i^-\}\|_1. \qquad (49) $$

Now, based on the calculated concentration measures for variations of real and imaginary parts, concentration measure gradient $\nabla_p$ is calculated and used to determine the direction for the update of $\beta_{p0k}$

$$ \nabla_p = 8\Delta \frac{\mathcal{M}_r^+ - \mathcal{M}_r^-}{\mathcal{M}_{new}} + j8\Delta \frac{\mathcal{M}_i^+ - \mathcal{M}_i^-}{\mathcal{M}_{new}} \qquad (50) $$

where $\mathcal{M}_{new}$ used for scaling the gradient is calculated as concentration measure of

$$ \mathbf{y} = \sum_{p=1}^{P} \beta_{pk} \mathbf{q}_p, $$

scaled by its energy, before updates of coefficient $\beta_{pk}$; that is

$$ \mathcal{M}_{new} = \left\| \text{STFT} \left\{ \frac{\mathbf{y}}{\|\mathbf{y}\|_2} \right\} \right\|_1. \qquad (51) $$

For coefficients $\beta_{pk}$, $p \neq i$, the gradient is equal to zero; that is, $\nabla_p = 0$, meaning that these coefficients should not be updated.

Coefficient $\beta_{pk}$ is updated using the calculated gradient, in the steepest descent manner

$$ \beta_{pk} \leftarrow \beta_{pk} - \nabla_p, \quad \text{for } p = 1, 2, \dots, P. \qquad (52) $$

The process is repeated until $\sum_{p=1}^{P} |\nabla_p|^2$ becomes smaller than a predefined precision $\epsilon$.

## 5. Results

For the visual presentation of the results, the discrete Wigner distribution (pseudo-Wigner distribution) will be used in our numerical examples. For a discrete signal $x(n)$, this second-order time-frequency representation is calculated according to

$$ WD(n, k) = \sum_{m=0}^{N_w-1} w(m)w(-m)x(n+m)x^*(n-m)e^{-j4\pi \frac{mk}{N_w}}, \qquad (53) $$

where $w(n)$ is a window function of length $N_w$.

For examples 1, 2, and 3, the quality of the decomposition will be determined based on two criteria:

* WD calculated for the $p$th original component (signal is given analytically), denoted by $WD_p^o(n, k) = WD\{\mathbf{s}_p\}$, is compared with $WD_p^e(n, k) = WD\{\mathbf{s}_p\}$, being the WD calculated
calculated for pth extracted component, for $p = 1, 2, \dots, P$. Here, $\hat{s}$ denotes the vector of the pth extracted component, whereas $s_p$ is the actual (original) pth signal component.

* Estimation results for the discrete IFs obtained from the two previous WDs are compared by the means of mean squared error (MSE) for each pair of components. The IF estimate based on the WD of the original pth component, $WD_p^o(n, k)$, $p = 1, 2, \dots, P$ is calculated as [3]

$$k_p^o(n) = \arg\max_k WD_p^o(n, k) \quad (54)$$

whereas the IF estimate based on the WD of the pth component extracted by the proposed approach is calculated as

$$k_p^e(n) = \arg\max_k WD_p^e(n, k). \quad (55)$$

Since the extracted components do not have any particular order after the decomposition is finished, the corresponding pairs of original and extracted components are automatically determined using the following procedure:

1. For $p$ in $1, 2, \dots, P$ repeat steps (a)–(f)

(a) Calculate $k_p^o(n)$ based on (54) for analytically defined component $s_p$.

(b) Run decomposition algorithm. Use only eigenvectors corresponding to the largest $P$ eigenvalues. Each of these eigenvectors, $q_p$, contain an extracted signal component. If $P$ is not given, estimate the number of components, $P$, as the number of eigenvalues, $\lambda_p$, of matrix $R$, larger than threshold $T = \sigma_e^2 + 10^{-4}$.

(c) Initialize set $\mathbb{E} \leftarrow \emptyset$, to store the errors between IFs estimated based on the given original component $s_p$ and extracted (unordered) components $q_i$, $i = 1, 2, \dots, P$, being the outputs of the decomposition procedure.

(d) For each extracted component, $q_1, q_2, \dots, q_P$ repeat steps i–iii:

i. Calculate the IF estimate $k_i^e(n)$ as:

$$k_i^e(n) \leftarrow \arg\max_k WD\{q_i\}.$$

ii. Calculate mean squared error (MSE) between $k_p^o(n)$ and $k_i^e(n)$ as

$$MSE(i) \leftarrow \frac{1}{N} \sum_{n=0}^{N-1} |k_p^o(n) - k_i^e(n)|^2.$$

iii. $\mathbb{E} \leftarrow \mathbb{E} \cup MSE(i)$.

(e) $\hat{p} \leftarrow \arg\min_i MSE(i)$

(f) $\hat{s}_p \leftarrow q_{\hat{p}}$ is the pth estimated component, corresponding to the original component $s_p$.

Upon determining pairs of original and estimated components, $(s_p, \hat{s}_p)$, respective IF estimation MSE is calculated for each pair

$$MSE_p = \frac{1}{N} \sum_{n=0}^{N-1} |k_p^o(n) - k_p^e(n)|^2, \quad p = 1, 2, \dots, P, \quad (56)$$

where $k_p^e(n) = \arg\max_k WD\{\hat{s}_p\}$.

It should also be noted that in Examples 1–3, in order to avoid IF estimation errors at the ending edges of components (since they are characterized by time-varying amplitudes), the IF estimation is based on the WD auto-term segments larger than 10% of the maximum absolute value of the WD corresponding to the given component (auto-term), i.e.,
$$ \hat{k}_p^o(n) = \begin{cases} k_p^o(n), & \text{for } |WD_p^o(n,k)| \ge T_{WD^o}, \\ 0, & \text{for } |WD_p^o(n,k)| < T_{WD^o}, \end{cases} \qquad (57) $$

where $T_{WD^o} = 0.1 \max\{|WD_p^o(n,k)|\}$ is a threshold used to determine whether a component is present at the considered instant $n$. If it is smaller than 10% of the maximal value of the WD, it indicates that the component is not present.

### Examples

**Example 1.** To evaluate the presented theory, we consider a general form of a multicomponent signal consisted of $P$ non-stationary components

$$ x_p^{(c)}(n) = \sum_{p=1}^{P} A_p \exp\left(-\frac{n^2}{L_p^2}\right) \exp\left(j\frac{1}{N^2}n^3 + j\frac{2\pi f_p}{N}n^2 + j\frac{2\pi \phi_p}{N}n + j\theta_c\right) + \varepsilon^{(c)}(n), \quad (58) $$

$-128 \le n \le 128$ and $N = 257$. Phases $\theta_c$, $c = 1, 2, \dots, C$, are random numbers with uniform distribution drawn from interval $[-\pi, \pi]$. The signal is available in the multivariate form $\mathbf{x}(n) = [x^{(1)}(n), x^{(2)}(n), \dots, x^{(C)}(n)]^T$, and is consisted of $C$ channels, since it is embedded in a complex-valued, zero-mean noise $\varepsilon^{(c)}(n)$ with a normal distribution of its real and imaginary part, $\mathcal{N}(0, \sigma_e^2)$. Noise variance is $\sigma_e^2$, whereas $A_p = 1.2$. Parameters $f_p$ and $\phi_p$ are FM parameters, while $L_p$ is used to define the effective width of the Gaussian amplitude modulation for each component.

We generate the signal of the form (58) with $P = 6$ components, whereas the noise variance is $\sigma_e = 1$. The respective number of channels is $C = 128$. The corresponding autocorrelation matrix, $\mathbf{R}$, is calculated, according to (20), and the presented decomposition approach is used to extract the components. Eigenvalues of matrix $\mathbf{R}$ are given in Figure 2a. Largest six eigenvalues correspond to signal components, and they are clearly separable from the remaining eigenvalues corresponding to the noise. WD and spectrogram of the given signal (from one of the channels) are given in Figure 2b,c, indicating that the signal is not suitable for the classical TF analysis, since the components are highly overlapped.

Each of eigenvectors of the matrix $\mathbf{R}$ is a linear combination of components, as shown in Figure 3. The presented decomposition approach is applied to extract the components by linearly combining the eigenvectors from Figure 3. The results are shown in Figure 4a-f. Although a small residual noise is present in the extracted components, they highly match the original components, presented in Figure 4g-l. The original components in Figure 4g-l are not corrupted by the noise.

As a measure of quality, we engage MSE$_p$ given by (56), which is the error between the IF estimation result based on the $p$th extracted signal component (shown in Figure 4a-f) versus the IF estimation calculated based on the WD of original, noise-free component (from Figure 4g-l). The IF estimates and the corresponding MSEs are, for each pair of components, presented in Figure 5, for standard deviation of the noise $\sigma_e = 1$, where the number of channels is $C = 128$.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Eigenvalues of autocorrelation matrix $\mathbf{R}$, (b) Wigner distribution of the signal from Example 1 and (c) Spectrogram of the signal from Example 1. Signal consists of $P = 6$ non-stationary components. The signal is embedded in an intensive complex, Gaussian, zero-mean noise with $\sigma_e = 1$. The number of channels is $C = 128$. The largest six eigenvalues correspond to signal components.</figcaption></figure>
Since MSE<sub>p</sub> given by (56) serves as a measure of the component extraction quality, we evaluate the decomposition performance for various standard deviations of the noise, σ<sub>ε</sub> ∈ {0.1, 0.4, 0.7, 1.0, 1.3, 1.9, 2.1}. Results are presented in Table 1. The presented MSEs are calculated by averaging the results obtained based on 10 realizations of multichannel signal of the form (58) with random phases θ<sub>c</sub>, c = 1, 2, ..., C and corrupted by random realizations of the noise ε<sup>(c)</sup>(n)r, for each observed variance (standard deviation) of the noise. Based on the results from Table 1, it can be concluded that each signal component is successfully extracted for noise characterized by standard deviation up to σ<sub>ε</sub> = 1.3. For stronger noise, only some components are successfully extracted. It shall be noted that the performance of the algorithm depends also on the number of channels, C. For the results from Table 1, the number of channels was set to C = 256. A larger value of C increases the probability of successful decomposition, as investigated in [31].

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time-frequency representations of eigenvectors corresponding to the largest six eigenvalues of the autocorrelation matrix **R** of the signal from Example 1. Each eigenvector represents a linear combination of non-stationary components with polynomial frequency modulation. Panels (a–f) show Wigner distribution of each eigenvector.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** *Cont.*</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Extracted and original signal components of the non-stationary multicomponent multi-channel signal considered in Example 1. Panels (a–f) present extracted using the proposed approach, whereas panels (g–l) show Wigner distributions calculated for individual components of the original, noise-free signal.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Instantaneous frequency estimation for individual signal components based on: extracted signal components (dashed black) and original signal components (solid white). MSEs between the two IF estimates is provided for each component of the signal from Example 1. The noise variance is $\sigma_{\epsilon}^2 = 1$. Decomposition is based on $C = 128$ channels.</figcaption></figure>
<table><caption>**Table 1.** Mean squared errors (MSEs) between IF estimations based on extracted and original components, for signal from Example 1 with *P* = 6 components. *MSE*<sub>*p*</sub>, *p* = 1, 2, ..., 6 corresponds to the *p*th component. The results are presented for various values of the standard deviation of the noise, *σ*<sub>ε</sub>. The results are averaged based on 10 random realizations of signals with random phases and noise, for each considered value of *σ*<sub>ε</sub>.</caption><thead><tr><th>σ<sub>ε</sub></th><th>MSE<sub>1</sub></th><th>MSE<sub>2</sub></th><th>MSE<sub>3</sub></th><th>MSE<sub>4</sub></th><th>MSE<sub>5</sub></th><th>MSE<sub>6</sub></th></tr></thead><tbody><tr><td>0.1</td><td>-20.89 dB</td><td>-16.12 dB</td><td>-18.67 dB</td><td>-15.66 dB</td><td>-11.86 dB</td><td>-22.65 dB</td></tr><tr><td>0.4</td><td>-16.63 dB</td><td>-14.52 dB</td><td>-11.19 dB</td><td>-12.44 dB</td><td>-10.22 dB</td><td>-17.21 dB</td></tr><tr><td>0.7</td><td>-20.89 dB</td><td>-12.23 dB</td><td>-12.04 dB</td><td>-12.04 dB</td><td>-9.13 dB</td><td>-15.66 dB</td></tr><tr><td>1.0</td><td>-13.62 dB</td><td>-10.89 dB</td><td>-7.27 dB</td><td>-10.22 dB</td><td>-6.21 dB</td><td>-12.04 dB</td></tr><tr><td>1.3</td><td>-12.65 dB</td><td>-9.86 dB</td><td>-7.46 dB</td><td>-9.53 dB</td><td>-3.99 dB</td><td>-13.62 dB</td></tr><tr><td>1.6</td><td>-9.63 dB</td><td>16.67 dB</td><td>35.04 dB</td><td>-10.74 dB</td><td>27.28 dB</td><td>36.01 dB</td></tr><tr><td>1.9</td><td>-9.75 dB</td><td>32.50 dB</td><td>39.85 dB</td><td>-12.87 dB</td><td>30.28 dB</td><td>34.61 dB</td></tr><tr><td>2.1</td><td>-9.86 dB</td><td>-7.33 dB</td><td>-8.42 dB</td><td>-9.64 dB</td><td>7.03 dB</td><td>-7.46 dB</td></tr></tbody></table>

**Example 2.** The decomposition algorithm is tested on a more complex signal of the form (58), with *P* = 8 components, whereas the standard deviation of the noise is now *σ*<sub>ε</sub> = 0.1. The number of channels is *C* = 128. After the input autocorrelation matrix, **R**, is calculated, according to (20), eigendecomposition produced the eigenvalues given in Figure 6a. Signal components overlap in the time-frequency domain and, therefore, the corresponding Wigner distribution and spectrogram shown in Figure 6b,c cannot be used as adequate tools for their analysis. Figure 7 indicates that the components are neither visible in the time-frequency representation of any eigenvector corresponding to the largest eigenvalues. This is in accordance with the fact that eigenvectors contain signal components in the form of their linear combinations. Upon applying the presented multivariate decomposition procedure on this set of eigenvectors, we obtain results presented in Figure 8. By comparing the results with Wigner distributions of individual, noise-free components, shown in Figure 9, comprising the considered multicomponent signal, it can be concluded that the components are successfully extracted with preserved integrity. This is additionally confirmed by the IF estimation results shown in Figure 10, where even lower MSE values for each component can be explained by the lower noise level, as compared with results from the previous example.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** (a) Eigenvalues of autocorrelation matrix **R**, (b) Wigner distribution of the signal from Example 2, and (c) spectrogram of the signal from Example 2. Signal consists of *P* = 8 non-stationary components. The signal is embedded in an intensive complex, Gaussian, zero-mean noise with *σ*<sub>ε</sub> = 1. The number of channels is *C* = 128. The largest eight eigenvalues correspond to signal components.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** (a–h) Time-frequency representations of eigenvectors corresponding to the largest eight eigenvalues of autocorrelation matrix **R** of the signal from Example 2. Each eigenvector represents a linear combination of non-stationary components with polynomial frequency modulation.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** (a–h) Extracted signal components of the non-stationary multicomponent multichannel signal considered in Example 2. The decomposition is performed using the presented multivariate approach. The number of components is *P* = 8.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 9. (a–h) Original signal components of the non-stationary multicomponent multichannel signal considered in Example 2. Wigner distributions are calculated, each individual, noise free component.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10. Instantaneous frequency estimation for individual signal components based on the extracted signal components (dashed black) and the original signal components (solid white). MSEs between the two IF estimates is provided for each component of the signal from Example 2. The noise variance is $\sigma_{\epsilon}^2 = 0.1$. Decomposition is based on $C = 128$ channels.</figcaption></figure>
**Example 3.** To illustrate the applicability of the presented approach in decomposition of components with faster or progressive frequency variations over time, we observe a signal consisted of P = 6 components, three of which have polynomial modulations as components in model (58), whereas three other components have frequency modulations of sinusoidal nature. The first three components are defined as:

$$s_{1}^{(c)}(n) = \exp\left(-\left(\frac{n}{128}\right)^2\right) \exp\left(j40.5 \cos\left(2.34\pi n/N\right) + j10\pi n/N + j\theta_c\right) \quad (59)$$

$$s_{2}^{(c)}(n) = \begin{cases} \exp(j16 \sin(9\pi n/N) + j\theta_c), & \text{for } n \ge 0 \\ \exp(j85.33 \sin(2\pi(n+128)/N) + j\theta_c), & \text{for } n < 0 \end{cases} \quad (60)$$

$$s_{3}^{(c)}(n) = \exp\left(-\left(\frac{n}{128}\right)^2\right) \exp\left(j30.5 \sin\left(5.47\pi n/N\right) + j10\pi n/N + j\theta_c\right) \quad (61)$$

The remaining components have polynomial frequency modulation, as in previous examples:

$$s_{p}^{(c)}(n) = A_{p} \exp\left(-\frac{n^2}{L_p^2}\right) \exp\left(j\frac{1}{N^2}n^3 + j\frac{2\pi f_p}{N}n^2 + j\frac{2\pi\phi_p}{N}n + \theta_c\right) + \varepsilon^{(c)}(n) \quad (62)$$

for p = 4, 5, 6. Again, the signal is defined for discrete indices $-128 \le n \le 128$ and $N = 257$ phases $\theta_c$, $c = 1, 2, \dots, C$, are random numbers with uniform distribution drawn from interval $[-\pi, \pi]$. The resulting multicomponent signal is formed in cth channel as:

$$x_{p}^{(c)}(n) = \sum_{p=1}^{6} s_{p}^{(c)}(n) + \varepsilon^{(c)}(n), \quad (63)$$

and is, as in previous examples, embedded in additive, white, complex-valued Gaussian noise, now with variance $\sigma_{\epsilon} = 1$. The number of channels is $C = 256$. Eigenvalues of the autocorrelation matrix **R**, WD and spectrogram are given in Figure 11, again proving that the considered signal with heavily overlapped components cannot be analyzed with these tools. Eigenvectors corresponding to the largest six eigenvalues are given in Figure 12. Extracted and original components can be visually compared in Figure 13, again proving the efficiency of the approach, even in the case for components with a faster varying frequency content. This is additionally confirmed by IF estimation results in Figure 14. Larger estimation errors when faster sinusoidal frequency modulations are present are related to poorer concentration of the WD in these cases [3].

<figure><img src="image_15.png" /><figcaption>**Figure 11.** (a) Eigenvalues of autocorrelation matrix **R**, (b) Wigner distribution of the signal from Example 3, and (c) Spectrogram of the signal from Example 3. Signal consists of P = 6 non-stationary components. The signal is embedded in an intensive complex, Gaussian, zero-mean noise with $\sigma_{\epsilon} = 1$. The number of channels is C = 256. The largest six eigenvalues correspond to the signal components.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** (a–f) Time-frequency representations of eigenvectors corresponding to the largest six eigenvalues of autocorrelation matrix **R** of signal from Example 3. Each eigenvector represents a linear combination of non-stationary components with polynomial frequency modulation.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Extracted and original signal components of the non-stationary multicomponent multichannel signal considered in Example 3. Panels (a–f) present the extracted using the proposed approach, whereas panels (g–l) show Wigner distributions calculated for individual components of the original, noise-free signal.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Instantaneous frequency estimation for individual signal components based on the extracted signal components (dashed black) and the original signal components (solid white). MSEs between the two IF estimates is provided for each component of the signal from Example 3. The noise variance is $\sigma_e^2 = 1$. Decomposition is based on $C = 256$ channels.</figcaption></figure>

**Example 4.** In this example, we consider the dispersive environment setup described in Section 2.2, with the transmitter located in the water at the depth of $z_t$. However, to obtain a multivariate signal, instead of one sensor, $K = 25$ sensors are placed at the depth $z_r$, comprising the receiver at distances $r + \delta_c$, $c = 1, 2, \dots, C$, from the transmitter. Moreover, the mutual sensor distances are negligible compared with their distance from the transmitter, $r = 2000$ m; that is $\delta_c \ll r$. This further implies that the range direction remains unchanged in our model. As a response to a monochromatic signal $s(n) = \exp(j\omega_0 n)$, at sensor c, the linear combination of modes $s_p^{(c)}(n) = A_t(m, \omega_0) \exp(j\omega_0 n - jk_c(p, \omega_0)r)$ is received:

$$x^{(c)}(n) = \sum_{p=1}^{P} s_p^{(c)}(n) = \sum_{p=1}^{P} A_t(m, \omega_0) \exp(j\omega_0 n - jk_c(p, \omega_0)r), \quad (64)$$

where $c = 1, 2, \dots, C$, and wavenumbers are modeled as $k_c(m, \omega)$ [55]

$$k_r^2(m, \omega) = \left(\frac{\omega}{c}\right)^2 - \left((m - 0.5)\frac{\pi}{D + \theta_c}\right)^2, \quad (65)$$

$D = 20$ m and $\theta_c$ is a random variable drawn from interval $[-0.25, 0.25]$ with uniform distribution; therefore, corresponding to depth variations of $\pm 0.25$ cm, modeling channel depth changes due to surface waves or uneven seabed. The speed of sound propagation underwater is $c = 1500$ m/s. The same results in this example are obtained for a more precise speed, i.e., at $c = 1480$ m/s. The received multichannel signal is of the form

$$\mathbf{x}(n) = \begin{bmatrix} x^{(1)}(n) \\ x^{(2)}(n) \\ \vdots \\ x^{(C)}(n) \end{bmatrix}. \quad (66)$$

Upon performing the eigenvalue decomposition of autocorrelation matrix $\mathbf{R}$, eigenvalues shown in Figure 15a are obtained. The Wigner distribution of the received signal is shown in Figure 15b, with very close and partially overlapped nodes. Wigner distributions of individual eigenvectors are shown in Figure 16a–e. The presented procedure for the decomposition of multicomponent signals successfully extracted the individual acous-
tic modes, as presented in Figure 17a–e. Such separated acoustic modes can be further analyzed; for example, their IF can be estimated and characterized.

<figure><img src="image_19.png" /><figcaption>**Figure 15.** (a) Eigenvalues of autocorrelation matrix **R**, (b) Wigner distribution and (c) spectrogram of the considered acoustic signal from the dispersive environment. The number of sensors is *C* = 25, whereas the number of modes is *P* = 5.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** (a–e) Time-frequency representations of eigenvectors corresponding to the largest eigenvalues of autocorrelation matrix **R**. Each eigenvector represents a linear combination of acoustic modes of the signal received from the dispersive environment.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** (a–e) Extracted modes of the signal received from the dispersive acoustic environment. The decomposition is performed using the presented multivariate approach.</figcaption></figure>
**6. Discussion**

Decomposition of non-stationary multicomponent signals has been a long-term, challenging topic in time-frequency signal analysis [19–34]. Although decomposition of non-overlapping components can be done using the S-method relations with the WD [26], this approach cannot be applied when the components partially overlap, i.e., share the same domain of support in the time-frequency plane.

Other alternative methodologies are specialized for some specific signal classes, and are efficient in the case of partially overlapped components [20,25,27]. In this sense, chirplet and Radon transform-based decomposition is applicable for linear frequency modulated signals [20,25]. Inverse Radon transform has produced excellent results in separation of micro-Doppler appearing in radar signal processing, characterized by a sinusoidal frequency modulation and periodicity [27]. However, outside the scope of their predefined signal models, these techniques are inefficient in separation of non-stationary signals characterized by some different laws of non-stationarity. Another very popular concept, namely the EMD, has been also applied multivariate data [39–43]. However, successful EMD-based multicomponent signal decomposition is possible only for signals having non-overlapping components in the TF plane. Amplitude variations of components pose an additional challenge to the EMD-based decomposition. The efficiency of the proposed method does not depend on the considered frequency range, but only on the ability of a time-frequency representation to concentrate signal components in the time-frequency plane. We use the STFT in concentration measure (31) due to its ability to concentrate signal energy at the instantaneous frequency of individual signal components. The decomposition approach studied in this paper successfully extracts components highly overlapped in the time-frequency plane. The method is not sensitive to the extent of overlap of the signal components.

Since the modes appearing in the considered acoustic dispersive environment framework are characterized by a non-linear (and non-sinusoidal) law of frequency variations and have a partially overlapped support, neither of the mentioned univariate techniques can produce acceptable decomposition results.

**7. Conclusions**

Characterization of modes in the acoustic dispersive environment is an ongoing research topic. As modes are non-stationary and present in a multicomponent form in the received signals, their separation (extraction) has been a challenging task. In this paper, we have shown that the modes can be successfully extracted based on a multivariate decomposition technique that exploits the eigenanalysis of the autocorrelation matrix of the received signal. This method, which utilizes concentration measures calculated based on time-frequency representations, separates the modes while completely preserving their integrity, thus opening the possibility for their individual analysis. IF estimations based on extracted components were highly accurate, even for a high level of noise. Results indicate that the efficiency of the method is increased with the larger number of sensors (channels). Our future work will be oriented towards the analysis of the separated components. Instantaneous frequency estimation techniques developed within the time-frequency signal analysis field can be applied directly on separated modes, providing new insights and tools for the analysis of dispersive channels.

**Author Contributions:** Conceptualization, M.B. and I.S.; methodology, M.B. and I.S.; validation, M.B.; writing—original draft preparation, M.B. and I.S.; writing—review and editing, J.L., C.I., E.Z. and M.D.; visualization, I.S.; supervision, C.I. and M.D.; project administration, J.L.; funding acquisition, J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by COST action CA17137—a network for gravitational waves, geophysics, and machine learning.

**Institutional Review Board Statement:** Not applicable.
**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Boashash, B. *Time-Frequency Signal Analysis and Processing—A Comprehensive Reference*; Elsevier Science: Oxford, UK, 2003.

2. Flandrin, P. *Time-Frequency/Time-Scale Analysis*; Academic Press: San Diego, CA, USA, 1998; Volume 10.

3. Stanković, L.; Daković, M.; Thayaparan, T. *Time-Frequency Signal Analysis with Applications*: Artech House: Norwood, NJ, USA, 2013.

4. Ouahabi, A. *Signal and Image Multiresolution Analysis*; John Wiley & Sons: Hoboken, NJ, USA, 2012.

5. Boashash, B. Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals. *Proc. IEEE* **1992**, *80*, 520–538.
[CrossRef]

6. Stanković, S.; Orović, I.; Sejdić, E. *Multimedia Signals and Systems: Basic and Advance Algorithms for Signal Processing*; Springer: New York, NY, USA, 2015

7. Akan, A.; Cura, O.K. Time-frequency signal processing: Today and future. *Digit. Signal Process.* **2021**, 103216. [CrossRef]

8. Hussain, Z.M.; Boashash, B. Adaptive instantaneous frequency estimation of multicomponent FM signals using quadratic time-frequency distributions. *IEEE Trans. Signal Process.* **2002**, *50*, 1866–1876. [CrossRef]

9. Shui, P.L.; Shang, H.Y.; Zhao, Y.B. Instantaneous frequency estimation based on directionally smoothed pseudo-Wigner-Ville distribution bank. *IET Radar Sonar Navig.* **2007**, *1*, 317–325. [CrossRef]

10. Lerga, J.; Sucic, V. An Instantaneous Frequency Estimation Method Based on the Improved Sliding Pair-Wise ICI Rule. In Proceedings of the 10th International Conference on Information Science, Signal Processing and Their Applications ISSPA, Kuala Lumpur, Malaysia, 10–13 May 2010.

11. Lerga, J.; Sucic, V. Nonlinear IF Estimation Based on the Pseudo WVD Adapted Using the Improved Sliding Pairwise ICI Rule. *IEEE Signal Process. Lett.* **2009**, *16*, 953–956. [CrossRef]

12. Barkat, B.; Boashash, B. Instantaneous frequency estimation of polynomial FM signals using the peak of the PWVD: Statistical performance in the presence of additive gaussian noise. *IEEE Trans. Signal Process.* **1999**, *47*, 2480–2490. [CrossRef]

13. Sekhar, S.C.; Sreenivas, T.V. Auditory motivated level-crossing approach to instantaneous frequency estimation. *IEEE Trans. Signal Process.* **2005**, *53*, 1450–1462. [CrossRef]

14. Lerga, J.; Sucic, V.; Boashash, B. Multicomponent Noisy Signal Adaptive Instantaneous Frequency Estimation Using Components Time Support Information. *IET Signal Process.* **2014**, *8*, 277–284.

15. Lerga, J.; Sucic, V.; Boashash, B. An Efficient Algorithm for Instantaneous Frequency Estimation of Nonstationary Multicomponent Signals in Low SNR. *EURASIP J. Adv. Signal Process.* **2011**, *2011*, 16. [CrossRef]

16. Katkovnik, V.; Stanković, L. Instantaneous frequency estimation using the Wigner distribution with varying and data driven window length. *IEEE Trans. Signal Process.* **1998**, *46*, 2315–2325. [CrossRef]

17. Ivanović, V.N.; Daković, M.; Stanković, L. Performance of Quadratic Time-Frequency Distributions as Instantaneous Frequency Estimators. *IEEE Trans. Signal Process.* **2003**, *51*, 77–89. [CrossRef]

18. Stanković, L. A measure of some time-frequency distributions concentration. *Signal Process.* **2001**, *81*, 621–631. [CrossRef]

19. Stanković, L. A method for time-frequency signal analysis. *IEEE Trans. Signal Process.* **1994**, *42*, 225–229. [CrossRef]

20. Lopez-Risueno, G.; Grajal, J.; Yeste-Ojeda, O. Atomic decomposition-based radar complex signal interception. *IEEE Proc. Radar Sonar Navig.* **2003**, *150*, 323–331. [CrossRef]

21. Wei, Y.; Tan, S. Signal decomposition by the S-method with general window functions. *Signal Process.* **2012**, *92*, 288–293. [CrossRef]

22. Yang, Y.; Dong, X.; Peng, Z.; Zhang, W.; Meng, G. Component extraction for non-stationary multicomponent signal using parameterized de-chirping and band-pass filter. *IEEE Signal Process. Lett.* **2015**, *22*, 1373–1377. [CrossRef]

23. Wang, Y.; Jiang, Y. ISAR Imaging of Maneuvering Target Based on the L-Class of Fourth-Order Complex-Lag PWVD. *IEEE Trans. Geosci. Remote Sens.* **2010**, *48*, 1518–1527. [CrossRef]

24. Orović, I.; Stanković, S.; Draganić, A. Time-Frequency Analysis and Singular Value Decomposition Applied to the Highly Multicomponent Musical Signals. *Acta Acust. United Acust.* **2014**, *100*, 93–101. [CrossRef]

25. Wood, J.C.; Barry, D.T. Radon transformation of time-frequency distributions for analysis of multicomponent signals. *IEEE Trans. Signal Process.* **1994**, *42*, 3166–3177. [CrossRef]

26. Stanković, L.; Thayaparan, T.; Daković, M. Signal Decomposition by Using the S-Method with Application to the Analysis of HF Radar Signals in Sea-Clutter. *IEEE Trans. Signal Process.* **2006**, *54*, 4332–4342. [CrossRef]

27. Stanković, L.; Daković, M.; Thayaparan, T.; Popović-Bugarin, V. Inverse Radon Transform Based on Micro-Doppler Analysis from a Reduced Set of Observations. *IEEE Trans. AES* **2015**, *51*, 1155–1169. [CrossRef]

28. Stanković, L.; Mandic, D.; Daković, M.; Brajović, M. Time-frequency decomposition of multivariate multicomponent signals. *Signal Process.* **2018**, *142*, 468–479. [CrossRef]

29. Stanković, L.; Brajović, M.; Daković, M.; Mandic, D. Two-component Bivariate Signal Decomposition Based on Time-Frequency Analysis. In Proceedings of the 22nd International Conference on Digital Signal Processing IEEE DSP, London, UK, 23–25 August 2017.
30. Brajović, M.; Stanković, L.; Daković, M.; Mandic, D. Additive Noise Influence on the Bivariate Two-Component Signal Decomposition. In Proceedings of the 7th Mediterranean Conference on Embedded Computing, MECO, Budva, Montenegro, 11–14 June 2018.

31. Stanković, L.; Brajović, M.; Daković, M.; Mandic, D. On the Decomposition of Multichannel Nonstationary Multicomponent Signals. *Signal Process.* **2020**, *167*, 107261. [CrossRef]

32. Brajović, M.; Stanković, I.; Daković, M.; Mandic, D.; Stanković, L. On the Number of Channels in Multicomponent Nonstationary Noisy Signal Decomposition. In Proceedings of the 5th International Conference on Information Technology (IT 2021), Zabljak, Montenegro, 16–20 February 2021.

33. Brajović, M.; Stanković, L.; Daković, M. Decomposition of Multichannel Multicomponent Nonstationary Signals by Combining the Eigenvectors of Autocorrelation Matrix Using Genetic Algorithm. *Digit. Signal Process.* **2020**, *102*. [CrossRef]

34. Brajović, M.; Stanković, I.; Stanković, L.; Daković, M. Decomposition of Two-Component Multivariate Signals with Overlapped Domains of Support. In Proceedings of the 11th Int’l Symposium on Image and Signal Processing and Analysis (ISPA 2019), Dubrovnik, Croatia, 23–25 September 2019.

35. Ahrabian, A.; Looney, D.; Stanković, L.; Mandic, D. Synchrosqueezing-Based Time-Frequency Analysis of Multivariate Data. *Signal Process.* **2015**, *106*, 331–341. [CrossRef]

36. Lilly, J.M.; Olhede, S.C. Analysis of Modulated Multivariate Oscillations. *IEEE Trans. Signal Process.* **2012**, *60*, 600–612. [CrossRef]

37. Omidvarnia, A.; Boashash, B.; Azemi, G.; Colditz, P.; Vanhatalo, S. Generalised phase synchrony within multivariate signals: An emerging concept in time-frequency analysis. In Proceedings of the 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, Japan, 25–30 March 2012; pp. 3417–3420.

38. Lilly, J.M.; Olhede, S.C. Bivariate Instantaneous Frequency and Bandwidth. *IEEE Trans. Signal Process.* **2010**, *58*, 591–603. [CrossRef]

39. Mandic, D.P.; Rehman, N.U.; Wu, Z.; Huang, N.E. Empirical Mode Decomposition-Based Time-Frequency Analysis of Multivariate Signals: The Power of Adaptive Data Analysis. *IEEE Signal Process. Mag.* **2013**, *30*, 74–86. [CrossRef]

40. Abdullah, S.M.U.; Rehman, N.U.; Khan, M.M.; Mandic, D.P. A Multivariate Empirical Mode Decomposition Based Approach to Pansharpening. *IEEE Trans. Geosci. Remote Sens.* **2015**, *53*, 3974–3984. [CrossRef]

41. Hemakom, A.; Ahrabian, A.; Looney, D.; Rehman, N.U.; Mandic, D.P. Nonuniformly sampled trivariate empirical mode decomposition. In Proceedings of the 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), South Brisbane, QLD, Australia, 19–24 April 2015; pp. 3691–3695.

42. Wang, G.; Teng, C.; Li, K.; Zhang, Z.; Yan, X. The Removal of EOG Artifacts From EEG Signals Using Independent Component Analysis and Multivariate Empirical Mode Decomposition. *IEEE J. Biomed. Health Inform.* **2016**, *20*, 1301–1308. [CrossRef]

43. Tavildar, S.; Ashrafi, A. Application of multivariate empirical mode decomposition and canonical correlation analysis for EEG motion artifact removal. In Proceedings of the 2016 Conference on Advances in Signal Processing (CASP), Pune, India, 9–11 June 2016; pp. 150–154.

44. Omidvarnia, A.; Azemi, G.; Colditz, P.B.; Boashash, B. A time-frequency based approach for generalized phase synchrony assessment in nonstationary multivariate signals. *Digit. Signal Process.* **2013**, *23*, 780–790. [CrossRef]

45. Cobos, M.; López, J.J. Stereo audio source separation based on time-frequency masking and multilevel thresholding. *Digit. Signal Process.* **2008**, *18*, 960–976 [CrossRef]

46. Belouchrani, A.; Abed-Meraim, K.; Cardoso, J.; Moulines, E. A blind source separation technique using second-order statistics. *IEEE Trans. Signal Process.* **1997**, *45*, 434–444. [CrossRef]

47. Belouchrani, A.; Amin, M.G. Blind source separation based on time-frequency signal representations. *IEEE Trans. Signal Process.* **1998**, *46*, 2888–2897. [CrossRef]

48. Aissa-El-Bey, A.; Linh-Trung, N.; Abed-Meraim, K.; Belouchrani, A.; Grenier, Y. Underdetermined Blind Separation of Nondisjoint Sources in the Time-Frequency Domain. *IEEE Trans. Signal Process.* **2007**, *55*, 897–907. [CrossRef]

49. Liu, S.; Yu, K. Successive multivariate variational mode decomposition based on instantaneous linear mixing model. *Signal Process.* **2021**, *190*, 108311. [CrossRef]

50. Sadhu, A.; Sony, S.; Friesen, P. Evaluation of progressive damage in structures using tensor decomposition-based wavelet analysis. *J. Vib. Control.* **2019**, *25*, 2595–2610. [CrossRef]

51. Labat, V.V.; Remenieras, J.P.; Matar, O.B.; Ouahabi, A.; Patat, F. Harmonic propagation of finite amplitude sound beams: experimental determination of the nonlinearity parameter B/A. *Ultrasonics* **2000**, *38*, 292–296. [CrossRef]

52. Girault, J.M.; Kouamé, D.; Ouahabi, A.; Patat, F. Estimation of the blood Doppler frequency shift by a time-varying parametric approach. *Ultrasonics* **2000**, *38*, 682–697. [CrossRef]

53. Zhang, Y.; Amin, M.G.; Obeidat, B.A. Polarimetric Array Processing for Nonstationary Signals. In *Adaptive Antenna Arrays: Trends and Applications*; Chandran, S., Ed.; Springer: Berlin/Heidelberg, Germany, 2004; pp. 205–218.

54. Su, Z.; Ye, L. Processing of Lamb Wave Signals. In *Identification of Damage Using Lamb Waves*; Lecture Notes in Applied and Computational Mechanics; Springer: London, UK, 2009; Volume 48.

55. Ioana, C.; Jarrot, A.; Gervaise, C.; Stèphan, Y.; Quinquis, A. Localization in under water dispersive channels using the time-frequency-phase continuity of signals. *IEEE Trans. Signal Process.* **2010**, *58*, 4093–4107. [CrossRef]

56. Zhang, J.J.; Papandreou-Suppappola, A.; Gottin, B.; Ioana, C. Time-frequency characterization and receiver waveform design for shallow water environments. *IEEE Trans. Signal Process.* **2009**, *57*, 2973–2985. [CrossRef]

57. Tolstoy, I.; Clay, C.S. *Ocean Acoustics*; McGraw-Hill: New York, NY, USA, 1966; Volume 293.
58. Westwood, E.K.; Tindle, C.T.; Chapman, N.R. A normal mode model for acousto-elastic ocean environments. *J. Acoust. Soc. Am.* **1996**, *100*, 3631–3645. [CrossRef]

59. Jensen, F.B.; Kuperman, W.A.; Porter, M.B.; Schmidt, H. *Computational Ocean Acoustics*; Springer Science & Business Media: New York, NY, USA, 2011.

60. Kuperman, W.A.; Lynch, J.F. Shallow-water acoustics. *Phys. Today* **2004**, *57*, 55–61. [CrossRef]

61. Stojanovic, M. Underwater acoustic communication. In *Wiley Encyclopedia of Electrical and Electronics Engineering*; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2001.

62. Stojanovic, M.; Preisig, J. Underwater acoustic communication channels: Propagation models and statistical characterization. *IEEE Commun. Mag.* **2009**, *47*, 84–89. [CrossRef]

63. Ioana, C.; Josso, N.; Gervaise, C.; Mars, J.; Stèphan, Y. Signal analysis approach for passive tomography: Applications for dispersive channels and moving configuration. In Proceedings of the 3rd International Conference and Exhibition on Underwater Acoustic Measurements: Technologies and Results, Napflion, Greece, 21–26 June 2009.

64. de Sousa Costa, E.; Medeiros, E.B.; Filardi, J.B.C. Underwater Acoustics Modeling in Finite Depth Shallow Waters. In *Modeling and Measurement Methods for Acoustic Waves and for Acoustic Microdevices*; IntechOpen: London, UK, 2013. Available online: https://www.intechopen.com/chapters/45579 (accessed on 10 May 2021).

65. Frisk, G.V. *Ocean and Seabed Acoustics: A Theory of Wave Propagation*; Pearson Education: London, UK, 1994.

66. Zhang, J.; Papandreou-Suppappola, A. Time-frequency based waveform and receiver design for shallow water communications. In Proceedings of the 2007 IEEE International Conference on Acoustics, Speech and Signal Processing-ICASSP’07, Honolulu, HI, USA, 15–20 April 2007; Volume 3, p. III-1149.

67. Jiang, Y.; Papandreou-Suppappola, A. Discrete time-frequency characterizations of dispersive linear time-varying systems. *IEEE Trans. Signal Process.* **2007**, *55*, 2066–2076. [CrossRef]