1

Spectral analysis for non-stationary audio

arXiv:1712.10252v1 [eess.AS] 29 Dec 2017

Adrien Meynard and Bruno Torrésani

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

Index Terms
Non-stationary signals, deformation, wavelet analysis, local spectrum, Doppler effect

I. I NTRODUCTION
Non-stationarity is a key feature of acoustic signals, in particular audio signals. To
mention a few examples, a large part of information carried by musical and speech
signals is encoded by their non-stationary nature, as is the case for environment sounds
(think of car noises for example, where non-stationarity informs about speed variations),
and many animals (bats, dolphins,...) use non-stationary signals for localization and
communication. Beyond acoustics, amplitude and frequency modulation are of prime
importance in telecommunication.
A. Meynard and B. Torrésani are with Aix Marseille Univ, CNRS, Centrale Marseille, I2M, UMR 7373, Marseille,
France.
Part of this work was done when the two authors were at Centre de Recherches Mathématiques, UMI 3457, CNRS
and Université de Montréal, Canada.
Manuscript received ???, 2017; revised ???, 2017.

2

While stationarity can be given rigorous definitions, non-stationarity 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 non-stationary 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 non-stationarity.
In a deterministic setting, a popular non-stationarity model expresses the signal as
a sum of K sinusoidal components y(t) = ∑K
k=1 Ak (t) cos (2πφk (t)). This model has
been largely used in speech processing since early works by McAulay and Quatieri [11]
(see also [12] and references therein for more recent developments, and [13], [14] for
probabilistic approaches). The instantaneous frequencies φk′ of each mode give important information about the physical phenomenon. Under smoothness assumptions
on functions Ak and φk′ , techniques such as ridge/multiridge detection (see [15] and
references therein), synchrosqueezing or reassignment have been developed to extract
theses quantities from a single signal observation (see [16], [17] for recent accounts).
In sound processing, signals often possess an harmonic structure, which corresponds
to a special case of the above model where each instantaneous frequency φk′ is multiple
of a fundamental frequency φ0′ : φk′ (t) = λk φ0′ (t), λk ∈ R. If the amplitudes Ak are
such that Ak (t) = αk A0 (t), we can describe such signals as a stationary signal x (t) =
∑K
k=1 α k cos(2πλk t + ϕ k ) modified by time warping and amplitude modulation: y(t) =
A0 (t) x (φ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 which have been
deformed by a stationarity breaking operator (time warping, amplitude modulation,...).
For example sounds generated by a variable speed engine or any stationary sound

3

deformed by Doppler effect can be described as such. A stochastic time warping model
has been introduced in [18], [19], where wavelet-based approximation and estimation
techniques were developed. In [9], [20], 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], [20] which we extend and improve in
several ways. We develop an approximate maximum likelihood method for estimating
jointly time warping and amplitude modulation 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 [9], [20],
and allows computing a Cramé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 plan of the paper is as follows. After giving some definitions and notations
in Section II, we detail in Section III the non-stationary signal models we consider,
and specify the assumptions made on the underlying stationary signal. 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 both synthetic signals and
real sounds. We also shortly describe in this section an extension published in [21]
involving simultaneously time warping and frequency modulation. More mathematical
developments are postponed to the Appendix.
II. N OTATIONS

AND BACKGROUND

A. Random signals, 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 processes1 t ∈

R → Xt ∈ C. In this paper, the random processes will be denoted by upper-case letters
1 Signals

of interest are real-valued, however we will also use complex-valued functions since we will use complex-

valued wavelet transforms later on.

4

while their realizations will be denoted by lower-case letters. The random processes
will be assumed to zero mean (E { Xt = 0} for all t) and be second order, i.e. they

have well defined covariance kernel E Xt X s . A particularly interesting class of such

stochastic processes is the class of second order (or weakly) stationary processes, for

∆
which CX (t − s) = E Xt X s is a function of t − s only. Under these assumptions, the

Wiener-Khinchin theorem states that the covariance kernel may be expressed as the
inverse Fourier transform of a non-negative measure dηX , which we will assume to be

continuous with respect to the Lebesgue measure: dηX (ν) = SX (ν)dν, for some nonnegative L1 function SX called the power spectrum. We then write
CX ( t ) =

Z

SX (ν)e2iπνt dν .

We refer to textbooks such as [1], [2] for a more complete mathematical account of the
theory, and to [20] for an extension to distribution theory setting.
B. Elementary operators
Our approach rests on non-stationary models obtained by deformations of stationary
random signals. We will mainly use as elementary operators the amplitude modulation
Aα , translation Tτ , dilation Ds , and frequency modulation Mν defined as follows:
Aα x (t) = αx (t) ,
Ds x (t) = q−s/2 x (q−s t) ,

Tτ x (t) = x (t − τ ) ,
Mν x (t) = e2iπνt x (t) .

where α, τ, s, ν ∈ R and q > 0 is a fixed number. The amplitude modulation commutes
with the other three operators, which satisfy the commutation rules
Tτ Ds = Ds Tq−s τ , Tτ Mν = e−2iπντ Mν Tτ , Mν Ds = Ds Mνqs .
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 ∈ R → Xt
is defined as:

W X (s, τ ) = hX, ψsτ i , with ψsτ = Tτ Ds ψ .

(1)

where ψ 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 ψ the wavelet transform is invertible

5

(see [15]), we will not use that property here. Notice that when X is a realization of
a continuous time random process, it does not need to decay at infinity. However,
for a suitably smooth and localized wavelet ψ, the wavelet transform can still be well
defined (see [15], [20] 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 ψ is complex valued and belongs to the set H 2 (R ) =

ψ ∈ L2 (R ) : supp(ψ̂ ) ⊂ R + . In that framework, a useful property is that if X is a

real zero mean Gaussian random process then, W X is a complex zero mean circular
Gaussian random field.

Classical choices of wavelets in H 2 (R ) are (analytic) derivative of Gaussian ψk (which
has k vanishing moments), and the sharp wavelet ψ♯ (with infinitely many vanishing
moments) introduced in [21]. These can be defined in the positive Fourier domain by
ψ̂k (ν) = νk e−kν

2 /2ν2
0

,

δ ( ν,ν0 )

ψ̂♯ (ν) = ǫ δ(ν1,ν0 ) , ν > 0

(2)

and vanish on the negative Fourier half axis. Here ν0 is the mode of ψ̂. In the expression
of ψ̂♯ , ν1 is chosen so that ψ̂♯ (ν1 ) = ǫ (a prescribed numerical tolerance at cutoff


frequency ν1 ), and the divergence δ is defined by δ(a, b) = 12 ba + ba − 1.
III. J OINT

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.
As said above, the non-stationary signals of interest are obtained as linear deformations
of stationary random signals. The deformations of interest here are amplitude modulations and time warpings. Amplitude modulations are pointwise multiplications by
smooth functions, defined as
Aa :

Aa x (t) = a(t) x (t) ,

(3)

where a ∈ C1 is a real valued function, such that
0 < c a ≤ a(t) ≤ Ca < ∞,

∀t ,

(4)

for some constants c a , Ca ∈ R ∗+ . Time warpings are compositions with smooth and
monotonic functions,
Dγ :

Dγ x ( t ) =

q

γ′ (t) x (γ(t))

(5)

6

where γ ∈ C2 is a strictly increasing smooth function, satisfying the control condition [20]
0 < cγ ≤ γ′ (t) ≤ Cγ < ∞,

∀t ,

(6)

for some constants cγ , Cγ ∈ R ∗+ .

Assume one is given a (unique) realization of a random signal of the form
Y = A a Dγ X

(7)

where X is a stationary zero mean real random process with (unknown) power spectrum
SX . The goal is to estimate the deformation functions a and γ from this realization of
Y, exploiting the assumed stationarity of X.
Remark 1: Clearly enough, the stationarity assumption is not sufficient to yield unambiguous estimates, since affine functions γ(t) = λt + µ do not break stationarity: for
any stationary X, Dγ X is stationary too. Therefore, the warping function γ can only be
estimated up to an affine function, as analyzed in [19] and [20]. Similarly, the amplitude
function a can only be estimated up to a constant factor.
A key ingredient of our approach is the smoothness of the deformation functions a
and γ, and their slow variations. This allows us to perform a local analysis using smooth
and localized test functions, on which the action of Aa and Dγ can be approximated
τ
g
fτ
by their so-called tangent operators A
a and Dγ (see [19], [9], [20], [22]). More precisely,

given a test function g located near t = τ (i.e. decaying fast enough as a function of

|t − τ |), Taylor expansions near t = τ yield

τ
g
fτ ∆
Aa g(t) ≈ A
a g(t) , with Aa = A a(τ ) ,

∆
fγτ g(t) , with D
fγτ =
Dγ g ( t ) ≈ D
Tτ D− logq(γ′ (τ )) T−γ(τ ) .

(8)
(9)

∆
fY (s, τ ) =
Therefore, the wavelet transform of Y will be approximated by WY (s, τ ) ≈ W
E
D
τD
τ X, T D ψ , i.e.
g
f
A
τ
s
a
γ


′
f
(10)
WY (s, τ ) = a(τ )W X s + logq (γ (τ )), γ(τ ) .

Here we have used the standard commutation rules of translation and dilation operators
given in Section II-B.

The result below provides a quantitative assessment of the quality of the approximation. There, we denote by k f k∞ = ess supt | f (t)| the essential absolute supremum of a
function f .

7

Theorem 1: Let X be a second order zero mean stationary random process, let Y be the
non-stationary process defined in (7). Let ψ be a smooth test function, localized in such
a way that |ψ(t)| ≤ 1/(1 + |t| β ) for some β > 2. Let WY denote the wavelet transform of

fY its approximation given in (10), and let ε = WY − W
fY denote the approximation
Y, W
error. Assume ψ and SX are such that
rZ
∞

(ρ) ∆

IX =

0

ξ 2ρ SX (ξ )dξ < ∞ , where ρ =

β−1
.
β+2

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

where

n
o
ρ
E |ε(s,τ )|2 ≤ Ca2 q3s K1kγ′′ k∞ + K2 qµs kγ′′ k∞ + K3 a′
βσX
√ ,
2( β − 2) c γ
p
Cγ βσX
,
K3 =
( β − 2) c a

(ρ)

K2 = I X

K1 =

µ=

 π ρ q
2

∞

Cγ

2

4
,
3ρ

β−4
,
β+2

σX2 being the variance of X.

The proof, which is an extension of the one given in [20], is given in appendix A.
Remark 2: The assumption on β ensures that the parameters belong to the following
intervals: 1/4 < ρ < 1 and −1/2 < µ < 1. Therefore, the variance of the approximation

error tends to zero when the scales are small (i.e. s → −∞). Besides, the error is
inversely proportional to the speed of variations of γ′ 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
fY , which
approximate random fields. The problem is then to estimate jointly a, γ from W
is a zero mean random field with covariance
o
n
fY (s′ , τ ′ ) = C(s, s′ , τ, τ ′ )
fY (s, τ )W
E W

(11)

where

C(s, s′, τ, τ ′ ) = a(τ )a(τ ′ )q

s+s′
2

q

Z ∞

SX (ξ )
γ′ (τ )γ′ (τ ′ )
0
  s′ ′ ′  2iπξ (γ(τ )−γ(τ ′ ))
s ′
dξ.
×ψ̂ q γ (τ )ξ ψ̂ q γ (τ )ξ e

(12)

8

B. Estimation
1) Estimation procedure: Our goal is to estimate both deformation functions γ and a
fy of a realization y of Y, assuming the
from the approximated wavelet transform W

latter is a reliable approximation of the true wavelet transform. From now on, we
fY is a zero mean
additionally assume that X is a Gaussian random process. Therefore, W
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 SX of the underlying stationary

signal X, which is unknown too. Therefore, the evaluation of the maximum likelihood
estimate for a and γ requires a guess for SX . This constraint naturally brings the
estimation strategy to an alternate algorithm. In [21], an estimate for the power spectrum
was obtained at each iteration by computing a Welch periodogram on a “stationarized”
signal Aã−1 Dγ̃−1Y, ã and γ̃ being the current estimates for the deformation functions a
and γ. We use here a simpler estimate, computed directly from the wavelet coefficients.
The two steps of the estimation algorithm are detailed below.
Remark 3: 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 nice 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 SX is known (in fact, only
an estimate 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 γ and a.
The wavelet transform (1) is computed on a regular time-scale grid Λ = s × τ, δs
being the scale sampling step and Fs the time sampling frequency. The sizes of s and τ
are respectively denoted by Ms and Nτ .
Considering the covariance expression (12) we want to estimate the vector of pa∆
fy (Λ ) denote the
rameters Θ = (θ1 , θ2 , θ3 ) = (a(τ )2 , logq (γ′ (τ )) , γ(τ )). Let Wy = W

9

discretized transform and let CW (Θ ) be the corresponding covariance matrix. The related log-likelihood is
1
1
L (Θ ) = − ln |det(CW (Θ ))| − CW (Θ )−1 Wy · Wy .
2
2

(13)

The matrix CW (Θ) is a matrix of size Ms Nτ × Ms Nτ , which is generally huge. For

instance, when the signal is a sound of 5 seconds sampled at frequency Fs = 44.1
kHz and the wavelet transform is computed on 8 scales, the matrix CW (Θ ) has about
3.1 trillion elements which makes it numerically intractable. In addition, due to the
redundancy of the wavelet transform, CW (Θ ) turns out to be singular, which makes
the evaluation of the likelihood 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
fy (s, τn ) as independent circular Gaussian vectors
consider fixed time vector wy,τn = W

with zero mean and covariance matrix

C(Θn )ij = θn,1 C0 (θn,2 )ij , 1 ≤ i, j ≤ Ms ,
where
C0 (θn,2 )ij = q(si +s j )/2

Z ∞
0

SX (q−θn,2 ξ )ψ̂ (qsi ξ )ψ̂ (qs j ξ )dξ .

(14)

(15)

In this situation, the regularized log likelihood L r splits into a sum of independent
terms
L r (Θ ) = ∑ L (Θn ) ,
n
∆

where Θn = (θn,1 , θn,2 ) = (θ1 (n), θ2 (n)) corresponds to the amplitude and warping
parameters at fixed time τn = τ (n). Notice that in such a formalism, θn,3 = γ(τn )
does not appear any more in the covariance expression. Thus we are led to maximize
independently for each n
1
1
L (Θn ) = − ln |det(C(Θn ))| − C(Θn )−1 wy,τn · wy,τn .
2
2

(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. On the one hand, the loglikelihood is maximized with respect to θn,2 using a gradient ascent method, for a fixed

10

value of θn,1 . On the other hand, for a fixed θn,2 , an estimate for θn,1 is directly obtained
which reads
θ̃n,1 =

1 −1
C (θn,2 )wy,τn · wy,τn ,
Ms 0

(17)

(b) Spectrum estimation. Assume the amplitude modulation and time-warping parameters θ1 and θ2 are known (in fact, only estimates θ̃1 and θ̃2 are known). For any n we
can compute the wavelet transform
1 f
Wy (s − θn,2 , τn ) = W x (s, γ (τn )) ,

(18)

1/2
θn,1

For fixed scale sm , w x,sm = W x (sm , γ(τ )) ∈ C Nτ 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
(
)
Z ∞
1
1
2
2
E
kwx,sm k =
SX (ξ )qsm ψ̂ (qsm ξ ) dξ
2
2
Nτ kψk2
k ψ k2 0
∆

= SX,ψ (q−sm ω0 ) ,

(19)

where ω0 is the central frequency of |ψ̂ |2 . SX,ψ is a band-pass filtered version of SX

centered around frequency νm = q−sm ω0 . Besides, the bandwidth of the filter is propor-

tional to the frequency νm . This motivates the introduction of the following estimator
S˜X of SX
S˜X (q−sm ω0 ) =

1
kwx,sm k2 .
Nτ kψk22

(20)

Finally, the estimate SX is extended to all ξ ∈ [0, Fs /2] by linear interpolation.

2) Algorithm: The estimation procedure is implemented in an iterative alternate op-

timization algorithm, whose pseudo-code is given as Algorithm 1. The initialization
requires an initial guess for the power spectrum SX of X. We use for this the spectrum
estimator (20) applied to the original observations Y.
(k)

After k iterations of the algorithm, estimates Θ̃n

(k)
and S˜X for Θn and SX are
(k)

available. Hence we can only evaluate the plug-in estimate C̃0

of C0 , obtained by

replacing the power spectrum with its estimate in the covariance matrix (15). This yields
an approximate expression L (k) for the log-likelihood, which is used in place of L
in (16) for maximum likelihood estimation. The influence of such approximations on
the performances of the algorithm are discussed in section III-C.

11

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

( k − 1) 2

θ̃j − θ̃j

2

( k − 1) 2
θ̃j
2

< T , for j = 1 and 2 ,

(21)

where 0 < T < 1 is a given threshold.
Finally, after convergence of the algorithm to the estimated value Θ̃

(k)

, logq (γ′ ) and a2

are estimated for all time by cubic spline interpolation. Besides, γ is given by numerical
integration assuming that γ(0) = 0.
Algorithm 1 Joint spectrum and deformations estimation
Initialization: Compute an estimate S˜Y of the power spectrum of Y as an initial
(0)
guess S˜X for SX . Initialize the estimator of the squared amplitude modulation with
(0)

θ̃n,1 = 1, ∀n.
Compute the wavelet transform Wy of y.
k := 1
while criterion (21) is false and k ≤ kmax do

( k + 1)

• For each n, subsample wy,τn on scales s p , and estimate θ̃n,2


(k)
approximate log-likelihood L (k) θ̃n,1 , θn,2 in (16).

by maximizing the

( k + 1)

• For each n, estimate θ̃n,1
by maximizing the approximate log-likelihood


( k + 1)
L (k) θn,1 , θ̃n,1
with respect to θn,1 in (16). Or, in absence of noise, directly apply

equation (17) using the regularized covariance matrix given by (22).

• Construct the estimated wavelet transform W x of the underlying stationary signal
(k)

by interpolation from Wy and θ̃

with equation (18). Estimate the corresponding

( k + 1)
with (20).
power spectrum S˜X

• k := k + 1
end while

• Compute ã and γ̃ by interpolation from Θ̃

(k)

.

Remark 4: In order to control the variances of the estimators, and the computational
cost, two different discretizations of the scale axis are used for θ̃1 or θ̃2 . Indeed, the
computation of the log-likelihood involves the evaluation of the inverse covariance
matrix. In [20], a sufficient condition for invertibility was given in the presence of

12

noise. The major consequence induced by this condition is that when δ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 θn,2 without impairing the performances of the estimator, it cannot be applied
to the estimation of θn,1 because of the influence of Ms on its Cramér-Rao bound (see
section III-C below). The choice we made is to maximize L (Θn ) for θn,2 with wy,τ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δs and ⌊ Ms /p⌋

for some p ∈ N ∗ . While L (Θn ) is maximized for θn,1 on the original fine sampling s,
a regularization of the covariance matrix has to be done to ensure its invertibility. The
regularized matrix is constructed by replacing covariance matrix C0 (θn,2 ) given by (15)
by its regularized version C0,r (θn,2 ), given by
C0,r (θn,2 ) = (1 − r )C0,r (θn,2 ) + rI ,

(22)

for some regularization parameter 0 ≤ r ≤ 1.
Remark 5: After convergence of the estimation algorithm, the estimated functions ã
and γ̃ allow constructing a “stationarized” signal
x̃ = Dγ̃−1 Aã−1 y .
x̃ is an estimation of the original underlying stationary signal x. Furthermore, the Welch
periodogram [23] may be computed from x̃ to obtain an estimator of SX whose bias is
not depending on frequency (unlike the estimator used within the iterative algorithm).
Remark 6: In order to accelerate the speed of the algorithm, the estimation can be done
only on a subsampled time. The main effect of this choice on the algorithm concerns
the final estimation of a and γ 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 of the algorithm. The reader
who is not directly interested in the statistical background may skip these section and
jump directly to the numerical results in part IV.

13

C. Performances of the estimators and the algorithm
(a) Bias. For θn,1 , the estimator is unbiased when the actual values of θn,2 and SX are
n
o
(k)
(k)
known. In our case, the bias bn,1 (θn,1 ) = E θ̃n,1 − θn,1 is written as




θn,1
(k) (k) −1
(k)
C0 (θn,2 ) − I .
(23)
Trace C̃0 θ̃n,2
bn,1 (θn,1 ) =
Ms
(k)

As expected, the better the covariance matrix estimation, the lower the bias bn,1 .
For θ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 θ̃n,2 is
asymptotically unbiased (i.e. θ̃n,2 → θn,2 when Ms → ∞).
Regarding SX , equation (19) shows that the estimator yields a smoothed, thus biased
version of the spectrum. Besides, 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 ψ ∈ H 2 (R ) be an analytic wavelet such that ψ̂ is bounded and

ψ̂(u) = Ou→∞ (u−η ) with η > 2. Let ϕ1 and ϕ2 be bounded functions defined on R +

by ϕ1 (u) = u ψ̂(u)

2

and ϕ2 (u) = u2 ψ̂ (u) . Assume SX is such that
JX =

(k)

Let SX

Z ∞
0

ξ −1 SX (ξ )dξ < ∞.

denote the estimation of the spectrum after k iterations of the algorithm. Let

(k)

bSX denote the bias defined for all m ∈ [[1, Ms ]] by
n
o
(k)
(k) − sm
˜
bSX (m) = E SX (q
ω0 ) − SX,ψ (q−sm ω0 ) .
(k)

Assume there exists a constant cθ1 > 0 such that θn,1 > cθ1 , ∀n, k. Then

J 
(k)
(k)
(k)
≤ X 2 K1′ θ1 − θ̃1
bSX
+ K2′ θ̃2 − θ2
,
∞
∞
∞
k ψ k2
where
K1′ =

k ϕ1 k ∞
<∞,
c θ1


K2′ = ln(q) k ϕ1 k∞ + 2kψ̂′ k∞ k ϕ2 k∞ < ∞ .

The proof is given in appendix B.

(24)

14

Remark 7: If

(k)
θ1

→ θ1 and

(k)
θ2

n

(k)
→ θ2 as k → ∞, we have E S˜X (νm )

which is the expected property.

o

→ SX,ψ (νm )
k

Formula (24) enables the control of the bias of the spectrum at frequencies νm =
q−sm ω0 only. We can also notice the required property JX < ∞ forces SX to vanish at
zero frequency.
(b) Variance. The Cramér-Rao lower bound (CRLB) gives the minimum variance that
can be attained by unbiased estimators. The Slepian-Bangs formula (see [24]) directly
gives the following CRLB for component θn,i
(
CRLB(θn,i ) = 2 Trace

C(Θn )

−1 ∂C(Θn )

∂θn,i

2 )!−1

.

This bound gives information about the variance of the estimator at convergence of the
algorithm, i.e. when both SX and the other parameters are well estimated.
Applying this formula to θn,1 gives
2
n
2θn,1
2 o

≥ CRLB(θn,1 ) =
E θ̃n,1 − E θ̃n,1
.
Ms
This implies that the number of scales Ms of the wavelet transform must be large enough
to yield an estimator with sufficiently small variance.
For θ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 θ̃n,2 can only be
based on numerical results, see section IV.
(c) Robustness to noise. Assume now the observations are corrupted by a random
2 (supposed to be known). The model becomes
Gaussian white noise W with variance σW

Y = A a Dγ X + W .

(25)

The estimator θ̃n,1 is not robust to noise. Indeed, if the maximum likelihood estimator
(k)

of model (7) in the presence of such white noise, a new term bn,1| (θn,1 ) must be added
W

to the bias expression (23), which becomes




1
(k) −1
(k)
(k)
Cwn ,
Trace C̃0 θ̃n,2
bn,1| (θn,1 ) =
W
Ms
R
2 q(si + s j )/2 ∞ ψ̂ (qsi ξ )ψ̂ (qs j ξ )dξ. In practice, this term can take large
where (Cwn )ij = σW
0

values, therefore noise has to be taken into account. To do so, the covariance matrix is
now written as
C(Θn)ij = q

si +s jZ
2

∞
0

(θn,2 SX (q−θn,1 ξ ) + σW2 )ψ̂ (qsi ξ )ψ̂ (qs j ξ )dξ ,

(26)

15

and the likelihood is modified accordingly. Formula (17) is no longer true and no closedform expression can be derived any more, the maximum likelihood estimate θ̃n,1 must
be computed by a numerical scheme (here we use a simple gradient ascent).
The estimator θ̃n,2 is very robust to noise. Indeed, equation (26) shows that the only
change in the covariance matrix formula is to replace the power spectrum SX by SZ =
σ2

W
. The additive constant term is not impairing the estimator as long as it is small
SX + θn,2

in comparison with the maximum values of SX .
Moreover, the estimator S˜X is modified because when computing
on scale sm , we compute:
wz,sm = w x,sm + ww∗ ,sm ,
1 f
Ww (sm − θ2 , τ ) is the wavelet transform of a
θ1/2
1
by a−1 . Thus a constant term σ̃W independent of

1 f
1/2 W y
θn,1

(s − θn,2 , τn )

where ww∗ ,sm =

white noise modulated

in amplitude

frequency is added to

the new spectrum estimator S˜Z , so that

2
E S˜Z = SX,ψ + σ̃W

where

2
2
σ̃W
= σW

1
Nτ

Nτ

∑

1

θ
n =1 n,1

.

D. Extension: estimation of other deformations
To describe other non-stationary 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 [21], we shortly account for this case here for the sake
of completeness. Let α ∈ C2 be a smooth function, and set
Mα :

Mα x (t) = e2iπα(t) x (t) ,

(27)

The deformation model in [21] is of the form
Y = A a M α Dγ X .

(28)

To perform joint estimation of amplitude and frequency modulation and time warping
for each time, a suitable time-scale-frequency transform V is introduced, defined as

V X (s, ν, τ ) = hX, ψsντ i, with ψsντ = Tτ Mν Ds ψ. In that case, an approximation theorem

similar to 1 can be obtained from which the corresponding log-likelihood can be written.
At fixed time τ, the strategy of estimation is the same as before, but the parameter
space is of higher dimension, and the extra parameter θ3 = α′ (τ ) complicates the loglikelihood maximization. In particular, the choice of the discretization of the two scale

16

and frequency variables s and ν influences performances of the estimator, in particular
the Cramér-Rao bound.
IV. N UMERICAL

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 ψ♯ defined
in (2) and set the scale constant q to q = 2.
We will systematically compare our approach to simple estimators for amplitude
modulation and time warping, commonly used in applications, defined below.
•

Amplitude modulation: we use as baseline estimator of a(τn )2 the average energy
(B)

θ̃n,1 defined as follows:
(B)

θ̃n,1 =

1
kwy,τn k2 .
Ms

This amounts to replace the estimated covariance matrix in (17) by the identity
(B)

matrix. Notice that θ̃n,1 does not depend on the time warping estimator, and can
be computed directly on the observation.
•

(B)

Time warping: the baseline estimator θ̃n,2 is the scalogram scale center of mass
defined as follows:
(B)

θ̃n,2 = C0 +

1

kwy,τn k2

Ms

∑

s[m]|wy,τn [m]|2 .

m=1

(B)

C0 is chosen such that θ̃n,2 is a zero mean vector.
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τ = 216 samples, sampled at Fs = 8 kHz (meaning the signal duration is t F = ( Nτ − 1)/Fs ≈ 8.2 s). The spectrum SX is written


(l )
(l )
(l )
(l )
if |ν − ν0 | < ∆ν /2
as SX = S1 + S2 where Sl (ν) = 1 + cos 2π (ν − ν0 )/∆ν

and vanishes elsewhere (for l ∈ {1, 2}). The amplitude modulation a is a sine wave
R
1 tF 2
a(t) = a0 (1 + a1 cos(2πt/T1 )), where a0 is chosen such that t−
F
0 a (t)dt = 1. The

17

3
2.5
2
1.5
1
0.5
0

0

1

2

3

4

5

6

7

8

5

6

7

8

Time (s)
1

0.5

0

-0.5

0

1

2

3

4

Time (s)

Fig. 1. Joint amplitude modulation/time warping estimation on a synthetic signal. Top: amplitude modulation
estimation (a1 = 0.4 and T1 = t F /3). Bottom: warping estimation (T2 = t F /2 and T3 = t F /2).

Estimation

Amplitude

Time

method

modulation

warping

Baseline

2.01 × 10−1

2.32 × 10−2

Algorithm 1

7.01 × 10−2

4.91 × 10−4

TABLE I
M EAN SQUARE ERROR OF THE ESTIMATION METHODS FOR BOTH DEFORMATIONS

time warping function γ is such that logq (γ′ (t)) = Γ + cos(2πt/T2 )e−t/T3 , where Γ is
R
1 tF ′
chosen such that t−
F
0 γ (t)dt = 1.
The estimation algorithm is implemented with Ms = 106 and p = 7. 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. The proposed algorithm is clearly more precise than the baseline
algorithm, furthermore its precision is well accounted for by the Cramér-Rao bound: in
Fig. 1, the estimate is essentially contained within the 95 % confidence interval provided
by the CRLB (assuming gaussianity and unbiasedness).
The left hand side of Fig. 2 displays the estimated spectrum given by the algorithm
with formula (20). The agreement with the actual spectrum is very good, with a slight
enlargement effect due to to the filtering by |ψ̂|2 . The right hand side of Fig. 2 gives

18

10 0

2
1.8

10 -1

1.6

Relative update

Spectrum

1.4
1.2
1
0.8
0.6

10 -2

10

-3

10

-4

0.4
0.2
0

400

600

800

1000

1200

1400

1600

1

T =0.1 %
Time warping
Amplitude modulation
Spectrum
2

3

Frequency (Hz)

( 1)

Fig. 2. Left: Spectrum estimation (ν0

4

5

6

7

Iteration

( 1)

= 600 Hz, ∆ ν

( 2)

= 200 Hz, ν0

( 2)

= 1.2 kHz, ∆ ν

= 400 Hz): actual (dash-dot

blue line) and estimated (solid red line) spectra. Right: Stopping criterion evolution.

the evolution of the stopping criterion (21) 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 the algorithm to converge.

B. Application to Doppler estimation
After studying the influence of the various parameters, let us now turn to a real life
audio example. The analyzed sound 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 the sound emitted by the car. Besides, as the car is
moving, 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. 3. 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

19

5513
2756
1378
689
345
172
0

0.5

1

1.5

Time (s)

2

0

0.5

1

1.5

2

Time (s)

1.15
1.1
1.05
1
0.95
0.9
0.85
0.2

0.4

0.6

0.8

1

Fig. 3. Doppler estimation. Top left: Scalogram of the original signal. Top right: Scalogram of the unwarped and
unmodulated signal. Bottom left: Scalogram of the unwarped signal. Bottom right: Estimated time warping compared
with the theoretical value given in (29).

warping function due to Doppler effect can be shown to be
c2
γ′ (t) = 2
c − V2

1− p

V 2t
d2 (c2 − V 2 ) + (cVt)2

!

(29)
.

We plot in Fig. 3 (bottom right) the estimation γ̃′ 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. 3 shows that the amplitude correction
is not perfect, due to the presence of noise, and the fact that the model is still too simple:
the amplitude modulation actually depends on frequency, which is not accounted for.

V. C ONCLUSIONS
We have discussed in this paper extensions of methods and algorithms described
earlier in [9], [22], [20], [21] for the joint estimation of deformation operator and power
spectrum for deformed stationary signals, a problem already addressed in [19] with
a different approach. 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;
2) a statistical study of the estimators and of the performances of the algorithm.

20

The proposed approach was validated on numerical simulations and an application to
Doppler estimation.
The results presented here show that the proposed extensions result in 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 [20] where the parameter space was discrete and the estimation based on exhaustive
search. Regarding the approach of [19], its domain of validity seems to be limited to
small scales (i.e. high frequency) signals, which is not the case here.
Contrary to [19] our approach is based on (approximate) maximum likelihood estimation 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 would avoid such filtering steps when necessary.
The code and datasets used to produce the numerical results of this paper are available
at the web site
https://cv.archives-ouvertes.fr/bruno-torresani

A PPENDIX A
P ROOF

OF THEOREM

1

√
To simplify notations, let Bγ denote the operator Dγ / γ′ . We split the approximation
error as
τ fτ
g
ε(s, τ ) = hAa Dγ X, ψsτ i − hA
a Dγ X, ψsτ i


q
(2)
(3)
(1)
′
= a(τ ) ε (s, τ ) + γ (τ ) ε (s, τ ) + ε (s, τ ) ,

where
(1)



p

q





γ′ − γ′ (τ ) Bγ X, ψsτ




q
p
′
′
γ − γ (τ ) ψsτ ,
= X, Bγ−1
D

E
(2)
τ
f
ε (s, τ ) =
Bγ − Bγ X, ψsτ

E
D 
−1
τ
f
ψsτ ,
= X, Bγ−1 − Bγ

ε

(s, τ ) =

21

ε

(3)


E
−1
τ
g
− 1 Dγ X, ψsτ
(s, τ ) = Aa Aa

E
D

−1
τ
g
− 1 ψsτ .
= X, Dγ−1 Aa Aa
D

Furthermore, the triangle inequality gives:
r n
o
n
E |ε(s,τ)|

2

≤

Ca2

2
ε(1) (s,τ)

E

o

r

n
o
2
+ Cγ E ε(2) (s,τ)

!
r n
o 2
2
.
+ E ε(3) (s, τ )

(30)

Let us now determine an upper bound for each error term. To this end, they are
written as follows:
E



ε

(k)

(s, τ )

2



=

Z ∞
0

2

(k)
SX (ξ ) fˆsτ (ξ ) dξ ,

with k ∈ {1, 2, 3}.

√
Concerning the first error term, a Taylor expansion of γ′ around τ gives

Z q
q
(
1
)
fˆsτ (ξ ) =
γ′ (t) − γ′ (τ ) ψsτ (t)e−2iπγ(t)ξ dt
R

where Iψ =
write

R

≤
R

Z

R

γ′′
√ ′
2 γ

3s

∞

|t − τ | |ψsτ (t)|dt ≤ q 2

kγ′′ k∞
Iψ ,
√
2 cγ

|tψ(t)|dt. Furthermore, the localization assumption on ψ allows us to
Iψ ≤ 2

Z ∞
0

t
dt ≤ 2
1 + tβ

Z

1
0

tdt +

Z ∞
1

1
t β−1

dt



=

β
.
β−2

Finally, we can control the first error term as follows:
2
 

′′
2
βσX
3s/2 kγ k∞
(1)
.
≤ q
E ε (s, τ )
√
c γ 2( β − 2)

Concerning the second error term, we have
Z

′
(2)
ˆf sτ
e−2iπγ(t)ξ − e−2iπ (γ(τ )+(t−τ )γ (τ ))ξ ψsτ (t)dt
(ξ ) =
Z

R

′

1 − e−2iπ (γ(τ )+(t−τ )γ (τ )−γ(t))ξ |ψsτ (t)| dt
R
Z
π

≤
2 sin
ξ (t − τ )2 γ′′ (t∗ ) |ψsτ (t)| dt ,
2
R

≤

for some t∗ between t and τ. Besides, we have | sin(u)| ≤ |u| and | sin(u)| ≤ 1 so that:

Z
Z
π
(
2
)
′′
2s
2
s/2
ξ kγ k∞ q t |ψ(t)| dt +
fˆsτ (ξ ) ≤ 2q
| ψ(t)| dt .
R\ J
J2

22

where J = [− T, T ]. One can prove that the value of T minimizing the the right-hand
−1/( β+2)
. Therefore, we have:
side of the latter equation is T = π2 ξ kγ′′ k∞ q2s
4( β + 2)
(2)
fˆsτ (ξ ) ≤ q 2( β+2)
3( β − 1)
5β−2

π
2

′′

ξ kγ k∞

Finally, we can control the second error term as follows:
E



ε

(2)

(s, τ )

2



≤ q

5β−2
2( β +2)

 β +2
β −1

.

β +2
4( β + 2)π ′′  β−1 (ρ)
kγ k∞
IX
3( β − 1) 2

!2

.

Concerning the third error term, we have


Z q
a
(
t
)
(
3
)
fˆsτ (ξ ) =
γ′ (t)
− 1 ψsτ (t)e−2iπγ(t)ξ dt
a
(
τ
)
R
Z
1 k a′ k
1
3s
k a′ k∞
∞
2
|t − τ | |ψsτ (t)|dt = q 2 Cγ2
Iψ ,
≤ Cγ
ca
R ca
Finally, we can control the third error term as follows:
!2
p


2
C
γ βσX
E ε(3) (s, τ )
≤ q3s/2
k a′ k∞
.
ca β − 2
To conclude the proof, the three errors terms in equation (30) are replaced by their
upper bounds to obtain the approximation error given in the theorem.



A PPENDIX B
P ROOF

OF PROPOSITION

1

(k)

Let w̃ x,sm ∈ C Nτ denote the estimation of w x,sm after k iterations of the algorithm.


(k)
(k)
1 f
q
Considering equation (18), we have w̃ x,s =
Wy sm − θ̃ , τ , thus
m

E

n

(k)
S˜X (q−sm ω0 )

o

2

(k)

θ̃1

1
=
E
Nτ kψk22



(k) 2
w̃ x,sm



.

(ξ )

To simplify notations, let us introduce some variables. We define sm = sm + logq (ξ )
and h( x ) = ϕ1 (qx ) = qx |ψ̂ (qx )|2 for x ∈ R.
By means of the covariance expression given in equation (12) we can write
 
 Nτ


 
1
(
k
)
(
k
)
(k) 2
fy sm − θ̃ , τn
fy sm − θ̃ , τn W
E w̃x,sm =∑ (k)E W
n,2
n,2
n =1 θ̃n,1
Nτ

=

∑

θn,1

(k)
n =1 θ̃n,1

Z ∞
SX (ξ )  (ξ )
0

ξ


(k)
h sm + θn,2 − θ̃n,2 dξ

23



Z ∞
SX (ξ ) Nτ θn,1  (ξ )

=

ξ

0

∑

(k)
n =1 θ̃n,1

h sm

(k)



(k)
+ θn,2 − θ̃ dξ.
n,2

Let us now split the bias into two terms such that bSX (m) = g1 (m) + g2 (m), where g1
and g2 are defined as

 





Z

N −1 ∞SX (ξ ) Nτθn,1   (ξ )
(k) 
−
1
+
θ
−
θ̃
h
s
g1 ( m ) = τ 2
n,2
m
n,2 dξ,
(k)
ξ
k ψ k2 0
θ̃
n =1 n,1

∑

N −1
g2 ( m ) = τ 2
k ψ k2

Z ∞
SX (ξ )
0

ξ

Nτ  

∑

n =1



(k)

(ξ )

h sm + θn,2 − θ̃n,2

!
 
(ξ )
− h sm dξ.

Regarding the first term, we directly have
(k) Z
∞ S (ξ )
khk∞ 1 Nτ θn,1 − θ̃n,1
X
| g1 (m)| ≤
dξ
∑
2
(
k
)
ξ
kψk2 Nτ n=1
0
θ̃n,1

≤

khk∞
(k)
θ1 − θ̃1
2
k ψ k 2 c θ1

∞

JX .

Besides, we have khk∞ = k ϕ1 k∞ and the smoothness and decay assumptions on ψ̂ allow

us to write ϕ1 (u) = Ou→∞ (u1−2η ) → 0. Then ϕ1 is bounded and K1′ < ∞. This yields
u→∞
JX K1′
(k)
| g1 (m)| ≤
.
θ1 − θ̃1
2
∞
k ψ k2
(ξ )

Regarding the second term, a Taylor expansion of h around sm gives
Z
  
N −1 ∞SX(ξ ) Nτ  (ξ )
(ξ )
(k)
+
θ
−
θ̃
| g2 (m)| ≤ τ 2
h
s
n,2
m
∑
n,2 − h s m dξ
ξ
k ψ k2 0
n =1

≤
≤
Furthermore, ∀ x ∈ R

kh′ k∞ 1
kψk22 Nτ
kh′ k

∞
kψk22

Nτ

∑
n =1

(k)

θn,2 − θ̃n,2
(k)

θ2 − θ̃2

∞

Z ∞
SX (ξ )

ξ

0

dξ

JX .

2

|h′( x )|ln (q)−1 = qxϕ1′ (qx ) ≤ qx ψ̂(qx ) + 2q2x ψ̂ (qx )ψ̂′(qx )
≤ k ϕ1 k∞ + 2kψ̂′ k∞ k ϕ2 k∞ = ln(q)−1 K2′ .

Besides, the decay assumption on ψ̂ gives | ϕ2 (u)| = Ou→∞ (u2−η ) → 0 because η > 2.

Then ϕ2 is bounded and K2′ < ∞. This yields
JX K2′
(k)
| g2 (m)| ≤
θ2 − θ̃2
2
k ψ k2

u→∞

∞

.

The proof is concluded by summing up the upper bounds of | g1 | and | g2 | to obtain the
(k)

upper bound of bSX . Notice that this bound does not depends on m.



24

R EFERENCES
[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, 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. [Online]. Available: http://dx.doi.org/10.1214/aos/1030563977
[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. [Online]. Available: http://www.jstor.org/stable/3088859
[6] R. Dahlhaus, “A likelihood approximation for locally stationary processes,” Ann. Statist., vol. 28, no. 6, pp.
1762–1794, 12 2000. [Online]. Available: http://dx.doi.org/10.1214/aos/1015957480
[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. [Online]. Available: http://dx.doi.org/10.1111/1467-9892.00200
[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. Torré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. [Online]. Available:

http://hal.archives-ouvertes.fr/hal-00822186
[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. Torré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.

25

[18] 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.
[19] ——, “Estimating deformations of stationary processes,” Ann. Statist., vol. 31, no. 6, pp. 1772–1821, 12 2003.
[Online]. Available: http://dx.doi.org/10.1214/aos/1074290327
[20] H. Omer and B. Torré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. [Online]. Available: https://hal.archives-ouvertes.fr/hal-01094835
[21] A. Meynard and B. Torrésani, “Spectral estimation for non-stationary signal classes,” in Sampling
Theory and Applications, ser. Proceedings of SampTA17, Tallinn, Estonia, Jul. 2017. [Online]. Available:
https://hal.archives-ouvertes.fr/hal-01475534
[22] H. Omer, “Modèles de déformation de processus stochastiques généralisés. application à l’estimation des non
stationnarités dans les signaux audio.” Ph.D. dissertation, Aix-Marseille Université, 2015.
[23] 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.
[24] P. Stoica and R. Moses, Introduction to Spectral Analysis. Prentice Hall, 1997.

