Adrien Meynard

Bruno Torrésani

Spectral analysis for nonstationary audio

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

A new approach for the analysis of nonstationary signals is proposed, with a focus on audio applications. Following earlier contributions, nonstationarity is modeled via stationaritybreaking 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.

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 [START_REF] Koopmans | The spectral analysis of time series[END_REF], [START_REF] Priestley | Spectral analysis and time series, ser. Probability and mathematical statistics[END_REF] and references therein) gives a clear meaning to the notion of stationarity. In the context of time series analysis, Priestley [START_REF] Priestley | Spectral analysis and time series, ser. Probability and mathematical statistics[END_REF], [START_REF]Non-linear and non-stationary time series analysis[END_REF] was A. Meynard and B. Torrésani are with Aix Marseille Univ, CNRS, Centrale Marseille, I2M, Marseille, France. 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 [START_REF] Mallat | Adaptive covariance estimation of locally stationary processes[END_REF], who proposed a wavelet-based approach to covariance estimation for locally stationary processes (see also [START_REF] Nason | Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum[END_REF]). An alternate theory of locally stationary time series was developed by Dahlhaus [START_REF] Dahlhaus | A likelihood approximation for locally stationary processes[END_REF] (see also [START_REF] Sachs | A wavelet-based test for stationarity[END_REF] for a corresponding stationarity test). In a different context, frequency-modulated stationary signal were considered in [START_REF] Wisdom | Extending coherence time for analysis of modulated random processes[END_REF], [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF], and time warping models were analyzed in [START_REF] Wisdom | Improved statistical signal processing of nonstationary random processes using time-warping[END_REF]. 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) = ∑ K k=1 A k (t) cos(2πφ k (t)). This model has been largely used in speech processing since early works by McAulay and Quatieri [START_REF] Mcaulay | Speech analysis/synthesis based on a sinusoidal representation[END_REF] (see [START_REF] Lagrange | Enhancing the tracking of partials for the sinusoidal modeling of polyphonic sounds[END_REF] and references therein for more recent developments, and [START_REF] Depalle | Tracking of partials for additive sound synthesis using hidden markov models[END_REF], [START_REF] Godsill | Bayesian harmonic models for musical pitch estimation and analysis[END_REF] for probabilistic approaches). The instantaneous frequencies φ ′ k of each mode give important information about the physical phenomenon. Under smoothness assumptions on functions A k and φ ′ k , techniques such as ridge/multiridge detection (see [START_REF] Carmona | Practical time-frequency analysis: Gabor and Wavelet Transforms With an Implementation[END_REF] and references therein), synchrosqueezing or reassignment have been developed to extract theses quantities from a single signal observation (see [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], [START_REF] Wu | Instantaneous frequency and wave shape function (i)[END_REF], [START_REF] Lin | Wave-shape function analysis[END_REF] 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 φ ′ k is multiple of a fundamental frequency φ ′ 0 : φ ′ k (t) = (k + 1)φ ′ 0 (t). In the special case A k (t) = α k A 0 (t), we can describe such signals as a stationary signal x(t) = ∑ K k=1 α k cos(2πkt + ϕ k ) modified by time warping and amplitude modulation: y(t) = A 0 (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, deformed by a stationarity-breaking operator (e.g. time warping, amplitude modulation). For example, sounds generated by a variablespeed engine or any stationary sound deformed by Doppler effect can be described as such. A stochastic time warping model has been introduced in [START_REF] Clerc | The texture gradient equation for recovering shape from texture[END_REF], [START_REF]Estimating deformations of stationary processes[END_REF], where wavelet-based approximation and estimation techniques were developed. In [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF], 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 [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF] 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 [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF]) 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 [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF] 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 ∈ R → X t ∈ C. Signals of interest are real-valued, however we will work with complexvalued 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 (E {X t = 0} for all t) and be second-order, i.e. they have a well-defined covariance kernel E X t 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 ′ ) ∆ = E X t X t ′
is a function of tt ′ 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η X , which we will assume to be continuous with respect to the Lebesgue measure: dη X (ν) = S X (ν) dν, for some nonnegative L 1 function S X called the power spectrum. We then write

C X (t) = S X (ν)e 2iπνt dν .
We refer to textbooks such as [START_REF] Koopmans | The spectral analysis of time series[END_REF], [START_REF] Priestley | Spectral analysis and time series, ser. Probability and mathematical statistics[END_REF] for a more complete mathematical account of the theory, and to [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF] 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 α , translation T τ , dilation D s , and frequency modulation M ν defined as follows:

A α x(t) = αx(t) , T τ x(t) = x(t -τ) , D s x(t) = q -s/2 x(q -s t) , M ν x(t) = e 2iπνt x(t) .
where α, τ, s, ν ∈ R and q > 1 is a fixed number.

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

T τ D s = D s T q -s τ , T τ M ν = e -2iπντ M ν T τ , M ν D s = D s M ν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 ∈ R → X t is defined as:

W X (s, τ) ∆ = X, ψ sτ , with ψ sτ = T τ D s ψ . ( 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 (see [START_REF] Carmona | Practical time-frequency analysis: Gabor and Wavelet Transforms With an Implementation[END_REF]), 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 ψ, the wavelet transform can still be well defined (see [START_REF] Carmona | Practical time-frequency analysis: Gabor and Wavelet Transforms With an Implementation[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF] 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 space

H 2 (R) = ψ ∈ L 2 (R) : supp( ψ) ⊂ 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 (R) are (analytic) derivative of Gaussian ψ k (which has k vanishing moments), and the sharp wavelet ψ ♯ (with infinitely many vanishing moments) introduced in [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF]. 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 ψ and ν 1 is a cutoff frequency. ǫ is a prescribed numerical tolerance chosen so that ψ♯ (ν 1 ) = ǫ, and

the divergence δ is defined by δ(a, b) = 1 2 a b + b a -1.
The quality factor of ψ ♯ , i.e. the center frequency to bandwidth ratio, can be expressed as Q = 1/ C(C + 4), where C = -δ(ν 1 , ν 0 ) ln 2/ ln ǫ > 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,

A a : A a x(t) = a(t)x(t) , (3)
where a ∈ C 1 is a real valued function, such that

0 < c a ≤ a(t) ≤ C a < ∞, ∀t , (4)
for some constants c a , C a ∈ R * + . Time warpings are compositions with smooth and monotonic functions,

D γ : D γ x(t) = γ ′ (t)x(γ(t)) , (5)
where γ ∈ C 2 is a strictly increasing smooth function, satisfying the control condition [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF] 0

< c γ ≤ γ ′ (t) ≤ C γ < ∞, ∀t , (6)
for some constants c γ , C γ ∈ 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 [START_REF] Wisdom | Improved statistical signal processing of nonstationary random processes using time-warping[END_REF], 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 [START_REF] Stowell | Computational Analysis of Sound Scenes and Events[END_REF] 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 [START_REF] Senin | Dynamic time warping algorithm review[END_REF] 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 [START_REF] Zhang | Harmonic transform[END_REF] or the Fan-Chirp transform [START_REF] Weruaga | The fan-chirp transform for non-stationary harmonic signals[END_REF] (see [START_REF] Dunn | Sinewave analysis/synthesis based on the fan-chirp transform[END_REF] 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 = A a D γ X ( 7
)
where X is a stationary zero-mean real random process with (unknown) power spectrum S X . The goal is to estimate the deformation functions a and γ 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 γ(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 [START_REF]Estimating deformations of stationary processes[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF]. 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 γ, and their slow variations. This allows us to perform a local analysis using smooth and localized test functions, on which the action of A a and D γ can be approximated by their so-called tangent operators A τ a and D τ γ (see [START_REF]Estimating deformations of stationary processes[END_REF], [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF], [START_REF] Omer | Modèles de déformation de processus stochastiques généralisés. application à l'estimation des non stationnarités dans les signaux audio[END_REF]). Given a test function g located near t = τ (i.e. decaying fast enough as a function of |t -τ|), Taylor expansions near

t = τ yield A a g(t) ≈ A τ a g(t) , with A τ a ∆ = A a(τ) , (8)
D γ g(t) ≈ D τ γ g(t) , with D τ γ ∆ = T τ D -log q (γ ′ (τ)) T -γ(τ) . (9)
Therefore, we approximate the wavelet transform of

Y by W Y (s, τ) ≈ W Y (s, τ) ∆ = A τ a D τ γ X, T τ D s ψ , i.e. W Y (s, τ) = a(τ)W X s + log q (γ ′ (τ)), γ(τ) . ( 10
)
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. Hereafter, we denote by f ∞ = 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 [START_REF] Sachs | A wavelet-based test for stationarity[END_REF]. Let ψ be a smooth test function, localized in such

a way that |ψ(t)| ≤ 1/(1+|t| β ) for some β > 2. Let W Y be the wavelet transform of Y,
W Y its approximation given in [START_REF] Wisdom | Improved statistical signal processing of nonstationary random processes using time-warping[END_REF], and let ε = W Y -W Y denote the approximation error. Assume ψ and S X are such that

I (ρ) X ∆ = ∞ 0 ξ 2ρ S X (ξ) dξ < ∞ , where ρ = β -1 β + 2 .
Then the approximation error ε is a second order, two-dimensional complex random field, and

E |ε(s,τ)| 2 ≤ C 2 a q 3s K 1 γ ′′ ∞ +K 2 q µs γ ′′ ρ ∞ +K 3 a ′ ∞ 2
where

K 1 = βσ X 2(β -2) √ c γ , K 2 = I (ρ) X π 2 ρ C γ 4 3ρ , K 3 = C γ βσ X (β -2)c a , µ = β -4 β + 2 , σ 2 X being the variance of X.
The proof, which is an extension of the one given in [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF], rests on Taylor approximations of γ ′ and a in the neighborhood of t = τ and subsequent integral majorizations, is given as supplementary material.

Remark 3:

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 [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF].

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, γ from W Y , which is a zero-mean random field with covariance

E W Y (s, τ) W Y (s ′ , τ ′ ) = C(s, s ′ , τ, τ ′ ) ( 11
)
where the kernel C reads

C(s, s ′ , τ,τ ′ ) = a(τ)a(τ ′ )q s+s ′ 2 γ ′ (τ)γ ′ (τ ′ ) ∞ 0 S X (ξ) × ψ q s γ ′ (τ)ξ ψ q s ′ γ ′ (τ ′ )ξ e 2iπξ(γ(τ)-γ(τ ′ )) dξ. ( 12
)
B. Estimation 1) Estimation procedure: Our goal is to estimate both deformation functions γ and a from the approximated wavelet transform 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, W Y is a zero-mean circular Gaussian random field and its probability density function is characterized by the covariance matrix. However, equation [START_REF] Lagrange | Enhancing the tracking of partials for the sinusoidal modeling of polyphonic sounds[END_REF] shows that besides deformation functions the covariance also depends on the power spectrum S X 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 S X . This constraint naturally brings the estimation strategy to an alternate algorithm. In [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF], an estimate for S X was obtained at each iteration using a Welch periodogram on a "stationarized" signal

A -1 ã D -1 γ Y, ã
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 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 S X is known (in fact, only an estimate SX 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 F s the time sampling frequency. The sizes of s and τ are respectively denoted by M s and N τ .

Considering the covariance expression [START_REF] Lagrange | Enhancing the tracking of partials for the sinusoidal modeling of polyphonic sounds[END_REF] we want to estimate the vector of pa-

rameters Θ = (θ 1 , θ 2 , θ 3 ) ∆ = (a(τ) 2 , log q (γ ′ (τ)) , γ(τ)). Let W y = W y (Λ) denote the
discretized transform and let C W (Θ) be the corresponding covariance matrix. The related log-likelihood is

L (Θ) = - 1 2 ln |det(C W (Θ))| - 1 2 C W (Θ) -1 W y • W y . ( 13
)
The matrix C W (Θ) is a matrix of size M s N τ × M s N τ , 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 (Θ) has about 3.1 trillion elements which makes it numerically intractable. In addition, due to the redundancy of the wavelet transform, C W (Θ) 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,τ n = W y (s, τ n ) as independent circular Gaussian vectors with zero-mean and covariance matrix

C(Θ n ) ij = θ n,1 C 0 (θ n,2 ) ij , 1 ≤ i, j ≤ M s , (14)
where

C 0 (θ n,2 ) ij = q (s i +s j )/2 ∞ 0 S X (q -θ n,2 ξ) ψ(q s i ξ) ψ(q s j ξ) dξ . ( 15
)
In this situation, the regularized log likelihood L r splits into a sum of independent terms

L r (Θ) = ∑ n L (Θ 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 anymore in the covariance expression. Thus, we are led to maximize independently for each n

L (Θ n ) = - 1 2 ln |det(C(Θ n ))| - 1 2 C(Θ n ) -1 w y,τ n • w y,τ n . ( 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 θ n,2 using a gradient ascent method, for a fixed value of θ n,1 .

Second, for a fixed θ n,2 , an estimate for θ n,1 is directly obtained which reads

θn,1 = 1 M s C -1 0 (θ n,2 )w y,τ n • w y,τ n . ( 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

θ 1/2 n,1 W y (s -θ n,2 , τ n ) = W x (s, γ (τ n )) , (18)
For fixed scale s m , w x,s m ∆ = W x (s m , γ(τ)) ∈ 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 the so-called wavelet spectrum

S X,ψ (q -s m ω 0 ) ∆ = E 1 N τ ψ 2 2 w x,s m 2 (19) = 1 ψ 2 2 ∞ 0 S X (ξ)q s m ψ (q s m ξ) 2 dξ , ( 20
)
where ω 0 is the central frequency of | ψ| 2 . S X,ψ is a narrowband version of S X centered around frequency ν m = q -s m ω 0 . Besides, the bandwidth of the filter is proportional to the frequency ν m . This motivates the introduction of the following estimator SX of S X SX (q -s m ω 0 )

∆ = 1 N τ ψ 2 2 w x,s m 2 . ( 21
)
Finally, the estimate SX is extended to all ξ ∈ [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 [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF] applied to the observation Y.

After k iterations of the algorithm, estimates Θ(k) n and S (k) X for Θ n and S X are available. Hence we can only evaluate the plug-in estimate C(k) 0 of C 0 , obtained by replacing the power spectrum with its estimate in the covariance matrix [START_REF] Carmona | Practical time-frequency analysis: Gabor and Wavelet Transforms With an Implementation[END_REF]. 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.

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

θ(k) j - θ(k-1) j 2 2 θ(k-1) j 2 2 < Λ , for j = 1, 2 , (22)
where 0 < Λ < 1 is a user defined threshold.

Finally, after convergence of the algorithm to the estimated value Θ(k) , log q (γ ′ ) and a 2 are estimated through time by cubic spline interpolation. Besides, γ is given by numerical integration assuming that γ(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 θ1 or θ2 . Indeed, the computation of the log-likelihood involves the evaluation of the inverse covariance matrix. In [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF],

a sufficient condition for invertibility was given in the presence of 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 M s on its Cramér-Rao bound (see section III-C below).

The choice we made is to maximize L (Θ n ) for θ n,2 with w y,τ 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 ⌊M s /p⌋ for some p ∈ N * . by maximizing the approximate log-likelihood

L (k) θ n,1 , θ(k+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 ( 23).

• Construct the estimated wavelet transform W x of the underlying stationary signal by interpolation from W y and θ(k) with equation [START_REF] Lin | Wave-shape function analysis[END_REF]. Estimate the corresponding power spectrum S (k+1)

X with (21).
• k := k + 1 end while

• Compute ã and γ by interpolation from Θ(k) .

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 invertibility. The regularized matrix is constructed by replacing covariance matrix C 0 (θ n,2 ) in ( 15) by its regularized version C 0,r (θ n,2 ), given by

C 0,r (θ n,2 ) = (1 -r)C 0,r (θ n,2 ) + rI , ( 23
)
for some regularization parameter 0 ≤ r ≤ 1.

Remark 6: After convergence of the estimation algorithm, the estimated functions ã and γ allow constructing a "stationarized" signal

x = D γ-1 A ã-1 y .
x is an estimation of the original underlying stationary signal x. Furthermore, the Welch periodogram [START_REF] Welch | The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF] may be computed from x to obtain an estimator of 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 γ 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 θ n,1 , the estimator is unbiased when the actual values of θ n,2 and S X are known. In our case, the bias b

(k) n,1 (θ n,1 ) = E θ(k) n,1 -θ n,1 is written as b (k) n,1 (θ n,1 ) = θ n,1 M s Trace C(k) 0 θ(k) n,2 -1 C 0 (θ n,2 )-I . ( 24
)
As expected, the better the covariance matrix estimation, the lower the bias b

(k)
n,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 M s → ∞).

Regarding S X , equation [START_REF] Clerc | The texture gradient equation for recovering shape from texture[END_REF] 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 ψ ∈ H 2 (R) be an analytic wavelet such that ψ is bounded and ψ(u) = O u→∞ (u -η ) with η > 2. Let ϕ 1 and ϕ 2 be bounded functions defined on R + by ϕ 1 (u) = u ψ(u) 2 and ϕ 2 (u) = u 2 ψ(u) . Assume S X is such that

J X = ∞ 0 ξ -1 S X (ξ) dξ < ∞.

Let S (k)

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

(k) S X denote the bias defined for all m ∈ [[1, M s ]] by b (k) S X (m) = E S (k) X (q -s m ω 0 ) -S X,ψ (q -s m ω 0 ) .
Assume there exists a constant c θ 1 > 0 such that θ

(k) n,1 > c θ 1 , ∀n, k. Then b (k) S X ∞ ≤ J X ψ 2 2 K ′ 1 θ 1 - θ(k) 1 ∞ + K ′ 2 θ(k) 2 -θ 2 ∞ , ( 25
)
where

K ′ 1 = ϕ 1 ∞ c θ 1 < ∞ , K ′ 2 = ln(q) ϕ 1 ∞ + 2 ψ′ ∞ ϕ 2 ∞ < ∞ .
The proof of the Proposition is given in supplementary materials.

Remark 8: If θ (k) 1 → θ 1 and θ (k) 2 → θ 2 as k → ∞, we have E S (k) X (ν m ) → k S X,ψ (ν m ),
as expected.

Formula ( 25) enables the control of the spectrum bias at frequencies ν m = q -s m ω 0 only. Notice also that the requirement J X < ∞ forces 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 [START_REF] Stoica | Introduction to Spectral Analysis[END_REF]) 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 S X and the other parameters are well estimated.

Applying this formula to θ n,1 gives

E θn,1 -E θn,1 2 ≥ CRLB(θ n,1 ) = 2θ 2 n,1 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 θ 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 observations are corrupted by a random Gaussian white noise W with variance σ 2 W (supposed to be known):

Y = A a D γ X + W . ( 26
)
The estimator θn,1 is not robust to noise. Indeed, if the maximum likelihood estimator of model [START_REF] Sachs | A wavelet-based test for stationarity[END_REF] in the presence of such white noise, a new term b

(k) n,1| W (θ n,1
) must be added to the bias expression [START_REF] Senin | Dynamic time warping algorithm review[END_REF], which becomes

b (k) n,1| W (θ n,1 ) = 1 M s Trace C(k) 0 θ(k) n,2 -1 C wn ,
where (C wn ) ij = σ 2 W q (s i +s j )/2 ∞ 0 ψ(q s i ξ) ψ(q s j ξ)dξ. 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

C(Θ n ) ij = q s i +s j 2 ∞ 0 (θ n,2 S X (q -θ n,1 ξ) + σ 2 W ) ψ(q s i ξ) ψ(q s j ξ) dξ (27)
and the likelihood is modified accordingly. Formula ( 17) is no longer true and no closedform expression can be derived anymore, 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 [START_REF] Dunn | Sinewave analysis/synthesis based on the fan-chirp transform[END_REF] shows that the only change in the covariance matrix formula is to replace the power spectrum S X by S Z =

S X + σ 2 W θ n,2 .
The additive constant term does not impair the estimator as long as it is small in comparison with the maximum values of S X . Moreover, the estimator SX is modified because when computing 1

θ 1/2 n,1 W y (s -θ n,2 , τ n )
on scale s m , we compute:

w z,s m = w x,s m + w w * ,s m , where w w * ,s m = 1 θ 1/2 1 W w (s m -θ 2 , τ)
is the wavelet transform of a white noise modulated in amplitude by a -1 . Thus a constant term σW independent of frequency is added to the new spectrum estimator SZ , so that [START_REF] Koopmans | The spectral analysis of time series[END_REF] .

E SZ = S X,ψ + σ2 W where σ2 W = σ 2 W 1 N τ N τ ∑ n=1 1 θ n,

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 [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF], we shortly account for this case here for the sake of completeness. Let α ∈ C 2 be a smooth function, and set

M α : M α x(t) = e 2iπα(t) x(t) , (28)
The deformation model in [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF] is of the form

Y = A a M α D γ X . (29)
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, ν, τ) = X, ψ sντ , with ψ sντ = T τ M ν D s ψ.
In that case, approximation results similar to Theorem 1 can be obtained from which the corresponding log-likelihood can be written. At fixed time τ, the estimation strategy 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 and frequency variables s and ν 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 ψ ♯ 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 [START_REF]Estimating deformations of stationary processes[END_REF] was also implemented, but we couldn't get satisfactory results with that approach.

• Amplitude modulation: we use as baseline estimator of a(τ n ) 2 the average energy θ(B) n,1 defined as follows:

θ(B) n,1 = 1 M s w y,τ n 2 .
This amounts to replace the estimated covariance matrix in [START_REF] Wu | Instantaneous frequency and wave shape function (i)[END_REF] by the identity matrix. Notice that θ(B) n,1 does not depend on the time warping estimator, and can be computed directly on the observation.

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

θ(B) n,2 = C 0 + 1 w y,τ n 2 M s ∑ m=1 s[m]|w y,τ n [m]| 2 . C 0 is chosen such that θ(B)
2 is a zero-mean vector.

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 ≈ 8.2 s). The spectrum S X is written as S X = S 1 + S 2 where S l (ν) = 1 + cos 2π(νν

(l) 0 )/∆ (l) ν if |ν -ν (l) 0 | < ∆ (l)
ν /2 and vanishes elsewhere (for l ∈ {1, 2}). The amplitude modulation a is a sine wave a(t) = a 0 (1 + a 1 cos(2πt/T 1 )), where a 0 is chosen such that t -1 F t F 0 a 2 (t)dt = 1. The time warping function γ is such that log q (γ ′ (t)) = Γ + cos(2πt/T 2 )e -t/T 3 , where Γ is chosen such that t -1

F t F 0 γ ′ (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(ǫ) = -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). The left hand side of Fig. 2 displays the estimated spectrum given by formula [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF].

The agreement with the actual spectrum is very good, with a slight enlargement effect due to filtering by | ψ| 2 . The right hand side of Fig. 2 gives the evolution of the stopping criterion [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF] 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 [START_REF] Stowell | Computational Analysis of Sound Scenes and Events[END_REF]. The wavelet transform of this signal in Fig. 3 shows that the warping model [START_REF] Sachs | A wavelet-based test for stationarity[END_REF]

ν = 200 Hz, ν (1)
0 = 1.2 kHz, ∆ (2)
ν = 400 Hz): actual (dash-dot blue line) and estimated (solid red line) spectra. Right: Relative update evolution. 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,

γ ′ (t) = c 2 c 2 -V 2 1 - V 2 t d 2 (c 2 -V 2 ) + (cVt) 2 . ( 30
)
We plot in Fig. 4 (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. 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 [START_REF] Omer | Estimation of frequency modulations on wideband signals; applications to audio signal analysis[END_REF], [START_REF] Omer | Modèles de déformation de processus stochastiques généralisés. application à l'estimation des non stationnarités dans les signaux audio[END_REF], [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF], [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF] for the joint estimation of deformation operator and power spectrum for deformed stationary signals, a problem already addressed in [START_REF]Estimating deformations of stationary processes[END_REF] 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 ( [START_REF] Meynard | Spectral estimation for non-stationary signal classes[END_REF] was limited to time warping and combinations of time warping with frequency modulation, and investigated generalized wavelet transforms); 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 [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling[END_REF] 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 [START_REF]Estimating deformations of stationary processes[END_REF], its domain of validity seems to be limited to small-scale (i.e. high-frequency) signals, which is not the case here.

Contrary to [START_REF]Estimating deformations of stationary processes[END_REF], 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 [START_REF] Podos | A fine-scale, broadly applicable index of vocal performance: frequency excursion[END_REF] 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).

Fig. 1 .

 1 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).

Fig. 2 .

 2 Fig. 2. Left: Spectrum estimation (ν (1) 0 = 600 Hz, ∆

Fig. 3 .

 3 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.

Fig. 4 .

 4 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).

  Algorithm 1 JEFAS (Joint Estimation of Frequency, Amplitude and Spectrum)

	Initialization: Compute an estimate SY of the power spectrum of Y as an initial
	guess	S (0) X for S X . Initialize the estimator of the squared amplitude modulation with
	θ(0) n,1 = 1, ∀n.
	Compute the wavelet transform W y of y.
	k := 1
	while criterion (22) is false and k ≤ k max do • For each n, subsample w y,τ n on scales s p , and estimate approximate log-likelihood L (k) θ(k) n,1 , θ n,2 in (16).	θ(k+1) n,2	by maximizing the
	• For each n, estimate	θ(k+1) n,1

  fits well this kind of signal, except for transient clicks that are not

		2								10 0
		1.8
		1.6								10	-1
	Spectrum	0.6 0.8 1 1.2 1.4							Relative update	10 10 -2 -3	T = 0.1 %
		0.4										Time warping
		0.2								10 -4	Amplitude modulation Spectrum
		0	400	600	800	1000	1200	1400	1600		1	2	3	4	5	6	7
					Frequency (Hz)							Iteration

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.


Alexis Fouineau

Marie-Ange Raulet

Bruno Lefebvre

Noël Burais

Fabien Sixdenier

Semi-analytical methods for calculation of leakage inductance and frequency-dependent resistance of windings in transformers

Keywords: Eddy currents, Electric resistance, Inductance, Power transformers

Leakage inductance and frequency-dependent resistance of windings are key parameters for the design of magnetic devices. This report focuses on analytical or semi-analytical calculation techniques to be integrated into a medium frequency transformer optimization process. To this extend, a review of currently most used method to calculate inductance and resistance of windings transformer is established. Based on this review, new calculation methods are developed. Magnetic field cartography is evaluated on the basis of two-dimensional structures thanks to a combination of Ampère's law and image method. From this magnetic field map, leakage inductance can be calculated. The obtained accuracy with this method is equivalent to three-dimensional finite element modeling. For frequency-dependent resistance evaluation, new methods based on two-dimensional or one-dimensional approach are compared to finite element modeling and reviewed models. An analysis of behavior of each model is done, showing that newly developed models are more suitable in the case of medium frequency transformer design in terms of accuracy, robustness and calculation time. This work could be useful to magnetic devices designers because newly developed models could easily be adapted to various geometries, such as toroidal structures or even inductors.

I. INTRODUCTION

In the field of energy conversion, power electronics will play a major role. For example, the use of solid-state transformers can provide new solutions to future electrical grids and traction networks [START_REF] Kolar | Solid-state-transformers: key components of future traction and smart grid systems[END_REF]. It allows better connection of new renewable energy source such as far offshore wind farms [START_REF] Prasai | A new architecture for offshore wind farms[END_REF]- [START_REF] Lagier | A comparison of insulated DC-DC converters for HVDC off-shore wind farms[END_REF], or reduction of energy consumption for traction systems thanks to higher efficiency and lower weight [START_REF] Dujic | Power electronic traction transformer technology[END_REF].

A key component of solid-state transformers is the medium frequency transformer, typically operating between 500 Hz and 20 kHz to reduce its size. Its main role is to provide a galvanic insulation but the converter also takes advantage of other intrinsic properties of transformers: transformation ratio, magnetizing and leakage inductances. The design of medium frequency transformer must take into account all these specifications, which involves a multi-objective optimization process [START_REF] Garcia-Bediaga | Multiobjective optimization of medium-frequency transformers for usolated soft-switching converters using a genetic algorithm[END_REF]- [START_REF] Leibl | Design and experimental analysis of a medium-frequency transformer for solid-state transformer applications[END_REF].

In this paper, we will focus on a crucial active part of the medium frequency transformers which is the windings. Two windings characteristics must be evaluated accurately during the design process: AC resistance and leakage inductance [START_REF] Mogorovic | Medium frequency transformer leakage inductance modeling and experimental verification[END_REF].

The design process usually need low calculation time models to allow a great number of calculation points in the shortest time possible. Therefore, most of the accurate numerical solutions, such as finite element modeling (FEM), leading to high calculation time are not suitable. Analytical or semi-analytical models are usually better solutions.

The purpose of this paper is to propose new models to calculate leakage inductance (considering steady current) and frequency-dependent resistance of windings and to compare them with existing ones. To this aim, part II will present existing and newly developed calculation methods for both inductance and resistance. Then, part III will evaluate all previously identified models on typical test cases. The accuracy and calculation time results given by the various models will be presented, and compared to numerical FEM calculation results in part III. Finally, part IV will analyze the performances of the different calculation methods and will list advantages, disadvantages and limitations of each one.

II. METHODS

In the case of medium frequency applications, it is important to note that the same magnetic fields are responsible for both stored energy and eddy currents, i.e. leakage inductance and frequency-dependent resistance. It is the reason why the calculation of magnetic field is a prerequisite.

A. Magnetic fields and leakage inductance 1) One-dimensional field of concentric windings

A classical geometry of transformer is the case of concentric primary and secondary windings inside a winding window as shown in Fig. 1. In this case, the magnetic field H is oriented along the winding axis y. The maximal value of H is located between the windings, and can be calculated thanks to Ampère's law. The expression is shown in [START_REF] Kolar | Solid-state-transformers: key components of future traction and smart grid systems[END_REF], where I is the total magnetomotive force and h is the windings' height. Finally, H increases and decreases linearly in both windings. This calculation method is widely used in transformer design because of its simplicity, but its limitations should not be forgotten. This model does not work if there is not a complete Ampère-Turns compensation between primary and secondary, such as in the case of an inductor. If the windings are not confined into a winding window or even if the winding height is inferior to the winding window height, then onedimensional field hypothesis becomes false. This is why a two-dimensional approach might become necessary.

2) Field of basic conductors

This paragraph explains how to calculate magnetic field produced by a single conductor carrying a current inside vacuum. Two basic shapes of conductor will be considered: round and rectangular. These basic conductors will be used extensively for magnetic field two-dimensional calculation in transformers.

a)

Round conductor

A round conductor of radius rc carrying a current I is considered. Magnetic field in this case is purely orthoradial and its expression inside and outside the conductor can be easily obtained from Ampère's law (2).

) 2 /( ) 2 /( ) ( 2 r I H r rI H c       if if c c r r r r   (2)
Fig. 2 shows a schematic of the configuration as well as an example of magnetic field cartography in this case.

b) Rectangular conductor

In the case of a rectangular conductor, H is not orthoradial, particularly inside and near the conductor. However, its component can be calculated by using (3) (development may be found in [START_REF] Binns | Analysis and Computation of Electric and Magnetic Field Problems[END_REF]).

                                                                                              
ab I H y x               (3)
Fig. 3 shows a schematic of the configuration with all the parameters mentioned in (3), as well as an example of magnetic field cartography in this case.

3) Conductors next to a column

In this configuration, one or several basic conductors are considered placed next to a column constituted of a high permeability material, the conductors themselves being in vacuum. This configuration might be representative of the case where the concentric windings are outside the winding window.

Using magnetic image method [START_REF] Hammond | Electric and magnetic images[END_REF], this situation is equivalent to having mirrored conductors as shown in Fig. 4. The schematic is just here to illustrate how the magnetic field can be calculated for any shape of conductors next to a magnetic circuit. Usually, transformer core material has a very high permeability with µr > 1000. But if not, equation ( 4) can be used to set the value of currents in image conductors.

I I r r image 1 1      (4)
The problem is then equivalent to multiple conductors in vacuum, and superposition principle can be used to determine total magnetic field as shown in [START_REF] Garcia-Bediaga | Multiobjective optimization of medium-frequency transformers for usolated soft-switching converters using a genetic algorithm[END_REF].

    c c conductor c c conductor y x H y x H H , ,    (5)
A similar approach has been developed in [START_REF] Duppalli | Computationally efficient leakage inductance calculation for a high-frequency core-type transformer[END_REF], by using Biot-Savart law and discretization of the windings in multiple source points. Both methods are equivalent, but the one developed here requires less discretization and should be easier to implement. Nevertheless, windings sections must be a combination of round and rectangular shapes for the method presented here (5) to be applicable.

4) Conductors inside a winding window

One-dimensional field assumption is not always valid inside a winding window, for example when windings heights are not the same or do not match the window height, which is often the case for power transformers. This is why a twodimensional magnetic field calculation might be necessary.

Applying the image method to conductors inside a winding window would theoretically create an infinity of image conductors. This is because the mirror effect of image method not only duplicates the current sources, but also the material interfaces [START_REF] Billig | The calculation of the magnetic field of rectangular conductors in a closed slot, and its application to the reactance of transformer windings[END_REF]. It creates mirrored version of magnetic mirrors. For a winding window, there are two pairs of magnetic mirrors which results in a sort of infinite tiling of the two-dimensional space as shown in Fig. 5. Mirrored cells may be indexed with m and n depending on their position relative to the original cell. A and B are respectively the winding window width and height. Nevertheless, superposition principle remains valid even for an infinity of conductors and total magnetic field may be calculated with [START_REF] Leibl | Design and experimental analysis of a medium-frequency transformer for solid-state transformer applications[END_REF], where the position of image conductor according to cell index is obtained from (7) considering geometric parameters defined in Fig. 5.

               m n conductor n b m a H H , (6)
              odd is if , even is if , ) ( odd is if , even is if , ) ( n b B nB n b nB n b m a A mA m a mA m a (7)
However, it is not possible numerically to compute the magnetic field of the theoretical infinity of cables. Therefore, in practice, a limited number of cells must be considered. From (2) and (3), it is clear that an image conductor far away from the original cell will have a low contribution to magnetic field value inside the original cell. In practice, setting an empirical limit such that only reflected conductors distant less than the maximal dimension of A and B from the original cell must be taken into account gives good results.

A similar approach has been established in [START_REF] Billig | The calculation of the magnetic field of rectangular conductors in a closed slot, and its application to the reactance of transformer windings[END_REF]. However, a preferential direction was supposed to take into account the infinity of conductors. This last assumption was too restrictive and consequently the calculated field was highly inaccurate in certain zones of the window area.

5) Core-type and shell-type leakage inductance

The previous models for calculation of two-dimensional magnetic field can now be applied to some transformer structures.

In particular, core-type and shell-type structures, which are the most common ones for power transformers, can be studied with a multiple 2D planes analysis. As shown in Fig. 6, each plane represents a portion of the geometry of the transformer and has its associated depth. The planes numbered "1" correspond to conductors inside the winding window, and the planes numbered "2" correspond to conductors next to a column but outside the winding window. Determining the total leakage inductance of a transformer can be carried out by setting perfectly opposed magnetomotive forces into the primary and secondary windings and calculating the total magnetic energy Emag. This energy will then correspond only to leakage flux because magnetizing flux will be equal to zero. To calculate this energy, firstly magnetic field is calculated for each 2D plane. Then this field is integrated over the whole plane to get the value of the plane energy E2D. This integration must be performed over a sufficiently large area for infinite or semi-infinite configurations (winding next to a column). Any numerical integration method is possible as long as a good precision and a low calculation time are guaranteed. Finally, the energies E2D are multiplied by their respective depths d2D and summed using ( 8) to obtain the total magnetic energy Emag.

   D D mag d E E 2 2 (8)
Leakage inductance can then be obtained from ( 9), where I is either primary or secondary current depending on which side the leakage inductance is considered.

2 2 1 LI E mag  (9)

B. Frequency-dependent resistance

At low frequency or direct current, winding resistance can be simply calculated from material electrical conductivity and geometric dimensions. However, when the frequency increases, magnetic field will create eddy currents in opposition to excitation current due to Faraday's law of induction.

These eddy currents can be categorized into two types representing skin effect and proximity effect [START_REF] Guillod | Litz wire losses: effects of twisting imperfections[END_REF]. Skin effect corresponds to eddy currents inside a conductor created by the current flowing in the same conductor. Proximity effect represents eddy currents inside a conductor created by currents flowing in nearby conductors. Because both phenomena originate from two independent sources, they are independent [START_REF] Acero | Frequency-dependent resistance in Litz-wire planar windings for domestic induction heating appliances[END_REF]. Total losses can be split as shown in [START_REF] Duppalli | Computationally efficient leakage inductance calculation for a high-frequency core-type transformer[END_REF], where I is the peak current inside the conductor, Hext is the magnetic field produced by nearby conductors and f is the frequency.

      f H P f I P I P P ext prox skin DC Joule , ,    (10)
A practical way to represent the increase of losses with the frequency is to use factors quantifying this increase compared to DC losses. In this paper, skin factor defined in [START_REF] Billig | The calculation of the magnetic field of rectangular conductors in a closed slot, and its application to the reactance of transformer windings[END_REF] and proximity factor defined in (12) will be used. It should be noted that skin effect losses Pskin are directly related to the current I flowing inside the considered conductor and this is why it makes sense to define skin factor in [START_REF] Billig | The calculation of the magnetic field of rectangular conductors in a closed slot, and its application to the reactance of transformer windings[END_REF]. However, proximity effect losses Pprox are indirectly related to the current I. In fact, proximity losses are due to the external magnetic field Hext. But in transformer windings, external magnetic field Hext is a function of the current I. This function depends on the transformer geometry (number of turns, number of independent conductors per turns and location of each independent conductor). So it also makes sense to define a proximity factor in [START_REF] Guillod | Litz wire losses: effects of twisting imperfections[END_REF].

Finally, the frequency dependent resistance RAC and the global resistance factor FR are defined thanks to [START_REF] Acero | Frequency-dependent resistance in Litz-wire planar windings for domestic induction heating appliances[END_REF].

  DC P S DC R AC R F F R F R    (13)
These factors will be used to consistently define the different models. Next paragraphs will present various existing models and the development of new ones.

1) Dowell's model

Dowell's model [START_REF] Dowell | Effects of eddy currents in transformer windings[END_REF] is probably the most commonly used one to calculate frequency-dependent resistance in transformer design phase. Main hypothesis of this model is the onedimensional magnetic field, as presented previously in II.A.1. Windings are supposed to be composed of m successive layers of thickness d parallel to the magnetic field. By definition, these layers must have a rectangular shape, which is very suitable for windings constituted of foil or flat conductors.

However, in the case of round wires cables or Litz cables, this model can be adapted by the introduction of a porosity factor ηw (ratio of winding height hw and window height h). Round conductors are transformed to square ones and empty space between them is taken into account through the porosity factor. In the following equations, δ represents the skin depth at the considered frequency.

     πf δ h h d d X w w 1 '    (14)         ' 2 cos ' 2 cosh ' 2 sin ' 2 sinh ' X X X X X F S    (15)           ' cos ' cosh ' sin ' sinh ' 1 3 2 2 X X X X X m F P     (16)

2) Ferreira's model

This model also considers the case of a one-dimensional field, but unlike Dowell's model, basic conductors are supposed round [START_REF] Ferreira | Improved analytical modeling of conductive losses in magnetic components[END_REF]. In the following equations, ds is the diameter of a basic conductor (strand or wire), and bern and bein are the Kelvin-Bessel functions of order n.

2   s d  (17)                    2 2 ' ' ' ' 2 bei ber ber bei bei ber F S    (18)

  P P

P P f m f K F 2 1 2     (19)                    2 2 2 2 ' ' bei ber ber bei bei ber f P     (20)
Proximity factor has been factorized into two terms in [START_REF] Sullivan | Optimal choice for number of strands in a Litz-wire transformer winding[END_REF]. KP represents the winding structure factor and is independent of the frequency, whereas fP represents the frequencydependent factor and depends only on frequency and conductor diameter.

3) Reatti and Kazimierczuk model

In [START_REF] Reatti | Comparison of various methods for calculating the AC resistance of inductors[END_REF], the authors have developed a model very similar to Ferreira's one. However, the porosity factor has been taken into account. All formulas are the same as Ferreira's model except for the winding structure factor defined by [START_REF] Prieto | Study of toroidal transformers by means of 2D approaches[END_REF].

             1 3 1 4 2 2 m K w P  (21)
Other models based on Ferreira's model exist, considering various structures for windings. For example, Litz cable structure is considered in [START_REF] Tourkhani | Accurate analytical model of winding losses in round Litz wire windings[END_REF] and [START_REF] Bartoli | Modeling Litz-wire winding losses in high-frequency power inductors[END_REF]. However, main hypotheses (one-dimensional field and layers inside winding) are the same therefore these models will not be considered in this paper.

4) Asymptotic model

Contrary to previously presented models, this one does not take into account skin effect and only represents impact of proximity losses due to external magnetic field. Therefore FS = 1 and FR = 1 + FP.

The basic formula is presented in [START_REF] Sullivan | Optimal choice for number of strands in a Litz-wire transformer winding[END_REF], and expresses the proximity losses inside a round conductor of diameter ds inside a homogeneous magnetic field H as shown in [START_REF] Wojda | Winding resistance of Litz-wire and multi-strand inductors[END_REF]. In this equation, ρ is the resistivity, µ the permeability, l the length and ω the angular frequency.

  

128 2 2 4 2 H ld P s prox  (22)
For a winding of width w, height h containing N turns placed inside a one-dimensional field such as the one described in II.A.1, a proximity factor can be expressed as shown in [START_REF] Wojda | Winding resistance and power loss of inductors with Litz and solid-round wires[END_REF]. Development of this equation may be found in Appendix A, and η corresponds to the proportion of conductive material inside the winding being considered as composed of evenly distributed round conductors with insulating material between them.

h Nw f d F s P 2 2 4 2 3 48    (23)
For a given geometry, the proximity factor is proportional to f 2 . Thus it creates a line in logarithmic scale and this is why this model is called "asymptotic".

5) Albach's model

In [START_REF] Albach | Two-dimensional calculation of winding losses in transformers[END_REF], the author has developed formulas for both skin and proximity losses inside a round conductor subject to any external magnetic field distribution.

Skin factor is expressed as shown in (25) where α is defined in [START_REF] Phung | High-frequency proximity losses determination for rectangular cross-section conductors[END_REF]. In is the modified Bessel function of first kind and a is the radius of the conductor.

  i   1 (24)              a I a I a F S    1 0 Re 2 1 ( 25
)
Proximity losses are calculated thanks to spatial Fourier series decomposition. In (26), ak and bk (resp. ck and dk) are the decomposition coefficients of Hr (resp. Hφ) on the surface of the conductor in polar coordinates.

                        a I a I a d a c b l P k k k k k k k prox      1 1 2 2
Re 2 (26) However, windings composed of Litz wires have strands whose size is usually very small compared to the size of the winding itself. Therefore external magnetic field can be considered constant over a strand. In this case, previous equation can be simplified into (27) as shown in Appendix B.

              a I a I a H l P prox      0 1 2 Re 2 (27)
Finally, a proximity factor can be expressed under the assumption of a one-dimensional field with the same method used for asymptotic model, as shown in (28). The development of this equation may be found in Appendix A.

                a I a I a h Nw f K F P P P     0 1 Re 3 4 ( 28
)
This proximity factor was split into a frequencyindependent factor KP and a frequency-dependent factor fP, as done previously with Ferreira's model. In fact, Albach and Ferreira models have the same skin factors FS and frequencydependent factor fP (demonstration in Appendix C); the only difference between them lies in the expression of the winding structure factor KP.

III. RESULTS

A. Magnetic fields and leakage inductance

To test out the models of magnetic field and leakage inductance calculation, structures presented on Fig. 6 were used. The considered geometry has a winding window of 100 mm  200 mm. Windings are 150 mm high, the distance between the first winding and the core is 10 mm, and the distance between the second winding the and first winding is 20 mm. Both windings are comprised of only one turn for the purpose of calculation. For both core-type and shell-type structures, the same geometry is used but the winding window of the core-type structure contains a total of four windings whereas the shell-type one has only two windings. Fig. 7 shows the magnetic field obtained for the core-type structure inside and outside the winding window, with a current of 1 A flowing inside each winding. Magnetic fields obtained with 2DFEM or with analytical model look identical. To confirm this observation, total energies of these planes were calculated and compared in Table I. The maximum difference on the calculation results between the two different methods is below 1%. From these energy values, total leakage inductance was calculated analytically on each structure. Also, the newly developed "2x2D" model was compared to classical 1D method using magnetic field calculated as presented in II.A.1 on the one hand, and to 3DFEM results on the other hand. This comparison is available in Table II. The "2x2D" method gives results very close to 3DFEM ones with less than 2% deviation, whereas one-dimensional approach overestimates leakage inductance by more than 10% in this case.

B. Frequency-dependent resistance

To be able to compare the calculation results achieved by the analytical models and the 2D finite element modeling, a typical transformer structure was considered. This twodimensional structure is located inside a winding window and constituted of two identical windings. Basic conductors are strands with 0.5 mm diameter and 33 µm enamel thickness. A winding is made up of 1 000 strands organized as shown in Fig. 8, with 10 strands per horizontal row. Horizontal distances between core and windings and between windings are all equal to 5 mm. Finally, the winding height wh is equal to 56.6 mm and the vertical gap eh is fixed to obtain the desired porosity factor following (29) where Ds and ds are the strand diameter with and without enamel. Even if this geometry may not be representative of a Litz wire winding of a medium frequency transformer, it is a good test case to compare models with each other. Moreover, the limited number of strands allows a good meshing of the geometry for finite element modeling by using an eddy current solver (in this case, ten triangles along a strand diameter) with reasonable computation time. It should be noted that in this configuration, current was considered equally divided between all strands, which corresponds to an ideal twisting pattern. Fig. 9 presents the obtained results in all considered cases of porosity factors. The resistance factor was calculated with each of the models and FEM modeling for various frequencies (corresponding to ratio diameter to skin depth ds/δ) in the range of 100 Hz -1 MHz. For asymptotic and Albach, 2D refers to [START_REF] Wojda | Winding resistance of Litz-wire and multi-strand inductors[END_REF] and ( 27) being used with a two-dimensional magnetic field calculation method, whereas 1D refers to ( 23) and (28).

From these results, Dowell's model is always closer to FEM one than Ferreira's or Reatti's. Moreover, asymptotic model gives the same results as Albach's one at low frequencies whereas it is completely false at higher frequencies. Therefore, Dowell and Albach look as the best performing models.

To compare them in a better way, two main zones can be considered: before and after a ratio of diameter to skin depth of 2. The mean absolute error was calculated for each model in both zones and results are presented in Fig. 10. Dowell's model presents an error increasing consistently with the reduction of porosity factor in each zone. However, Albach's model has an error quite independent of the porosity factor. This error remains low at low frequencies, particularly for Albach 2D, but it is much higher at high frequencies.

To complement these results, it should be mentioned that the calculation time is about 0.5 ms for all one-dimensional models and about 25 ms in the case of two-dimensional models.

IV. DISCUSSION

A. Magnetic fields and leakage inductance

For magnetic field calculation, results show very good agreement between the newly developed "2x2D" method and 2D finite element modeling. Magnetic fields and plane energies are almost identical in the considered case.

Moreover, leakage inductance calculation based on this "2x2D" method presents a good accuracy when comparing with finite element modeling for both core-type and shell-type structures. It means that the decomposition of the structure in multiple characteristic planes is a good hypothesis.

However, the accuracy of this approach is limited inside the winding window because the calculation would require an infinity of images to make the error close to zero. This limitation is highlighted by a small difference of 0.41% between numerical simulation and plane energy in Table I, whereas leg side results are almost identical with a difference of 0.01%.

Winding window energy is not the only cause of difference between leakage inductance calculated analytically and numerically. The other one is the two-dimensional decomposition. In fact, the associated depth of each plane is inaccurate because the geometry varies along this depth, particularly in corner region of windings as it can be seen on Fig. 6. It explains why there is more difference for leakage inductance than for planes energies. Nevertheless, accuracy (below 2%) can be considered as very satisfying, especially when considering that calculation time has been reduced by more than 100 times. However, one-dimensional calculation is still faster because there is no discretization at all, but at the expense of a much lower accuracy. So the newly developed "2x2D" model is a good compromise between accuracy and calculation time.

This model is also a very good solution to compute magnetic field only in a specific area without having to calculate it over the whole geometry, as it is the case in finite element modeling.

Finally, maybe the most important advantage of this model is its adaptability. The method used is very general and applicable to almost any transformer or magnetic device geometry, provided that it can be decomposed into multiple two-dimensional structures (symmetries). There is no need for windings to be concentric, and it can be easily adapted to toroidal structures by using the two-dimensional decomposition presented in [START_REF] Prieto | Study of toroidal transformers by means of 2D approaches[END_REF].

B. Frequency-dependent resistance 1) General comments

In the case of frequency-dependent resistance calculation methods, the results show a great disparity between models on accuracy levels and frequency behavior. First of all, three main zones can be considered for resistance factor curve obtained by FEM, in black straight line on Fig. 9. First zone corresponds to a resistance factor close to one which means that there is no significant proximity or skin effect. Then, by increasing the frequency, a second zone appears with a first slope corresponding to proximity effect. This transition can occur at different ratio ds/δ depending on the geometry. Here, the geometry of the winding does not vary so much between the different cases so the transition always occur around ds/δ = 0.7. Finally, the last zone corresponds to another slope, lower than the previous one, combining proximity and skin effects. This transition will always occur around ds/δ = 2 because of skin effect. These three zones and their associated boundary frequencies have been identified in [START_REF] Wojda | Winding resistance of Litz-wire and multi-strand inductors[END_REF] and [START_REF] Wojda | Winding resistance and power loss of inductors with Litz and solid-round wires[END_REF]. Actually, skin effect is already present before this limit but is not high enough to have a significant impact on the resistance factor. The slope in this zone is lower because the skin effect acts as a shielding of the strands and prevents the penetration of external magnetic field inside them, thus confining proximity effect in a reduced portion of the conductors. In other words, magnetic field induced by eddy currents is not negligible anymore in this zone and reduce drastically the static magnetic fields inside conductors.

Before discussing the results obtained for each model, it should be noted that, in the case of medium frequency transformer, the accuracy level in the zones before ds/δ = 2 is much more important that the one obtained after. In fact, a proper design must have a low resistance factor, and this case usually corresponds to the situation where skin effect is negligible if an important number of conductors has to be considered. Therefore, the inaccuracy that can occur in the high frequency zone will mostly impact designs that wouldn't have been retained (in a first design step) because their losses would have been too high even in the low frequency region.

However, it is still relevant to be accurate enough in this zone (after ds/δ = 2), in a second design step, for possible highfrequency harmonics of the considered current, and ensure that these harmonics will not cause too many additional losses.

2) Models comparison

Ferreira's model always overestimates resistance factor. This is because the porosity factor is not considered for this model and therefore the gaps between conductors are ignored as if the vertical dimension was completely filled with conductors. It results in an overestimation of proximity losses. This is why this model will not be further investigated.

Dowell's model is almost exact for a porosity factor of 0.78. In this case, the distance between windings and core eh is zero and the porosity factor is only due to gaps between strands inside windings. Therefore, the one-dimensional hypothesis for magnetic field is valid. The use of equivalent foil to represent strands seems to be in accordance with Dowell's hypothesis regarding the level of accuracy obtained (less than 2%). The model remains exact even at high frequency where skin effect is important because Dowell's model is not influenced by magnetic field of nearby conductors. In fact, the equation is obtained by setting boundary conditions on both sides of the foil and solving Maxwell equations. The field on the side of the foil is independent of the frequency so boundary conditions remain valid on the whole frequency range and the equation of the resistance factor is also valid. This is why Dowell's model gives the best results in the case of one-dimensional field. However, when the distance between core and windings increases, one-dimensional hypothesis is no longer valid which explains why Dowell's model accuracy decreases in both zones before and after significant skin effect.

Reatti and Kazimierczuk model is similar to Ferreira's model except that it takes into account the porosity factor. This is why the model behaves like Dowell's model. However, its accuracy at high frequency is limited because the eddy currents due to magnetic field induced by eddy currents of nearby conductors is not taken into account. Thus this model is less accurate than Dowell's one and will not be further investigated.

The main disadvantage of the asymptotic model is its behavior in the high frequency zone. Even if the resistance factor is accurate at low frequency, skin effect is not considered at all. Therefore, no shielding effect is preventing eddy current development inside conductors at high frequency, and the losses are greatly overestimated. Everything happens as if proximity effect was independent of frequency, and this is why the asymptotic behavior occurs. The disadvantage is the same for one-dimensional or twodimensional versions, so this model will not be further investigated.

For a porosity factor of 0.78, Albach's model has the same behavior as Reatti and Kazimierczuk. This is because in this case, formulation of equations are almost equivalent. It means Albach's model suffers from the same limitation at high frequency, not taking into account mutual influence between conductors due to magnetic field produced by eddy currents. However, when the porosity factor decreases, Albach's model accuracy remains almost constant. This is because the porosity factor is not directly used, but the proportion of conductive material inside winding η is considered instead. Overall, this is why this model gives the best results in the cases of low porosity factor, with a very good accuracy at low frequency and a contained error at high frequency. The one-dimensional version of the model has an error of less than 10% at low frequency and below 100% at high frequency, whereas the two-dimensional version has an error below 2% at low frequency and 60% at high frequency. In fact, both versions have the same accuracy for a porosity factor of 0.78 because one-dimensional hypothesis is valid in this case. For lower porosity factor, the two-dimensional model is logically more accurate because the magnetic field is two-dimensional. However, if calculation time is critical, the one-dimensional version is still a solution for this type of geometry.

All these observations explain why Albach's model is considered as the best candidate for medium frequency transformer design. It provides a good accuracy and a great resilience to porosity factor variations from a design to another, with no additional calculation time in its onedimensional version. Actually, it is even faster than Dowell's model as Bessel functions inputs depend only on strand diameter and frequency, whereas the inputs of trigonometric functions of Dowell's model include the porosity factor and therefore depend on the transformer geometry. This allows a reduction of the calculation time when applied to a great number of geometries.

Moreover, the model is very adaptable. As long as the magnetic field is known for a given geometry, the proximity losses can be calculated and a resistance factor formula can be derived. This hypothesis is still satisfied for configurations without concentric windings, such as toroidal structures or even inductors where there is no secondary and therefore no Ampère-Turns compensation.

Nevertheless, this model has some limitations. First one, already discussed, is its reduced accuracy at higher frequency. This problem could be limited if equation (26) was used instead of its simplified version (27) and an iterative process was followed [START_REF] Albach | Two-dimensional calculation of winding losses in transformers[END_REF]. This last option would complicate the model implementation and it would have a high negative impact on the global calculation time in a first design step. A possible way would be to use the improved version only for a very limited number of designs (second design step). It leads us to another downside of this model which is the necessity of calculating magnetic field inside windings with enough accuracy. The method presented in this article to calculate magnetic field is very adaptable, but at high frequency, the windings cannot be homogenized as a macroscopic rectangular conductor and each strand contribution must be considered. Therefore the computation time would depend on the number of strands, which is a real inconvenience. A way to solve this last problem would be to combine Albach's model to determine equivalent complex magnetic permeability and electrical conductivity for a given winding structure, following a method such as the one presented in [START_REF] Phung | High-frequency proximity losses determination for rectangular cross-section conductors[END_REF]. It would correct the accuracy at higher frequency, but to the detriment of a more complex model and an important calculation time. This could be a solution to calculate a precise frequency-dependent resistance for specific designs, but again, it is not suitable for integration inside an optimization process.

V. CONCLUSION

In this article, a review of magnetic field, leakage inductance and frequency-dependent resistance calculation methods has been done. For each subject, new methods of calculation were developed. The two-dimensional ("2x2D") method for magnetic field calculation revealed its accuracy with low calculation time, and leads to precise leakage inductance evaluation in the case of usual core-type and shell-type transformer structures.

Concerning the frequency-dependent resistance calculation method, the newly developed model based on Albach's theory is more accurate, more robust and faster than reviewed models.

Overall, newly developed models are already an improvement from reviewed ones, and are very promising due to their possibility of further improvements: adaptability to new structures (toroidal, inductors) and increased accuracy for high frequency. In this sense, they are very good candidates for magnetic devices design and optimization. Implementation of these tools is in progress. From these relations and ( 18) and (25):

                                                                      
                                                                     

Fig. 1 .

 1 Fig. 1. Concentric windings inside a winding window and associated onedimensional magnetic field distribution.

Fig. 2 .

 2 Fig. 2. Schematic of a round conductor and magnetic field distribution of a 1 mm radius conductor carrying 1 A.

Fig. 3 .

 3 Fig. 3. Schematic of a rectangular conductor and magnetic field distribution of a 1 mm  2 mm conductor carrying 1 A.

Fig. 4 .

 4 Fig. 4. Image method for conductors next to a high permeability material surface.

Fig. 5 .

 5 Fig. 5. Image method for conductors inside a winding window.

Fig. 6 .

 6 Fig. 6. a) Core-type structure. b) Planes and associated depths for core-type structure (top view). c) Shell-type structure. d) Planes and associated depths for shell-type structure (top-view).

Fig. 7 .

 7 Fig. 7. Magnetic field distribution for: a) Winding window, analytical model; b) Leg side, analytical model; c) Winding window, 2DFEM; d) Leg side, 2DFEM.

Fig. 8 .

 8 Fig. 8. Geometry definition of windings and window used for test cases of frequency-dependent resistance calculation.

Fig. 9 .

 9 Fig. 9. Frequency-dependent resistance factor obtained for all tested models on test cases with different porosity factor (one per graph).

Fig. 10 .

 10 Fig.10. Mean Absolute Error (MAE) of three best-performing models. First graph is MAE only for frequencies where skin effect is negligible. Second graph is MAE only for frequencies where skin effect is not negligible.

  APPENDIX ADemonstration of one-dimensional proximity factor for asymptotic model.



  model in case of homogeneous magnetic field.Series Fourier of any magnetic field in polar coordinates on the surface of the conductor:If magnetic field is constant over the strand:Then by injecting this result in (26), one can obtain (27).APPENDIX CDemonstration of equality between FS and fP of Albach and Ferreira models.

  Same method is applied for (20) and (28):

TABLE II COMPARISON

 II OF LEAKAGE INDUCTANCE CALCULATION METHODS

	Structure	Method	Leakage Inductance	Difference Time
		Analytical 1D	1.6775.10 -7 H	+13.44%	0.3 ms
		Analytical 22D 1.4691.10 -7 H	-0.54%	230 ms
	Core -	FEM 3D	1.4896.10 -7 H	+0.85%	28 s
	type	8k tetrahedrons
		FEM 3D	1.4770.10 -7 H	0%	3066 s
		1M tet.
		Analytical 1D	8.3776.10 -8 H	+13.84%	0.3 ms
		Analytical 22D 7.2434.10 -8 H	-1.57%	100 ms
	Shell -	FEM 3D	7.4310.10 -8 H	+0.98%	21 s
	type	5k tet.
		FEM 3D	7.3590.10 -8 H	0%	1139 s
		0.7M tet.
	Comparison of leakage inductance values obtained with analytical and
	numerical methods. FEM 3D with 0.7 -1 million tetrahedrons is taken as
	reference.

ACKNOWLEDGMENT

This work was supported by a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" Program (ANE-ITE-002-01).