# Spectral Analysis for Nonstationary Audio

Adrien Meynard, Bruno Torrésani

► To cite this version:

Adrien Meynard, Bruno Torrésani. Spectral Analysis for Nonstationary Audio. IEEE/ACM Transactions on Audio, Speech and Language Processing, 2018, 26 (12), pp.2371 - 2380. 10.1109/TASLP.2018.2862353. hal-01670187v3

HAL Id: hal-01670187

https://hal.science/hal-01670187v3

Submitted on 22 Aug 2018

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
# Spectral analysis for nonstationary audio

Adrien Meynard and Bruno Torrésani

## Abstract

A new approach for the analysis of nonstationary signals is proposed, with a focus on audio applications. Following earlier contributions, nonstationarity is modeled via stationarity-breaking operators acting on Gaussian stationary random signals. The focus is on time warping and amplitude modulation, and an approximate maximum-likelihood approach based on suitable approximations in the wavelet transform domain is developed. This paper provides theoretical analysis of the approximations, and introduces JEFAS, a corresponding estimation algorithm. The latter is tested and validated on synthetic as well as real audio signal.

## Index Terms

Nonstationary signals, deformation, wavelet analysis, time warping, local spectrum, Doppler effect

## I. INTRODUCTION

Nonstationarity is a key feature of acoustic signals, in particular audio signals. For example, a large part of information carried by musical and speech signals is encoded by their nonstationary nature. This is also the case for environment sounds (for example, nonstationarity of car noises or wind informs about speed variations), and many animals (e.g. bats and dolphins) use nonstationary signals for localization and communication. Beyond acoustics, amplitude and frequency modulation are of prime importance in many domains such as telecommunication.

While stationarity can be given rigorous definitions, nonstationarity is a very wide concept, as there are infinitely many ways to depart from stationarity. The theory of random signals and processes (see [1], [2] and references therein) gives a clear meaning to the notion of stationarity. In the context of time series analysis, Priestley [2], [3] was
one of the first to develop a systematic theory of nonstationary processes, introducing the class of locally stationary processes and the notion of evolutionary spectrum. A similar approach was followed in [4], who proposed a wavelet-based approach to covariance estimation for locally stationary processes (see also [5]). An alternate theory of locally stationary time series was developed by Dahlhaus [6] (see also [7] for a corresponding stationarity test). In a different context, frequency-modulated stationary signal were considered in [8], [9], and time warping models were analyzed in [10]. In several of these approaches, wavelet, time-frequency and similar representations happen to play a key role for the characterization of nonstationarity.

In a deterministic setting, a popular nonstationarity model expresses the signal as a sum of $K$ sinusoidal components $y(t) = \sum_{k=1}^{K} A_k(t) \cos(2\pi\phi_k(t))$. This model has been largely used in speech processing since early works by McAulay and Quatieri [11] (see [12] and references therein for more recent developments, and [13], [14] for probabilistic approaches). The instantaneous frequencies $\phi'_k$ of each mode give important information about the physical phenomenon. Under smoothness assumptions on functions $A_k$ and $\phi'_k$, techniques such as ridge/multiridge detection (see [15] and references therein), synchrosqueezing or reassignment have been developed to extract these quantities from a single signal observation (see [16], [17], [18] for recent accounts).

In sound processing, signals often possess a harmonic structure, which corresponds to a special case of the above model where each instantaneous frequency $\phi'_k$ is multiple of a fundamental frequency $\phi'_0$: $\phi'_k(t) = (k + 1)\phi'_0(t)$. In the special case $A_k(t) = \alpha_k A_0(t)$, we can describe such signals as a stationary signal $x(t) = \sum_{k=1}^{K} \alpha_k \cos(2\pi k t + \phi_k)$ modified by time warping and amplitude modulation: $y(t) = A_0(t)x(\phi_0(t))$. A major limit of this model is that each component is purely sinusoidal while audio signals often contain broadband information. However, sounds originating from physical phenomena can often be modeled as stationary signals, deformed by a stationarity-breaking operator (e.g. time warping, amplitude modulation). For example, sounds generated by a variable-speed engine or any stationary sound deformed by Doppler effect can be described as such. A stochastic time warping model has been introduced in [19], [20], where wavelet-based approximation and estimation techniques were developed. In [9], [21], an approximate maximum-likelihood approach was proposed for the joint estimation of the time warping and power spectrum of the underlying Gaussian stationary signal,
exploiting similar approximations.

In this paper, we build on results of [9], [21] which we extend and improve in several ways. We develop an approximate maximum likelihood method for estimating jointly time warping and amplitude modulation (not present in [9], [21]) from a single realization. While the overall structure of the algorithm is similar, we formulate the problem as a continuous parameter estimation problem, which avoids quantization effects present in the earlier approaches, and allows computing a Cramér-Rao bound for assessing the precision of the estimate. After completing the estimation, the inverse deformation can be applied to the input signal, which yields an estimate for the power spectrum.

The outline of the paper is as follows. After giving some definitions and notations in Section II, we detail in Section III the nonstationary signal models we consider, and specify the assumptions made on underlying stationary signals. We also analyze the effect of time warping and amplitude modulation in the wavelet domain, which we exploit in designing the estimation procedure. We finally propose an alternate estimation algorithm and analyze the expected performances of the corresponding estimator. Section IV is devoted to numerical results on synthetic and real signals. We also shortly describe in this section an extension published in [22] involving simultaneously time warping and frequency modulation. Mathematical developments are given as supplementary material, together with additional examples.

## II. NOTATIONS AND BACKGROUND

### A. Random signals and stationarity

Throughout this paper, we will work in the framework of the theory of random signals. Signals of interest will be modeled as realizations of random processes $t \in \mathbb{R} \rightarrow X_t \in \mathbb{C}$. Signals of interest are real-valued, however we will work with complex-valued functions since complex-valued wavelet transforms will be used. In this paper, the random processes will be denoted by uppercase letters while their realizations will be denoted by lowercase letters. The random processes will be assumed to have null mean ($\mathbb{E}\{X_t = 0\}$ for all $t$) and be second-order, i.e. they have a well-defined covariance kernel $\mathbb{E}\{X_t \overline{X}_{t'}\}$. A particularly interesting class of such stochastic processes is the class of second order (or weakly) stationary processes, for which $C_X(t - t') \triangleq \mathbb{E}\{X_t \overline{X}_{t'}\}$
is a function of $t - t'$ only. Under these assumptions, the Wiener-Khinchin theorem states that the covariance kernel may be expressed as the inverse Fourier transform of a nonnegative measure $d\eta_X$, which we will assume to be continuous with respect to the Lebesgue measure: $d\eta_X(v) = \mathcal{S}_X(v) dv$, for some nonnegative $L^1$ function $\mathcal{S}_X$ called the power spectrum. We then write

$$C_X(t) = \int \mathcal{S}_X(v) e^{2i\pi vt} dv.$$

We refer to textbooks such as [1], [2] for a more complete mathematical account of the theory, and to [21] for an extension to the setting of distribution theory.

## B. Elementary operators

Our approach rests on nonstationary models obtained by deformations of stationary random signals. We will mainly use as elementary operators the pointwise multiplication $A_\alpha$, translation $T_\tau$, dilation $D_s$, and frequency modulation $M_\nu$ defined as follows:

$$A_\alpha x(t) = \alpha x(t),$$

$$T_\tau x(t) = x(t - \tau),$$

$$D_s x(t) = q^{-s/2} x(q^{-s}t),$$

$$M_\nu x(t) = e^{2i\pi\nu t} x(t).$$

where $\alpha, \tau, s, \nu \in \mathbb{R}$ and $q > 1$ is a fixed number.

The amplitude modulation commutes with the other three operators, which satisfy the commutation rules

$$T_\tau D_s = D_s T_{q^{-s}\tau}, \quad T_\tau M_\nu = e^{-2i\pi\nu\tau} M_\nu T_\tau, \quad M_\nu D_s = D_s M_{\nu q^s}.$$

## C. Wavelet transform

Our analysis relies heavily on transforms such as the continuous wavelet transform (and discretized versions). In particular, the wavelet transform of a signal $X: t \in \mathbb{R} \rightarrow X_t$ is defined as:

$$W_X(s, \tau) \stackrel{\Delta}{=} \langle X, \psi_{s\tau} \rangle, \text{ with } \psi_{s\tau} = T_\tau D_s \psi. \quad (1)$$

where $\psi$ is the analysis wavelet, i.e. a smooth function with fast decay away from the origin. It may be shown that, for suitable choices of $\psi$, the wavelet transform is invertible (see [15]), but we will not use that property here. Notice that realizations of a continuous time random process generally do not decay at infinity. However, for
a suitably smooth and localized wavelet $\psi$, the wavelet transform can still be well defined (see [15], [21] for more details). In such a situation the wavelet transform of $X$ is a two-dimensional random field, which we analyze in the next section. Besides, in this paper the analysis wavelet $\psi$ is complex-valued and belongs to the space $H^2(\mathbb{R}) = \{\psi \in L^2(\mathbb{R}) : \text{supp}(\hat{\psi}) \subset \mathbb{R}^+\}$. In that framework, a useful property is that, if $X$ is a real, zero-mean, Gaussian process, then $W_X$ is a complex, zero-mean, circular, Gaussian random field.

Classical choices of wavelets in $H^2(\mathbb{R})$ are (analytic) derivative of Gaussian $\psi_k$ (which has $k$ vanishing moments), and the sharp wavelet $\psi_\#$ (with infinitely many vanishing moments) introduced in [22]. These can be defined in the positive Fourier domain by

$$ \hat{\psi}_k(\nu) = \nu^k e^{-k\nu^2/2\nu_0^2}, \quad \hat{\psi}_\#(\nu) = \epsilon^{\frac{\delta(\nu, \nu_0)}{\delta(\nu_1, \nu_0)}}, \quad \nu > 0 \quad (2) $$

and vanish on the negative Fourier half axis. Here $\nu_0$ is the mode of $\hat{\psi}$ and $\nu_1$ is a cutoff frequency. $\epsilon$ is a prescribed numerical tolerance chosen so that $\hat{\psi}_\#(\nu_1) = \epsilon$, and the divergence $\delta$ is defined by $\delta(a, b) = \frac{1}{2}\left(\frac{a}{b} + \frac{b}{a}\right) - 1$. The quality factor of $\psi_\#$, i.e. the center frequency to bandwidth ratio, can be expressed as $Q = 1/\sqrt{C(C+4)}$, where $C = -\delta(\nu_1, \nu_0) \ln 2 / \ln \epsilon > 0$ (see supplementary material for more details).

**Remark 1:** In (1), the scale constant $q > 1$ acts a unit selector for the scale $s$. For example, in musical terminology, $q = 2$ means that $s$ is measured in octaves, whereas for $q = 2^{1/12}$, $s$ is measured in semitones.

### D. Amplitude modulation, time warping

The nonstationary signals under consideration are obtained as linear deformations of stationary random signals. Deformations of interest here are amplitude modulations and time warpings. Amplitude modulations are pointwise multiplications by smooth functions,

$$ \mathcal{A}_a : \quad \mathcal{A}_a x(t) = a(t)x(t), \quad (3) $$

where $a \in C^1$ is a real valued function, such that

$$ 0 < c_a \le a(t) \le C_a < \infty, \quad \forall t, \quad (4) $$

for some constants $c_a, C_a \in \mathbb{R}_+^*$. Time warpings are compositions with smooth and monotonic functions,

$$ \mathcal{D}_\gamma : \quad \mathcal{D}_\gamma x(t) = \sqrt{\gamma'(t)} x(\gamma(t)), \quad (5) $$
where $\gamma \in C^2$ is a strictly increasing smooth function, satisfying the control condition [21]

$$0 < c_{\gamma} \le \gamma'(t) \le C_{\gamma} < \infty, \quad \forall t, \qquad (6)$$

for some constants $c_{\gamma}, C_{\gamma} \in \mathbb{R}_{+}^*$.

Amplitude modulations constitute a simple model for nonstationarity. While there exists a well established state of the art for demodulation algorithms in deterministic settings (in particular in telecommunications), the stochastic case has attracted less attention. In the recent literature, one may mention [10], where the so-called DEMON spectrum is proposed for amplitude modulation estimation. This problem will be tackled here using wavelet transform.

Time warping is an important transformation which has been exploited in various contexts, starting from Doppler effect, which we briefly address at the end of this paper, but also speech processing, bioacoustics (see [23] and references therein) and more generally in diverse fields such as chemistry or bioinformatics. The reference algorithm for time warping estimation is DWT (Dynamic Time Warping, see [24] for a review), which has been successfully applied to speech processing, in particular speech recognition. However, DTW is essentially a template matching algorithm, and does not address the problem considered here, where no template is available. Closer to our point of view are approaches based upon transforms such as the Harmonic transform [25] or the Fan-Chirp transform [26] (see [27] for an application to speech analysis/synthesis). These mainly involve computing a time-frequency representation of a warped copy of the input signal. Even though this does not seem to be strictly necessary, the warping function generally belongs to a parametrized family (for example quadratic), and has to be estimated. The application domain seems to be limited so far to locally harmonic, deterministic signal models, and the estimation of the warping function strongly relies on these assumptions. Wavelet transform and scalogram are also natural representations for estimating warping. Actually, interpreting (normalized) time slices of the signal scalogram as probability distributions naturally suggests to compute a time-dependent average scale, directly related to the value of the warping function. This is the approach we will use as baseline approach. However, our approximations below allow us to give a more precise meaning to this remark; in addition, scalogram lacks the phase information
which turns out to be quite relevant and yield more precise estimations.

The approach we develop below exploits complex valued wavelet transform, and combines amplitude modulation and time warping in the framework of a generic stochastic signal model, without any harmonicity assumption. This allows us to set the corresponding estimation problems as statistical inference problems, and use tools from estimation theory.

### III. JOINT ESTIMATION OF TIME WARPING AND AMPLITUDE MODULATION

#### A. Model and approximations

Let us first describe the deformation model we will mainly be using in the following.

Assume one is given a (unique) realization of a random signal of the form

$$Y = \mathcal{A}_a \mathcal{D}_\gamma X \qquad (7)$$

where $X$ is a stationary zero-mean real random process with (unknown) power spectrum $\mathcal{S}_X$. The goal is to estimate the deformation functions $a$ and $\gamma$ from this realization of $Y$, exploiting the stationarity of $X$.

**Remark 2:** The stationarity assumption is not sufficient to yield unambiguous estimates, as affine functions $\gamma(t) = \lambda t + \mu$ do not break stationarity: for any stationary $X$, $\mathcal{D}_\gamma X$ is stationary too. Therefore, the warping function $\gamma$ can only be estimated up to an affine function, as analyzed in [20], [21]. Similarly, the amplitude function $a$ can only be estimated up to a constant factor.

Key ingredients here are the smoothness of the functions $a$ and $\gamma$, and their slow variations. This allows us to perform a local analysis using smooth and localized test functions, on which the action of $\mathcal{A}_a$ and $\mathcal{D}_\gamma$ can be approximated by their so-called tangent operators $\tilde{\mathcal{A}}_a^\tau$ and $\tilde{\mathcal{D}}_\gamma^\tau$ (see [20], [9], [21], [28]). Given a test function $g$ located near $t = \tau$ (i.e. decaying fast enough as a function of $|t - \tau|$), Taylor expansions near $t = \tau$ yield

$$\mathcal{A}_a g(t) \approx \tilde{\mathcal{A}}_a^\tau g(t), \text{ with } \tilde{\mathcal{A}}_a^\tau \stackrel{\Delta}{=} A_{a(\tau)}, \qquad (8)$$

$$\mathcal{D}_\gamma g(t) \approx \tilde{\mathcal{D}}_\gamma^\tau g(t), \text{ with } \tilde{\mathcal{D}}_\gamma^\tau \stackrel{\Delta}{=} T_\tau D_{-\log_q(\gamma'(\tau))} T_{-\gamma(\tau)}. \qquad (9)$$

Therefore, we approximate the wavelet transform of $Y$ by $W_Y(s, \tau) \approx \tilde{W}_Y(s, \tau) \stackrel{\Delta}{=} \langle \tilde{\mathcal{A}}_a^\tau \tilde{\mathcal{D}}_\gamma^\tau X, T_\tau D_s \psi \rangle$, i.e.

$$\tilde{W}_Y(s, \tau) = a(\tau) W_X \left( s + \log_q(\gamma'(\tau)), \gamma(\tau) \right). \qquad (10)$$
Here, we have used the standard commutation rules of translation and dilation opera-
tors given in Section II-B.

The result below provides a quantitative assessment of the quality of the approxima-
tion. Hereafter, we denote by $||f||_\infty = \ess\sup_t|f(t)|$ the essential absolute supremum of
a function $f$.

Theorem 1: Let X be a second order zero-mean stationary random process, let Y be the nonstationary process defined in (7). Let $\psi$ be a smooth test function, localized in such a way that $|\psi(t)| \le 1/(1+|t|^\beta)$ for some $\beta > 2$. Let $W_Y$ be the wavelet transform of Y, $\tilde{W}_Y$ its approximation given in (10), and let $\varepsilon = W_Y - \tilde{W}_Y$ denote the approximation error. Assume $\psi$ and $S_X$ are such that

$$
I_X^{(\rho)} \triangleq \sqrt{\int_0^\infty \zeta^{2\rho} S_X(\zeta) \, d\zeta} < \infty, \text{ where } \rho = \frac{\beta-1}{\beta+2}.
$$

Then the approximation error $\varepsilon$ is a second order, two-dimensional complex random
field, and

$$
\mathbb{E} \left\{ |\varepsilon(s, \tau)|^2 \right\} \leq C_a^2 q^{3s} (K_1 \|\gamma''\|_{\infty} + K_2 q^{\mu s} \|\gamma''\|_{\infty}^{\rho} + K_3 \|a'\|_{\infty}^{\rho})^2
$$

where

$$
K_1 = \frac{\beta \sigma_X}{2(\beta - 2) \sqrt{c_{\gamma}}},
$$

$$
K_2 = I_X^{(\rho)} \left( \frac{\pi}{2} \right)^{\rho} \sqrt{C_{\gamma} \frac{4}{3\rho}},
$$

$$
K_3 = \frac{\sqrt{C_\gamma \beta \sigma_X}}{(\beta - 2) c_a},
$$

$$
\mu = \frac{\beta - 4}{\beta + 2},
$$

$\sigma_X^2$ being the variance of X.

The proof, which is an extension of the one given in [21], rests on Taylor approxima-
tions of $\gamma'$ and $a$ in the neighborhood of $t = \tau$ and subsequent integral majorizations,
is given as supplementary material.

Remark 3: The assumption on $\beta$ ensures that the parameters belong to the following intervals: $1/4 < \rho < 1$ and $-1/2 < \mu < 1$. Therefore, the variance of the approximation error tends to zero when the scales are small (i.e. $s \to -\infty$). Besides, the error is inversely proportional to the speed of variations of $\gamma'$ and $a$. This is consistent with the approximations of the deformation operators by their tangent operators made in equations (8) and (9).
From now on, we will assume the above approximations are valid, and work on the approximate random fields. The problem is then to estimate jointly $a, \gamma$ from $\tilde{W}_Y$, which is a zero-mean random field with covariance

$$ \mathbb{E} \{ \tilde{W}_Y(s, \tau) \overline{\tilde{W}_Y(s', \tau')} \} = C(s, s', \tau, \tau') \quad (11) $$

where the kernel $C$ reads

$$ C(s, s', \tau, \tau') = a(\tau)a(\tau')q^{\frac{s+s'}{2}} \sqrt{\gamma'(\tau)\gamma'(\tau')} \int_0^\infty \mathcal{S}_X(\xi) \\ \times \overline{\hat{\psi}}(q^s\gamma'(\tau)\xi) \hat{\psi}(q^{s'}\gamma'(\tau')\xi) e^{2i\pi\xi(\gamma(\tau)-\gamma(\tau'))} d\xi. \quad (12) $$

## B. Estimation

1) **Estimation procedure:** Our goal is to estimate both deformation functions $\gamma$ and $a$ from the approximated wavelet transform $\tilde{W}_Y$ of a realization $y$ of $Y$, assuming the latter is a reliable approximation of the true wavelet transform. From now on, we additionally assume that $X$ is a Gaussian random process. Therefore, $\tilde{W}_Y$ is a zero-mean circular Gaussian random field and its probability density function is characterized by the covariance matrix. However, equation (12) shows that besides deformation functions the covariance also depends on the power spectrum $\mathcal{S}_X$ of the underlying stationary signal $X$, which is unknown too. Therefore, the evaluation of the maximum likelihood estimate for $a$ and $\gamma$ requires a guess for $\mathcal{S}_X$. This constraint naturally brings the estimation strategy to an alternate algorithm. In [22], an estimate for $\mathcal{S}_X$ was obtained at each iteration using a Welch periodogram on a "stationarized" signal $A_{\tilde{a}}^{-1}D_{\tilde{\gamma}}^{-1}Y$, $\tilde{a}$ and $\tilde{\gamma}$ being the current estimates for the deformation functions $a$ and $\gamma$. We use here a simpler estimate, computed directly from the wavelet coefficients. The two steps of the estimation algorithm are detailed below.

**Remark 4:** The alternate likelihood maximization strategy is reminiscent of the Expectation-Maximization (EM) algorithm, the power spectrum being the nuisance parameter. However, while it would be desirable to apply directly the EM paradigm (whose convergence is proven) to our problem, the dimensionality of the latter (and the corresponding size of covariance matrices) forces us to make additional simplifications that depart from the EM scheme. Therefore we turn to a simpler approach with several dimension reduction steps.
(a) Deformation estimation. Assume that the power spectrum $\mathcal{S}_X$ is known (in fact, only an estimate $\mathcal{S}_X$ is known). Thus, we are able to write the likelihood corresponding to the observations of the wavelet coefficients. Then the maximum likelihood estimator is implemented to determine the unknown functions $\gamma$ and $a$.

The wavelet transform (1) is computed on a regular time-scale grid $\Lambda = s \times \tau$, $\delta_s$ being the scale sampling step and $F_s$ the time sampling frequency. The sizes of $s$ and $\tau$ are respectively denoted by $M_s$ and $N_\tau$.

Considering the covariance expression (12) we want to estimate the vector of parameters $\Theta = (\theta_1, \theta_2, \theta_3) \stackrel{\Delta}{=} (a(\tau)^2, \log_q(\gamma'(\tau)), \gamma(\tau))$. Let $W_y = \tilde{W}_y(\Lambda)$ denote the discretized transform and let $C_W(\Theta)$ be the corresponding covariance matrix. The related log-likelihood is

$$ \mathcal{L}(\Theta) = -\frac{1}{2} \ln |\det(C_W(\Theta))| - \frac{1}{2} C_W(\Theta)^{-1} W_y \cdot W_y. \quad (13) $$

The matrix $C_W(\Theta)$ is a matrix of size $M_s N_\tau \times M_s N_\tau$, which is generally huge. For instance, for a 5 seconds long signal, sampled at frequency $F_s = 44.1$ kHz, when the wavelet transform is computed on 8 scales, the matrix $C_W(\Theta)$ has about 3.1 trillion elements which makes it numerically intractable. In addition, due to the redundancy of the wavelet transform, $C_W(\Theta)$ turns out to be singular, and likelihood evaluation is impossible.

To overcome these issues, we use a block-diagonal regularization of the covariance matrix, obtained by forcing to zeros entries corresponding to different time indices. In other words, we disregard time correlations in the wavelet domain, which amounts to considering fixed time vector $w_{y,\tau_n} = \tilde{W}_y(s, \tau_n)$ as independent circular Gaussian vectors with zero-mean and covariance matrix

$$ C(\Theta_n)_{ij} = \theta_{n,1} C_0(\theta_{n,2})_{ij}, \quad 1 \le i, j \le M_s, \quad (14) $$

where

$$ C_0(\theta_{n,2})_{ij} = q^{(s_i+s_j)/2} \int_0^\infty \mathcal{S}_X(q^{-\theta_{n,2}}\xi) \bar{\psi}(q^{s_i}\xi) \hat{\psi}(q^{s_j}\xi) \, d\xi. \quad (15) $$

In this situation, the regularized log likelihood $\mathcal{L}^r$ splits into a sum of independent terms

$$ \mathcal{L}^r(\Theta) = \sum_n \mathcal{L}(\Theta_n), $$
where $\Theta_n = (\theta_{n,1}, \theta_{n,2}) \stackrel{\Delta}{=} (\theta_1(n), \theta_2(n))$ corresponds to the amplitude and warping parameters at fixed time $\tau_n = \tau(n)$. Notice that, in such a formalism, $\theta_{n,3} = \gamma(\tau_n)$ does not appear anymore in the covariance expression. Thus, we are led to maximize independently for each $n$

$$ \mathcal{L}(\Theta_n) = -\frac{1}{2}\ln|\det(\mathbf{C}(\Theta_n))| - \frac{1}{2}\mathbf{C}(\Theta_n)^{-1}\mathbf{w}_{y,\tau_n} \cdot \mathbf{w}_{y,\tau_n}. \quad (16) $$

For simplicity, the estimation procedure is done by an iterative algorithm (given in more details in part III-B2), which rests on two main steps. First, the log-likelihood is maximized with respect to $\theta_{n,2}$ using a gradient ascent method, for a fixed value of $\theta_{n,1}$. Second, for a fixed $\theta_{n,2}$, an estimate for $\theta_{n,1}$ is directly obtained which reads

$$ \tilde{\theta}_{n,1} = \frac{1}{M_s} \mathbf{C}_0^{-1}(\theta_{n,2}) \mathbf{w}_{y,\tau_n} \cdot \mathbf{w}_{y,\tau_n}. \quad (17) $$

(b) Spectrum estimation. Assume the amplitude modulation and time-warping parameters $\theta_1$ and $\theta_2$ are known (in fact, only estimates $\tilde{\theta}_1$ and $\tilde{\theta}_2$ are known). For any $n$ we can compute the wavelet transform

$$ \frac{1}{\theta_{n,1}^{1/2}} \tilde{W}_y(s - \theta_{n,2}, \tau_n) = W_x(s, \gamma(\tau_n)), \quad (18) $$

For fixed scale $s_m$, $\mathbf{w}_{x,s_m} \stackrel{\Delta}{=} W_x(s_m, \gamma(\tau)) \in \mathbb{C}^{N_\tau}$ is a zero-mean random circular Gaussian vector with time independent variance (as a realization of the wavelet transform of a stationary process). Hence, the empirical variance is an unbiased estimator of the variance. We then obtain the so-called wavelet spectrum

$$ \mathcal{S}_{X,\psi}(q^{-s_m}\omega_0) \stackrel{\Delta}{=} \mathbb{E} \left\{ \frac{1}{N_\tau \|\psi\|_2^2} \|\mathbf{w}_{x,s_m}\|^2 \right\} \quad (19) $$

$$ = \frac{1}{\|\psi\|_2^2} \int_0^\infty \mathcal{S}_X(\xi) q^{s_m} |\hat{\psi}(q^{s_m} \xi)|^2 d\xi, \quad (20) $$

where $\omega_0$ is the central frequency of $|\hat{\psi}|^2$. $\mathcal{S}_{X,\psi}$ is a narrowband version of $\mathcal{S}_X$ centered around frequency $v_m = q^{-s_m}\omega_0$. Besides, the bandwidth of the filter is proportional to the frequency $v_m$. This motivates the introduction of the following estimator $\tilde{\mathcal{S}}_X$ of $\mathcal{S}_X$

$$ \tilde{\mathcal{S}}_X(q^{-s_m}\omega_0) \stackrel{\Delta}{=} \frac{1}{N_\tau \|\psi\|_2^2} \|\mathbf{w}_{x,s_m}\|^2. \quad (21) $$

Finally, the estimate $\tilde{\mathcal{S}}_X$ is extended to all $\xi \in [0, F_s/2]$ by linear interpolation.
2) Algorithm: The estimation procedure is implemented in an iterative alternate optimization algorithm. This algorithm whose pseudo-code is given as Algorithm 1 is named Joint Estimation of Frequency, Amplitude, and Spectrum (JEFAS). The initialization needs an initial guess for the power spectrum $S_X$ of X. We use the spectrum estimator (21) applied to the observation Y.

After $k$ iterations of the algorithm, estimates $\tilde{\Theta}_n^{(k)}$ and $\tilde{\mathcal{F}}_X^{(k)}$ for $\Theta_n$ and $S_X$ are available. Hence we can only evaluate the plug-in estimate $\tilde{C}_0^{(k)}$ of $C_0$, obtained by replacing the power spectrum with its estimate in the covariance matrix (15). This yields an approximate expression $\mathcal{L}^{(k)}$ for the log-likelihood, which is used in place of $\mathcal{L}$ in (16) for maximum likelihood estimation. The influence of such approximations on the performances of the algorithm are discussed in section III-C.

To assess the convergence of the algorithm, the relative update of the parameters is chosen as stopping criterion:

$$ \frac{\|\tilde{\theta}_j^{(k)} - \tilde{\theta}_j^{(k-1)}\|_2^2}{\|\tilde{\theta}_j^{(k-1)}\|_2^2} < \Lambda, \text{ for } j = 1, 2, \quad (22) $$

where $0 < \Lambda < 1$ is a user defined threshold.

Finally, after convergence of the algorithm to the estimated value $\tilde{\Theta}^{(k)}$, $\log_q(\gamma')$ and $a^2$ are estimated through time by cubic spline interpolation. Besides, $\gamma$ is given by numerical integration assuming that $\gamma(0) = 0$.

**Remark 5:** To control the variances of the estimators, and the computational cost, two different discretizations of the scale axis are used for $\tilde{\theta}_1$ or $\tilde{\theta}_2$. Indeed, the computation of the log-likelihood involves the evaluation of the inverse covariance matrix. In [21], a sufficient condition for invertibility was given in the presence of noise. The major consequence induced by this condition is that when $\delta_s$ is close to zero (i.e. the sampling period of scales is small), the covariance matrix could not be numerically invertible. The scale discretization must then be sufficiently coarse to ensure good conditioning for the matrix. While this condition can be reasonably fulfilled to estimate $\theta_{n,2}$ without impairing the performances of the estimator, it cannot be applied to the estimation of $\theta_{n,1}$ because of the influence of $M_s$ on its Cramér-Rao bound (see section III-C below). The choice we made is to maximize $\mathcal{L}(\Theta_n)$ for $\theta_{n,2}$ with $w_{y,\tau_n}$ corresponding to a coarse sampling $s_p$ which is a subsampled version of the original vector $s$, the scale sampling step and the size of $s_p$ being respectively $p\delta_s$ and $\lfloor M_s/p \rfloor$ for some $p \in \mathbb{N}^*$.
**Algorithm 1 JEFAS (Joint Estimation of Frequency, Amplitude and Spectrum)**

**Initialization:** Compute an estimate $\tilde{\mathcal{F}}_Y$ of the power spectrum of $Y$ as an initial guess $\tilde{\mathcal{F}}_X^{(0)}$ for $\mathcal{S}_X$. Initialize the estimator of the squared amplitude modulation with $\tilde{\theta}_{n,1}^{(0)} = 1, \forall n$.

Compute the wavelet transform $\mathbf{W}_y$ of $y$.

$k := 1$

**while** criterion (22) is false and $k \le k_{max}$ **do**

* For each $n$, subsample $\mathbf{w}_{y,\tau_n}$ on scales $\mathbf{s}_p$, and estimate $\tilde{\theta}_{n,2}^{(k+1)}$ by maximizing the approximate log-likelihood $\mathcal{L}^{(k)}(\tilde{\theta}_{n,1}^{(k)}, \theta_{n,2})$ in (16).

* For each $n$, estimate $\tilde{\theta}_{n,1}^{(k+1)}$ by maximizing the approximate log-likelihood $\mathcal{L}^{(k)}(\theta_{n,1}, \tilde{\theta}_{n,1}^{(k+1)})$ with respect to $\theta_{n,1}$ in (16). Or, in absence of noise, directly apply equation (17) using the regularized covariance matrix given by (23).

* Construct the estimated wavelet transform $\mathbf{W}_x$ of the underlying stationary signal by interpolation from $\mathbf{W}_y$ and $\tilde{\theta}^{(k)}$ with equation (18). Estimate the corresponding power spectrum $\tilde{\mathcal{F}}_X^{(k+1)}$ with (21).

* $k := k + 1$

**end while**

* Compute $\tilde{a}$ and $\tilde{\gamma}$ by interpolation from $\tilde{\mathcal{O}}^{(k)}$.

While $\mathcal{L}(\mathcal{O}_n)$ is maximized for $\theta_{n,1}$ on the original fine sampling $\mathbf{s}$, a regularization of the covariance matrix has to be done to ensure invertibility. The regularized matrix is constructed by replacing covariance matrix $\mathbf{C}_0(\theta_{n,2})$ in (15) by its regularized version $\mathbf{C}_{0,r}(\theta_{n,2})$, given by

$$ \mathbf{C}_{0,r}(\theta_{n,2}) = (1-r)\mathbf{C}_{0,r}(\theta_{n,2}) + r\mathbf{I}, \quad (23) $$

for some regularization parameter $0 \le r \le 1$.

**Remark 6:** After convergence of the estimation algorithm, the estimated functions $\tilde{a}$ and $\tilde{\gamma}$ allow constructing a “stationarized” signal

$$ \tilde{x} = \mathcal{D}_{\tilde{\gamma}-1} \mathcal{A}_{\tilde{a}-1} y. $$

$\tilde{x}$ is an estimation of the original underlying stationary signal $x$. Furthermore, the Welch periodogram [29] may be computed from $\tilde{x}$ to obtain an estimator of $\mathcal{S}_X$ whose bias does not depend on frequency (unlike the estimator used within the iterative algorithm).
Remark 7: In order to accelerate the speed of the algorithm, the estimation can be done only on a subsampled time grid. The main effect of this choice on the algorithm concerns the final estimation of a and $\gamma$ which is more sensitive to the interpolation operation.

In the following section, we analyze quantities that enable the evaluation of the expected performances of the estimators, and their influence on the algorithm. The reader who is not directly interested in the statistical background may skip this section and jump directly to the numerical results in part IV.

### C. Performances of the estimators and the algorithm

(a) Bias. For $\theta_{n,1}$, the estimator is unbiased when the actual values of $\theta_{n,2}$ and $\mathcal{S}_X$ are known. In our case, the bias $b_{n,1}^{(k)}(\theta_{n,1}) = \mathbb{E}\{\tilde{\theta}_{n,1}^{(k)}\} - \theta_{n,1}$ is written as

$$b_{n,1}^{(k)}(\theta_{n,1}) = \frac{\theta_{n,1}}{M_s} \text{Trace} \left\{ \tilde{\mathbf{C}}_0^{(k)} (\tilde{\theta}_{n,2}^{(k)})^{-1} \mathbf{C}_0(\theta_{n,2}) - \mathbf{I} \right\}. \quad (24)$$

As expected, the better the covariance matrix estimation, the lower the bias $|b_{n,1}^{(k)}|$.

For $\theta_{n,2}$, as we do not have a closed-form expression for the estimator we are not able to give an expression of the bias. Nevertheless, if we assume that the two other true variables are known, as a maximum likelihood estimator we make sure that $\tilde{\theta}_{n,2}$ is asymptotically unbiased (i.e. $\tilde{\theta}_{n,2} \to \theta_{n,2}$ when $M_s \to \infty$).

Regarding $\mathcal{S}_X$, equation (19) shows that the estimator yields a smoothed, thus biased version of the spectrum. Proposition 1 below shows that the estimated spectrum converges to this biased version when the deformation parameters converge to their actual values.

**Proposition 1:** Let $\psi \in H^2(\mathbb{R})$ be an analytic wavelet such that $\hat{\psi}$ is bounded and $|\hat{\psi}(u)| = \mathcal{O}_{u \to \infty}(u^{-\eta})$ with $\eta > 2$. Let $\varphi_1$ and $\varphi_2$ be bounded functions defined on $\mathbb{R}^+$ by $\varphi_1(u) = u|\hat{\psi}(u)|^2$ and $\varphi_2(u) = u^2|\hat{\psi}(u)|$. Assume $\mathcal{S}_X$ is such that

$$J_X = \int_0^\infty \zeta^{-1} \mathcal{S}_X(\zeta) \, d\zeta < \infty.$$

Let $\mathcal{S}_X^{(k)}$ denote the estimation of the spectrum after $k$ iterations of the algorithm. Let $b_{\mathcal{S}_X}^{(k)}$ denote the bias defined for all $m \in [1, M_s]$ by

$$b_{\mathcal{S}_X}^{(k)}(m) = \mathbb{E} \left\{ \mathcal{S}_X^{(k)}(q^{-s_m} \omega_0) \right\} - \mathcal{S}_{X, \psi}(q^{-s_m} \omega_0).$$
Assume there exists a constant $c_{\theta_1} > 0$ such that $\theta_{n,1}^{(k)} > c_{\theta_1}, \forall n, k$. Then

$$ \|b_{\mathcal{S}_X}^{(k)}\|_{\infty} \le \frac{J_X}{\|\psi\|_2^2} \left( K_1' \|\theta_1 - \tilde{\theta}_1^{(k)}\|_{\infty} + K_2' \|\tilde{\theta}_2^{(k)} - \theta_2\|_{\infty} \right), \quad (25) $$

where

$$ K_1' = \frac{\|\varphi_1\|_{\infty}}{c_{\theta_1}}, $$

$$ K_2' = \ln(q) (\|\varphi_1\|_{\infty} + 2\|\hat{\psi}'\|_{\infty}\|\varphi_2\|_{\infty}) < \infty. $$

The proof of the Proposition is given in supplementary materials.

*Remark 8:* If $\theta_1^{(k)} \to \theta_1$ and $\theta_2^{(k)} \to \theta_2$ as $k \to \infty$, we have $\mathbb{E}\{\mathcal{S}_X^{(k)}(v_m)\} \to_k \mathcal{S}_{X,\psi}(v_m)$, as expected.

Formula (25) enables the control of the spectrum bias at frequencies $v_m = q^{-sm}\omega_0$ only. Notice also that the requirement $J_X < \infty$ forces $\mathcal{S}_X$ to vanish at zero frequency.

(b) *Variance.* The Cramér-Rao lower bound (CRLB) gives the minimum variance that can be attained by unbiased estimators. The Slepian-Bangs formula (see [30]) directly gives the following CRLB for component $\theta_{n,i}$

$$ \text{CRLB}(\theta_{n,i}) = 2 \left( \text{Trace} \left\{ \left( \mathbf{C}(\Theta_n)^{-1} \frac{\partial \mathbf{C}(\Theta_n)}{\partial \theta_{n,i}} \right)^2 \right\} \right)^{-1}. $$

This bound gives information about the variance of the estimator at convergence of the algorithm, i.e. when both $\mathcal{S}_X$ and the other parameters are well estimated.

Applying this formula to $\theta_{n,1}$ gives

$$ \mathbb{E} \left\{ \left( \tilde{\theta}_{n,1} - \mathbb{E} \{ \tilde{\theta}_{n,1} \} \right)^2 \right\} \ge \text{CRLB}(\theta_{n,1}) = \frac{2\theta_{n,1}^2}{M_s}. $$

This implies that the number of scales $M_s$ of the wavelet transform must be large enough to yield an estimator with sufficiently small variance.

For $\theta_{n,2}$, no closed-form expression is available for the CRLB. Therefore, the evaluation of this bound and its comparison with the variance of the estimator $\tilde{\theta}_{n,2}$ can only be based on numerical results, see section IV.

(c) *Robustness to noise.* Assume now observations are corrupted by a random Gaussian white noise $W$ with variance $\sigma_W^2$ (supposed to be known):

$$ Y = \mathcal{A}_a \mathcal{D}_\gamma X + W. \quad (26) $$
The estimator $\tilde{\theta}_{n,1}$ is not robust to noise. Indeed, if the maximum likelihood estimator of model (7) in the presence of such white noise, a new term $b_{n,1|W}^{(k)}(\theta_{n,1})$ must be added to the bias expression (24), which becomes

$$b_{n,1|W}^{(k)}(\theta_{n,1}) = \frac{1}{M_s} \mathrm{Trace} \left\{ \tilde{\mathbf{C}}_{0}^{(k)} (\tilde{\theta}_{n,2}^{(k)})^{-1} \mathbf{C}_{wn} \right\},$$

where $(\mathbf{C}_{wn})_{ij} = \sigma_W^2 q^{(s_i+s_j)/2} \int_0^\infty \bar{\psi}(q^{s_i}\xi) \hat{\psi}(q^{s_j}\xi) d\xi$. In practice, this term can take large values, therefore noise has to be taken into account. To do so, the covariance matrix is now written as

$$\mathbf{C}(\Theta_n)_{ij} = q^{\frac{s_i+s_j}{2}} \int_0^\infty (\theta_{n,2} \mathcal{S}_X(q^{-\theta_{n,1}}\xi) + \sigma_W^2 \bar{\psi}(q^{s_i}\xi) \hat{\psi}(q^{s_j}\xi)) d\xi \quad (27)$$

and the likelihood is modified accordingly. Formula (17) is no longer true and no closed-form expression can be derived anymore, the maximum likelihood estimate $\tilde{\theta}_{n,1}$ must be computed by a numerical scheme (here we use a simple gradient ascent).

The estimator $\tilde{\theta}_{n,2}$ is very robust to noise. Indeed, equation (27) shows that the only change in the covariance matrix formula is to replace the power spectrum $\mathcal{S}_X$ by $\mathcal{S}_Z = \mathcal{S}_X + \frac{\sigma_W^2}{\theta_{n,2}}$. The additive constant term does not impair the estimator as long as it is small in comparison with the maximum values of $\mathcal{S}_X$.

Moreover, the estimator $\tilde{\mathcal{S}}_X$ is modified because when computing $\frac{1}{\theta_{n,1}^{1/2}}\tilde{\mathcal{W}}_y(s - \theta_{n,2}, \tau_n)$ on scale $s_m$, we compute:

$$\mathbf{w}_{z,s_m} = \mathbf{w}_{x,s_m} + \mathbf{w}_{w_*,s_m},$$

where $\mathbf{w}_{w_*,s_m} = \frac{1}{\theta_1^{1/2}}\tilde{\mathcal{W}}_w(s_m - \theta_2, \tau)$ is the wavelet transform of a white noise modulated in amplitude by $a^{-1}$. Thus a constant term $\tilde{\sigma}_W$ independent of frequency is added to the new spectrum estimator $\tilde{\mathcal{S}}_Z$, so that

$$\mathbb{E}\{\tilde{\mathcal{S}}_Z\} = \mathcal{S}_{X,\psi} + \tilde{\sigma}_W^2 \quad \text{where} \quad \tilde{\sigma}_W^2 = \sigma_W^2 \frac{1}{N_\tau} \sum_{n=1}^{N_\tau} \frac{1}{\theta_{n,1}}.$$

### D. Extension: estimation of other deformations

To describe other nonstationary behaviors of audio signals, other operators can be investigated. For example, combination of time warping and frequency modulation can be considered, as was done in [22], we shortly account for this case here for the sake of completeness. Let $\alpha \in C^2$ be a smooth function, and set

$$\mathcal{M}_{\alpha}: \quad \mathcal{M}_{\alpha}x(t) = e^{2i\pi\alpha(t)}x(t), \quad (28)$$
The deformation model in [22] is of the form

$$Y = \mathcal{A}_a \mathcal{M}_\alpha \mathcal{D}_\gamma X. \qquad (29)$$

To perform joint estimation of amplitude and frequency modulation and time warping for each time, a suitable time-scale-frequency transform $\mathcal{V}$ is introduced, defined as $\mathcal{V}_X(s, v, \tau) = \langle X, \psi_{sv\tau} \rangle$, with $\psi_{sv\tau} = T_\tau M_v D_s \psi$. In that case, approximation results similar to Theorem 1 can be obtained from which the corresponding log-likelihood can be written. At fixed time $\tau$, the estimation strategy is the same as before, but the parameter space is of higher dimension, and the extra parameter $\theta_3 = \alpha'(\tau)$ complicates the log-likelihood maximization. In particular, the choice of the discretization of the two scale and frequency variables $s$ and $v$ influences performances of the estimator, in particular the Cramér-Rao bound.

## IV. NUMERICAL RESULTS

We now turn to numerical simulations and applications. A main ingredient is the choice of the wavelet transform. Here we shall always use the sharp wavelet $\psi_{\#}$ defined in (2) and set the scale constant to $q = 2$.

We systematically compare our approach to simple estimators for amplitude modulation and time warping, commonly used in applications, defined below. The approach of [20] was also implemented, but we couldn't get satisfactory results with that approach.

* Amplitude modulation: we use as baseline estimator of $a(\tau_n)^2$ the average energy $\tilde{\theta}_{n,1}^{(B)}$ defined as follows:

$$\tilde{\theta}_{n,1}^{(B)} = \frac{1}{M_s} \|\mathbf{w}_{y,\tau_n}\|^2.$$

This amounts to replace the estimated covariance matrix in (17) by the identity matrix. Notice that $\tilde{\theta}_{n,1}^{(B)}$ does not depend on the time warping estimator, and can be computed directly on the observation.

* Time warping: the baseline estimator $\tilde{\theta}_{n,2}^{(B)}$ is the scalogram scale center of mass defined as follows:

$$\tilde{\theta}_{n,2}^{(B)} = C_0 + \frac{1}{\|\mathbf{w}_{y,\tau_n}\|^2} \sum_{m=1}^{M_s} s[m] |\mathbf{w}_{y,\tau_n}[m]|^2.$$

$C_0$ is chosen such that $\tilde{\theta}_2^{(B)}$ is a zero-mean vector.
<table><caption>TABLE I
ESTIMATION MEAN SQUARE ERRORS FOR BOTH DEFORMATIONS</caption><thead><tr><th>Estimation method</th><th>Amplitude modulation</th><th>Time warping</th></tr></thead><tbody><tr><td>Baseline</td><td>0.2015</td><td>0.0232</td></tr><tr><td>JEFAS</td><td>0.0701</td><td>0.0005</td></tr></tbody></table>

Numerical evaluation is performed on both synthetic signals and deformations and real audio signals.

## A. Synthetic signal

We first evaluate the performances of the algorithm on a synthetic signal. This allows us to compare variance and bias with their theoretical values.

The simulated signal has length $N_τ = 2^{16}$ samples, sampled at $F_s = 8$ kHz (meaning the signal duration is $t_F = (N_τ - 1)/F_s \approx 8.2$ s). The spectrum $S_X$ is written as $S_X = S_1 + S_2$ where $S_l(v) = 1 + \cos(2\pi(v - v_0^{(l)})/\Delta_v^{(l)})$ if $|v - v_0^{(l)}| < \Delta_v^{(l)}/2$ and vanishes elsewhere (for $l \in \{1, 2\}$). The amplitude modulation $a$ is a sine wave $a(t) = a_0 (1 + a_1 \cos(2\pi t/T_1))$, where $a_0$ is chosen such that $t_F^{-1} \int_0^{t_F} a^2(t) dt = 1$. The time warping function $\gamma$ is such that $\log_q(\gamma'(t)) = \Gamma + \cos(2\pi t/T_2)e^{-t/T_3}$, where $\Gamma$ is chosen such that $t_F^{-1} \int_0^{t_F} \gamma'(t) dt = 1$.

JEFAS is implemented in the MATLAB/Octave scientific environment. Dimensions were set as $M_s = 106$ and $p = 7$. The wavelet transform is computed using the sharp wavelet with $\ln(\epsilon) = -25$ corresponding to a quality factor $Q = 6$. In this problem, the algorithm took 67 seconds to converge on a standard desktop computer (CPU Intel Core @ 3.20 GHz × 4, 7.7 GB RAM). Results are shown in Fig. 1 and compared with baseline estimations. For the sake of visibility, the baseline estimator of the amplitude modulation (which is very oscillatory) is not displayed, but numerical assessments are provided in Table I, which gives MSEs for the different estimations. JEFAS is clearly more precise than the baseline algorithm, furthermore its precision is well accounted for by the Cramér-Rao bound: in Fig. 1, the estimate is essentially contained within the 95 % confidence interval provided by the CRLB (assuming Gaussianity and unbiasedness).
<figure><img src="image_2.png" /><figcaption>Fig. 1. Joint amplitude modulation/time warping estimation on a synthetic signal. Top: amplitude modulation estimation ($a_1 = 0.4$ and $T_1 = t_F/3$). Bottom: warping estimation ($T_2 = t_F/2$ and $T_3 = t_F/2$).</figcaption></figure>

The left hand side of Fig. 2 displays the estimated spectrum given by formula (21). The agreement with the actual spectrum is very good, with a slight enlargement effect due to filtering by $|\hat{\psi}|^2$. The right hand side of Fig. 2 gives the evolution of the stopping criterion (22) with iterations. Numerical results show that time warping estimation converges faster than amplitude modulation estimation. Nevertheless, when fixing a stopping criterion to 0.1 % only 7 iterations are necessary for JEFAS to converge.

## B. Application to dolphin sound spectral analysis

After studying the influence of the various parameters, we now turn to real-world audio examples. First, we analyze a recording of a two seconds long dolphin vocalization sound, described in [23]. The wavelet transform of this signal in Fig. 3 shows that the warping model (7) fits well this kind of signal, except for transient clicks that are not
<figure><img src="image_3.png" /><figcaption>Fig. 2. Left: Spectrum estimation ($\nu_0^{(1)}$ = 600 Hz, $\Delta_0^{(1)}$ = 200 Hz, $\nu_0^{(2)}$ = 1.2 kHz, $\Delta_0^{(2)}$ = 400 Hz): actual (dash-dot blue line) and estimated (solid red line) spectra. Right: Relative update evolution.</figcaption></figure>

accounted for. JEFAS allows the estimation of the spectrum of the underlying stationary signal.

On the top-right of Fig. 3, we display the wavelet transform of the signal obtained by application of the inverse deformations estimated by JEFAS. Notice that the presence of clicks slightly disturbs the stationarization process. Nonetheless, it makes sense to estimate a power spectrum from this signal, since the time dependence of its wavelet transform is negligible with respect to its scale dependence. The estimated spectra from the original signal and from the estimated underlying stationary signal are displayed on the middle and the bottom of Fig. 3. Thanks to JEFAS, the harmonic structure clearly appears (bottom plot). We believe the application of JEFAS to these types of sounds can potentially bring new insights in bioacoustic applications.

### C. Application to Doppler estimation

Finally, we analyze a sound which is a recording (from a fixed location) from a racing car, moving with constant speed. The car engine sound is then deformed by the Doppler effect, which results in time warping, as explained below. Besides, as the car is moving,
<figure><img src="image_4.png" /><figcaption>Fig. 3. Dolphin sound spectral analysis. Top left: log-scalogram of the original signal. Top right: log-scalograms of unwarped and unmodulated signals. Middle: estimated spectrum from the original signal. Bottom: spectrum from the estimated underlying stationary signal.</figcaption></figure>
the closer the car to the microphone, the larger the amplitude of the recorded sound. Thus, our model fits well this signal.

The wavelet transforms of the original signal and the two estimations of the underlying stationary signal are shown in Fig. 4. While the estimation of time warping only corrects the displacement of wavelet coefficients in the time-scale domain, the joint estimation of time warping and amplitude modulation also approximately corrects nonstationary variations of the amplitudes.

The physical relevance of the estimated time warping function can be verified. Indeed, denote by $V$ the (constant) speed of the car and by $c$ the sound velocity. Fixing the time origin to the time at which the car passes in front of the observer at distance $d$, the time warping function due to Doppler effect can be shown to be

$$ \gamma'(t) = \frac{c^2}{c^2 - V^2} \left( 1 - \frac{V^2 t}{\sqrt{d^2(c^2 - V^2) + (cVt)^2}} \right). \quad (30) $$

We plot in Fig. 4 (bottom right) the estimation $\tilde{\gamma}'$ compared with its theoretical value where $d = 5$ m and $V = 54$ m/s. Clearly the estimate is close to the corresponding theoretical curve obtained with these data, which are therefore realistic values.

Nevertheless, a closer look at scalograms in Fig. 4 shows that the amplitude correction is still not perfect, due to the presence of noise, and the fact that the model remains too simple: the amplitude modulation actually depends on frequency, which is not accounted for.

## V. CONCLUSIONS

We have discussed in this paper extensions of methods and algorithms described earlier in [9], [28], [21], [22] for the joint estimation of deformation operator and power spectrum for deformed stationary signals, a problem already addressed in [20] with a different approach. Besides some improvements on the estimation algorithm itself, the main improvements described in this paper concern the following two points

1) the extension of the algorithm to the joint estimation of deformations including amplitude modulation to the model and its estimation ([22] was limited to time warping and combinations of time warping with frequency modulation, and investigated generalized wavelet transforms);
<figure><img src="image_5.png" /><figcaption>Fig. 4. Doppler estimation. Top: log-scalograms of the original (left) and unwarped and unmodulated (right) signals. Bottom left: log-scalogram of the unwarped signal. Bottom right: Estimated time warping compared with the theoretical value given in (30).</figcaption></figure>

2) a statistical study of the estimators and of the performances of JEFAS algorithm, with precise mathematical statements.

The proposed approach was validated on numerical simulations and applications to two case studies: spectral estimation from non-stationary dolphin vocalization, and Doppler estimation.

The results presented here show that the proposed extensions yield a significant improvement in terms of precision, and a better theoretical control. In particular, the continuous parameter estimation procedure avoids quantization effects that were present in [21] where the parameter space was discrete and the estimation based on exhaustive search. It also allows the derivation of precision estimates, in particular a Cramér-Rao bound. Numerical results show that the introduction of amplitude modulation also improves results. Finally, regarding the approach of [20], its domain of validity seems to be limited to small-scale (i.e. high-frequency) signals, which is not the case here.

Contrary to [20], our approach is based on (approximate) maximum likelihood estima-
tion in the Gaussian framework. Because of our choice to disregard time correlations, the estimates obtained here generally present spurious fluctuations, which can be smoothed out by appropriate filtering. A natural extension of our approach would be to introduce a smoothness prior that should avoid such filtering steps when necessary.

We believe that being able to estimate precisely warping functions can be valuable in a variety of audio applications. Speech applications have already been studied, we may also mention bioacoustic signals, for example to refine the frequency excursion indices used to assess vocal performances of songbirds (see [31] and references therein). Quite obviously, controlling warping and spectrum opens new perspectives in sound design, for example for cross-synthesis. Future work will apply the aforementioned methods to a task of blind source separation of nonstationary signals.

The code and datasets used to produce the numerical results of this paper, and other audio examples (female voice, wind, etc.) are available at the web site

https://github.com/AdMeynard/JEFAS

More details on the sharp wavelet and proofs of Theorem 1 and Proposition 1 are given as supplementary material, that also includes another case study (application to wind sound).

## ACKNOWLEDGMENT

We wish to thank M. Kowalski and R. Kronland-Martinet for fruitful discussions. We also thank the anonymous reviewers for many valuable comments and suggestions, and for bringing very interesting references to our attention. We also thank the Centre de Recherches Mathématiques of Université de Montréal, where part of this research has been done, for kind hospitality.

## REFERENCES

[1] L. H. Koopmans, *The spectral analysis of time series*, 2nd ed., ser. Probability and mathematical statistics, Z. W. Birnbaum and E. Lukacs, Eds. Elsevier, 1995, vol. 22.

[2] M. B. Priestley, *Spectral analysis and time series*, ser. Probability and mathematical statistics. Academic Press, 1982, no. v. 1-2.

[3] —, *Non-linear and non-stationary time series analysis*. Academic Press, 1988.

[4] S. Mallat, G. Papanicolaou, and Z. Zhang, “Adaptive covariance estimation of locally stationary processes,” *Ann. Statist.*, vol. 26, no. 1, pp. 1–47, 02 1998.
[5] G. P. Nason, R. von Sachs, and G. Kroisandt, "Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum," *Journal of the Royal Statistical Society. Series B (Statistical Methodology)*, vol. 62, no. 2, pp. 271-292, 2000.

[6] R. Dahlhaus, "A likelihood approximation for locally stationary processes," *Ann. Statist.*, vol. 28, no. 6, pp. 1762-1794, 12 2000.

[7] R. von Sachs and M. H. Neumann, "A wavelet-based test for stationarity," *Journal of Time Series Analysis*, vol. 21, no. 5, pp. 597-613, 2000.

[8] S. Wisdom, L. Atlas, and J. Pitton, "Extending coherence time for analysis of modulated random processes," in *2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, May 2014, pp. 340-344.

[9] H. Omer and B. Torrésani, "Estimation of frequency modulations on wideband signals; applications to audio signal analysis," in *Proceedings of the 10th International Conference on Sampling Theory and Applications (SampTA)*, G. Pfander, Ed. Eurasip Open Library, 2013, pp. 29-32.

[10] S. T. Wisdom, "Improved statistical signal processing of nonstationary random processes using time-warping," Master's thesis, School of Electrical Engineering, University of Washington, USA, March 2014.

[11] R. J. McAulay and T. F. Quatieri, "Speech analysis/synthesis based on a sinusoidal representation," *IEEE Transactions on Acoustics, Speech, Signal Processing*, vol. 34, pp. 744-754, August 1986.

[12] M. Lagrange, S. Marchand, and J. B. Rault, "Enhancing the tracking of partials for the sinusoidal modeling of polyphonic sounds," *IEEE Transactions on Audio, Speech, and Language Processing*, vol. 15, no. 5, pp. 1625-1634, July 2007.

[13] P. Depalle, G. García, and X. Rodet, "Tracking of partials for additive sound synthesis using hidden markov models," in *Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP-93)*, vol. 1, April 1993, pp. 225-228.

[14] S. Godsill and M. Davy, "Bayesian harmonic models for musical pitch estimation and analysis," in *The IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 02)*, 2002, pp. 1769-1772.

[15] R. Carmona, W. L. Hwang, and B. Torrésani, *Practical time-frequency analysis: Gabor and Wavelet Transforms With an Implementation in S*, C. K. Chui, Ed. Academic Press, 1998.

[16] F. Auger, P. Flandrin, Y.-T. Lin, S. McLaughlin, S. Meignen, T. Oberlin, and H.-T. Wu, "Time-frequency reassignment and synchrosqueezing: An overview," *IEEE Signal Processing Magazine*, vol. 30, no. 6, pp. 32-41, 2013.

[17] H.-T. Wu, "Instantaneous frequency and wave shape function (i)," *Applied and Computational Harmonic Analysis*, vol. 35, pp. 181-199, 2013.

[18] C.-Y. Lin, L. Su, and H.-T. Wu, "Wave-shape function analysis," *Journal of Fourier Analysis and Applications*, vol. 24, no. 2, pp. 451-505, Apr 2018.

[19] M. Clerc and S. Mallat, "The texture gradient equation for recovering shape from texture," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 24, no. 4, pp. 536-549, 2002.

[20] —, "Estimating deformations of stationary processes," *Ann. Statist.*, vol. 31, no. 6, pp. 1772-1821, 12 2003.

[21] H. Omer and B. Torrésani, "Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling," *Applied and Computational Harmonic Analysis*, vol. 43, no. 1, pp. 1-22, 2017.

[22] A. Meynard and B. Torrésani, "Spectral estimation for non-stationary signal classes," in *Sampling Theory and Applications*, ser. Proceedings of SampTA17, Tallinn, Estonia, Jul. 2017.
[23] D. Stowell, *Computational Analysis of Sound Scenes and Events*. Springer, 2018, ch. Computational Bioacoustic Scene Analysis, pp. 303–333.

[24] P. Senin, "Dynamic time warping algorithm review," Information and Computer Science Department, University of Hawaii at Manoa, Honolulu, USA, Tech. Rep., 2008.

[25] F. Zhang, G. Bi, and Y.-Q. Chen, "Harmonic transform," *IEE Proceedings - Vision, Image and Signal Processing*, vol. 151, no. 4, pp. 257-263, Aug 2004.

[26] L. Weruaga and M. Képesi, "The fan-chirp transform for non-stationary harmonic signals," *Signal Processing*, vol. 87, p. 15041522, 2007.

[27] R. Dunn and T. F. Quatieri, "Sinewave analysis/synthesis based on the fan-chirp transform," in *Proceedings if the 2007 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics*, 2007, pp. 247–250.

[28] H. Omer, "Modèles de déformation de processus stochastiques généralisés. application à l'estimation des non stationnarités dans les signaux audio." Ph.D. dissertation, Aix-Marseille Université, 2015.

[29] P. Welch, "The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms," *IEEE Transactions on Audio and Electroacoustics*, vol. 15, no. 2, pp. 70-73, Jun 1967.

[30] P. Stoica and R. Moses, *Introduction to Spectral Analysis*. Prentice Hall, 1997.

[31] J. Podos, D. L. Moseley, S. E. Goodwin, J. McClure, B. N. Taft, A. V. H. Strauss, C. Rega-Brodsky, and D. C. Lahti, "A fine-scale, broadly applicable index of vocal performance: frequency excursion," *Animal Behaviour*, vol. 116, pp. 203-212, 2016.