Marc Rébillat
email: marc.rebillat@limsi.fr

Romain Hennequin

Etienne Corteel

Brian F G Katz

Étienne Corteel

Identification of cascade of Hammerstein models for the description of non-linearities in vibrating devices

Keywords: Non-linear system identification, Cascade of Hammerstein models, Exponential sine sweep method, Harmonic distortion

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Vibratory phenomena are usually assumed to be linear. However, many vibrating systems are subject to non-linear behaviours, such as loudspeakers [START_REF] Quaegebeur | Nonlinear vibrations of loudspeaker-like structures[END_REF], musical instruments [START_REF] Gilbert | A simulation tool for brassiness studies (l)[END_REF] and vibrating plates [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF]. Even wave propagation in air is not completely linear [START_REF] Gee | Measurement and prediction of nonlinearity in outdoor propagation of periodic signals[END_REF]. The study of these non-linearities is thus of great importance in order to model these devices and phenomena or to justify their "linearity".

Identification of non-linear systems requires measurements or estimation of model's structural elements from a finite set of input/output data [START_REF] Chen | Modeling and identification of parallel nonlinear systems: Structural classification and parameter estimation methods[END_REF]. Classical linear measurement methods [START_REF] Muller | Transfer-function measurement with sweeps[END_REF][START_REF] Stan | Comparison of different impulse response measurement techniques[END_REF] capture only the linear behaviour of the system under study. Traditional non-linear measurement methods [START_REF] Czerwinski | Multitone testing of sound system components some results and conclusions -part 2: Modeling and application[END_REF] give total harmonic distortion (THD), harmonic distortion of order n (HD n ) or inter-modulation products (IMP). These quantities are measured using pure tones at a given amplitude and frequency. They do not describe non-linearities themselves but only some of their effects for arbitrary excitations. Moreover, experimental processes involved in those methods are very time-consuming if a wide range of frequencies and amplitudes is to be considered. There is thus a real need for rapid model based procedures to measure non-linearities.

Non-linear systems can be classically represented by Volterra series [START_REF] Hasler | Phénomènes non linéaires, Chapitre 3 : Séries de Volterra[END_REF] or by "Sandwich" structures [START_REF] Chen | Modeling and identification of parallel nonlinear systems: Structural classification and parameter estimation methods[END_REF]. The cascade of Hammerstein models [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF] is a subclass of those models and can be used to exactly represent systems having diagonal Volterra Kernels. This model is composed of N branches in parallel. Each branch comprised a static polynomial non-linearity followed by a linear filter.

A simple method which makes it possible to quickly access the structural elements of a cascade of Hammerstein models is presented in this paper. This method is based on the method initially proposed by Farina [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF]. Exponential sine sweeps are used as input signals, and allow for the temporal separation of the different orders of non-linearity [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Craven | Practical adaptive room and loudspeaker equaliser for HI-FI use[END_REF][START_REF] Norcross | A survey on the effects on nonlinearity on various types of transfer function measurements[END_REF][START_REF] Griesinger | Beyond MLS -Occupied hall measurement with FFT techniques[END_REF]. Structural elements of the assumed model are then extracted from only one response of the system. The method is validated on a simulated system and the influence of the different parameters is shown.

As transducers are most often the least linear part in the audio chain, knowing their non-linear behaviour is very important. Two major types of transducers exists: electrodynamic loudspeakers [START_REF] Klippel | Tutorial: Loudspeaker nonlinearities -causes, parameters, symptoms[END_REF] and panel-based transducers (distributed mode loudspeakers [START_REF] Harris | Introduction to distributed mode loudspeakers (DML) with first-order behavioural modelling[END_REF] or multi-actuator panels [START_REF] Kuster | Structural and acoustic analysis of multiactuator panels[END_REF]). In electrodynamic loudspeakers, a motor converts the electrical signal into motion and makes a cone vibrate. The piston-like movement of the cone generates the sound field. In panel-based loudspeakers, a motor is also used but transmits its motion to a light and stiff panel. The flexural waves travelling in the panel then generate the sound field.

In both loudspeaker types, the motor induces non-linearities because of non-uniform magnetic field, Eddy currents and variations of the electrical inductance with displacements [START_REF] Klippel | Tutorial: Loudspeaker nonlinearities -causes, parameters, symptoms[END_REF][START_REF] Colloms | Distortion mechanisms of distributed mode (DM) panel loudspeakers[END_REF]. In electrodynamic loudspeakers, significant excursion can induce non-linear bending in the cone and a non-linear behaviour of the suspensions [START_REF] Quaegebeur | Nonlinear vibrations of loudspeaker-like structures[END_REF]. In panel-based loudspeakers, large amplitude displacements occur in the plate near the exciter position. In this case the propagation of flexural waves [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF] and the strain/stress relation of the material which compose the plate [START_REF] Daniel | Fabrication, testing and analysis of composite sandwich beams[END_REF] can be non-linear.

In the literature, electrodynamic loudspeakers have been greatly studied from a non-linear point of view. To represent their non-linear behaviour, different physical models have been built. Their formulation was either completely analytical [START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a Volterra series expansion[END_REF][START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF] or based on the finite element method [START_REF] Tsuchiya | Finite element simulation of non-linear acoustic generation in a horn loudspeaker[END_REF]. In Ref. [START_REF] Klippel | Nonlinear system identification for horn loudspeakers[END_REF], Klippel proposed to reduce the Volterra series expansion to a "Sandwich" model and identified its parameters from measurements using the method presented in [START_REF] Baumgartner | Complete identification of a class of nonlinear-systems from steady-state frequencyresponse[END_REF]. In [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF], it is suggested to use simpler models, i.e. cascade of Hammerstein models, to model and identify different audio systems, including acoustical transducers.

In the present paper, two different acoustical transducers (an electrodynamic one and a panel-based one) are studied experimentally under the assumption that they can be modelled using cascade of Hammerstein models. Their models are completely identified using the previously introduced method. THD and HD n at different frequencies and amplitudes of the input signal are evaluated for these transducers by classical means and compared to predictions made using the identified Kernels. The agreement between the results given by the two methods is very good for a wide range of amplitudes.

After reviewing how to model and to measure non-linearities (section 2), the mathematical foundations of the current method are presented (section 3). Implementation of the method is then described (section 4) and validated on a simulated system (section 5). Acoustical transducers are identified using the previous method (section 6) and the resulting Kernels are used successfully to predict the harmonic distortion generated by the two transducers (section 7).

Modelling and measurements of non-linearities

An overview of existing models and measurement methods of non-linearities is given in this section. Only Single-Input/Single-Output (SISO) time-invariant causal non-linear systems without continuous component will be considered here.

Volterra series

Volterra series [START_REF] Hasler | Phénomènes non linéaires, Chapitre 3 : Séries de Volterra[END_REF] enables one to express the relationship between the non-linear system input e(t) and output s(t) as a series of multiple convolution integrals :

s(t) = +∞ k=1 +∞ 0 ... +∞ 0 v k (τ 1 , ..., τ k )e(t -τ 1 )...e(t -τ k )dτ 1 ...dτ k (1)
The functions {v k (t 1 , .., t k )} k∈N * are called Volterra Kernels and completely characterize the system. Volterra models can then be seen as a generalization of the simple convolution operator used for linear systems. Such models represent exactly any non-linear "analytical" system [START_REF] Palm | On representation and approximation of nonlinear systems[END_REF], and approximate any non-linear system with a "fading memory" [START_REF] Boyd | Fading memory and the problem of approximating nonlinear operators with Volterra series[END_REF]. Measurement methods exists to identify the first two or three terms of Volterra series [START_REF] Boyd | Measuring Volterra Kernels[END_REF][START_REF] Bicken | Frequency domain Kernel estimation for 2nd-order Volterra models using random multi-tone excitation[END_REF][START_REF] Prazenica | Volterra Kernel identification using triangular wavelets[END_REF][START_REF] Kayran | Nonlinear system identification using deterministic multilevel sequences[END_REF]. These experimental methods are time consuming because they require many measurements. Moreover the difficult physical interpretation of the different terms of the Volterra series limits its use [START_REF] Voishvillo | Graphing, interpretation, and comparison of results of loudspeaker nonlinear distortion measurements[END_REF].

Sandwich approach

Another approach to non-linear system identification is to assume that systems have a given block-structure. Following the "Sandwich" approach [START_REF] Chen | Modeling and identification of parallel nonlinear systems: Structural classification and parameter estimation methods[END_REF], a non-linear system is represented as M parallel branches composed of three elements in series: a static non-linear part P m (.) sandwiched between two linear parts l l m (t) and l r m (t). Such systems are a subclass of Volterra systems. It can be shown that any continuous non-linear system can be approximated by such a model [START_REF] Palm | On representation and approximation of nonlinear systems, part II: Discrete time[END_REF].

To identify such structures, the form of the static non-linear part can be assumed and the two unknown linear parts can be estimated from measurements at different frequencies and amplitudes [START_REF] Baumgartner | Complete identification of a class of nonlinear-systems from steady-state frequencyresponse[END_REF]. This leads to only a discrete knowledge of L l m (f ) and L r m (f ) in the frequency domain and remains a long experimental task. In Ref. [START_REF] Abel | A technique for nonlinear system measurement[END_REF], Abel et al. proposed another method to identify such structures. Unfortunately, this general method cannot be used successfully in practice due to numerical instabilities.

Cascade of Hammerstein models

In a cascade of Hammerstein models [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF], each branch is composed of one non-linear static polynomial element followed by a linear one h n (t), as shown in Fig. 1. Mathematically, the relation between the input e(t) and the output s(t) of such a system is given by Eq. ( 2), where * denotes the convolution.

s(t) = N n=1 h n * e n (t) (2)
In this model, each impulse response h n (t) is convolved with the input signal elevated to its n th power and the output s(t) is the sum of these convolutions. The first impulse response h 1 (t) represents the linear response of the system. The other impulse responses {h n (t)} n∈{2...N } model the non-linearities.

The family {h n (t)} n∈{1...N } will be referred to as the Kernels of the model. These Kernels are assumed to be integrable. Any cascade of Hammerstein models is fully represented by its Kernels.

It can easily be shown from Eqs. ( 1) and ( 2) that cascade of Hammerstein models correspond to Volterra models having diagonal Kernels in the temporal domain, as in Eq. ( 3), where δ(t) represents the Dirac distribution. This non-linear model is thus referred in the literature as a diagonal Volterra model [START_REF] Farina | Non-linear convolution: a new approach for the auralization of distorting systems[END_REF], but also as a cascade of Hammerstein models [START_REF] Chen | Modeling and identification of parallel nonlinear systems: Structural classification and parameter estimation methods[END_REF] or Uryson model [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF].

∀(τ 1 , . . . , τ k ) v k (τ 1 , . . . , τ k ) = h k (τ 1 )δ(τ 1 ) . . . δ(τ k ) (3)
As can be seen in Eq. ( 2), cascade of Hammerstein models are linear in the parameters to be estimated, i.e. the output of the system is a linear combination of the Kernels {h n (t)} n∈{1...N } . A naive approach is to identify the model using a classical least square method, as proposed for general Volterra systems in [START_REF] Reed | Identification of discrete Volterra series using maximum length sequences[END_REF]. Thus the mean squared error between the actual output of the system y(t) and the output of the estimated model s(t) given in Eq. ( 2) can be minimized with respect to the coefficients of h 1 (t), h 2 (t), . . . , h N (t) and the solution is given by: argmin h1(t),h2(t),...,hN (t) t y(t) -s(t) 2 (4) However, the least square method requires the inversion of a M N × M N matrix, where N is the order of the system under test and M is the length of the impulse responses h n (t) in samples. This matrix can be very ill-conditioned since it is generated from the exponent (until N ) of the input signal. This results in important errors in parameters estimation especially in noisy conditions. Moreover the computation of the matrix from the input signal and of the inverse of the matrix is computationally costly and limits in practical case the memory M of the system. Some numerical methods are however available to limit these points (see for example [START_REF] Golub | Numerical methods for solving linear least squares problems[END_REF][START_REF] Chen | Orthogonal least-squares methods and their application to non-linear systemidentification[END_REF]), and to overcome these drawbacks alternative methods have also been developed.

Gallman [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF] and Hawksford [START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF] proposed a method to estimate the elements of a cascade of Hammerstein models using Gaussian noise at different amplitudes as inputs. The employed estimation procedures are strongly based on the knowledge of the order of non-linearity of the polynomial expansion, which is unknown in practical cases. Moreover, these methods to identify the Kernels from the measurements are also computationally costly.

Farina proposed another method using sine sweeps with frequency varying exponentially with time [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF]. An upper bound of the order of non-linearity of the model has to be assumed. This method allows only for the separation of the different orders of non-linearity and not for the complete identification of the Kernels of the system. Recently, Novák et al. [START_REF] Novák | Modeling of nonlinear audio systems using swept-sine signals: application to audio effects[END_REF] have identified Kernels from the contributions of the different orders of non-linearity using a least mean square minimization procedure. No results are provided to judge the influence of the different parameters on its performances.

The method proposed here gives direct mathematical access to all the Kernels {h n (t)} n∈{1...N } from the contributions of the different orders of non-linearity obtained as in Ref. [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF]. The main advantage of the proposed approach over the least squared based technique, besides conditioning and computational problems, is that it provides a direct evaluation of the N impulse responses h n (t) of the system. The foundations and the key implementation of this simple method are explained in detail. The whole procedure is validated on a simulated system and on two real systems.

Mathematical foundations of the method

Mathematical foundations of the method used for direct identification of the elements of a cascade of Hammerstein models are given in this section. This is based on the procedure initially proposed by Farina [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF].

A cascade of Hammerstein models fed with sine sweeps

To experimentally cover the frequency range on which the system under study is to be identified, cosines with time-varying frequencies are interesting signals. Eq. ( 5) defines such a signal.

∀t ∈ R e(t) = cos(Φ(t)) (5)
If e(t) is the input of the cascade of Hammerstein models, the output of the non-linear block of the i th branch will have the form of Eq. ( 6), as can be seen in Fig. 1.

e i (t) = cos i [Φ(t)] (6)
Using Chebyshev polynomials, e i (t) is rewritten in Eq. ( 7) as a linear function of {cos[kΦ(t)]} k∈ [1,i] . Details of the computation of the matrix C are provided in the appendix.

∀i ∈ {1 . . . N } cos i [Φ(t)] = i k=0 C(i, k) cos[kΦ(t)] (7)

Exponential sine sweeps

When the instantaneous frequency of e(t) is increasing exponentially from

f 1 to f 2 (f 1 , f 2 > 0)
in a time T , such a signal is referred to as an "exponential sine sweep" [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Norcross | A survey on the effects on nonlinearity on various types of transfer function measurements[END_REF] and its instantaneous phase is given by :

∀t ∈ R Φ(t) = 2π f 1 T ln f2 f1 (e t T ln f 2 f 1 -1) -π/2 (8)
The corresponding instantaneous frequency of e(t) is :

∀t ∈ R f (t) = Φ ′ (t) 2π = f 1 e t T ln f 2 f 1 (9)
Thus f (0) = f 1 and f (T ) = f 2 . The frequency range [f 1 , f 2 ] corresponds to the band of interest of the system under test.

Fundamental phase property

From Eq. ( 8), it can be shown that this type of signal exhibits the following phase property :

∀k ∈ N * , ∀t ∈ R kΦ(t) = Φ(t + T ln k ln f2 f1 ) -(k -1)( π 2 + 2π f 1 T ln f2 f1 ) (10)
By choosing

T m = (2mπ - π 2 ) ln f 2 /f 1 2πf 1
with m ∈ N * , the second term in Eq. ( 7) becomes a multiple of 2π and one obtains Eq. ( 11) which is another way to express the k th term of the linearisation presented in Eq. [START_REF] Stan | Comparison of different impulse response measurement techniques[END_REF].

∀k ∈ N * , cos(kΦ(t)) = cos(Φ(t + ∆t k )) with ∆t k = T m ln k ln (f 2 /f 1 ) (11)
For any T m -long logarithmic sweep, multiplying the phase by a factor k results in the same signal, but in advance in the time domain by ∆t k . As can be seen from Eq. ( 11), this time advance depends only on the sweep parameters T m , f 1 , f 2 and on k. In Refs. [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Morise | Warped-TSP: An acoustic measurement signal robust to background noise and harmonic distortion[END_REF], similar time advances were obtained using different arguments.

The fact that T m must take only a discrete set of values to ensure the fundamental phase property Eq. ( 11) has been first shown in Ref. [START_REF] Novák | Modeling of nonlinear audio systems using swept-sine signals: application to audio effects[END_REF] but is mathematically demonstrated here for the first time.

One should note that e(t) has been designed for all t with its instantaneous frequency increasing from f 1 to f 2 between t = 0 and t = T . In practice, signals are defined only on [0, T ]. Thus the phase property is not valid on the whole support of the function. The phase property becomes false when the instantaneous frequency of cos[kΦ(t)] is outside the frequency range of interest (i.e. [f 1 , f 2 ]).

Inverse convolution

Using the Eqs. ( 11) and ( 2), one obtains :

s(t) = N n=1 g n * e(t + ∆t n ) + K (12)
with:

g n (t) = N k=1 C(k, n)h k (t) and K = N n=1 C(n, 0) +∞ -∞ h n (t)dt (13)
In Eq. ( 13), g n (t) corresponds to the contribution of the different Kernels to the n th harmonic. K is the global continuous component resulting from the continuous components of the different Kernels. As the Kernels are assumed to be integrable, K is correctly defined.

In order to separately identify each Kernel of the cascade of Hammerstein models, a signal y(t) which looks like an inverse in the convolution sense of e(t) is needed. Unfortunately, such an inverse does not necessarily exist mathematically. However, a band-limited inverse y(t) can easily be defined such that it satisfies the relation [START_REF] Griesinger | Beyond MLS -Occupied hall measurement with FFT techniques[END_REF] with sinc(x) = sin(πx)/πx.

y * e(t) = sinc(2f 2 t) -sinc(2f 1 t) = d(t) (14)
d(t) can be seen as a band-limited Dirac Function, since its Fourier transform is

1 [-f2,-f1]∩[f1,f2] (f ).
Then, Y (f ), the Fourier transform of the inverse filter y(t) can be built in the frequency domain using Eq. [START_REF] Klippel | Tutorial: Loudspeaker nonlinearities -causes, parameters, symptoms[END_REF], where E * (f ) is the complex conjugate of E(f ).

Y (f ) = 1 E(f ) 1 [-f2,-f1]∩[f1,f2] (f ) ≃ E * (f ) |E(f )| 2 + ǫ(f ) (15)
In practice, the filter Y (f ) should be built by replacing the discontinuous function 1 [-f2,-f1]∩[f1,f2] (f ) by a function which ensures a smoother transition between the two frequency domains and thus generates less unwanted side effects in the time domain.

ǫ(f ) = β ×T (f ) is a frequency-dependent real parameter chosen as equal to 0 in the bandwidth and as having a large value β outside of it, with a continuous transition between the two domains. In the following, a weight

β = f s 2 0 |Y (f )| 2 df
, which corresponds to the energy of the signal to be inverted, has been chosen. In practice, transitions between the two domains can be simple linear functions or C ∞ Gevrey functions. An example of such a function defining a C ∞ transition between T (f a ) = 0 and T (f b ) = 1 is:

∀f ∈]f a , f b [ T (f ) = 1 2 1 + tanh 1 f a -f + 1 f b -f (16)
The application of this procedure leads to y(t) that can be considered as an inverse of e(t) in the sense of convolution in the frequency range [f 1 , f 2 ].

Kernel Identification in the temporal domain

After convolving the output of the cascade of Hammerstein models s(t) given in Eq. ( 12) with y(t), one obtains :

y * s(t) = N i=1 d * g i (t + ∆t n ) = N i=1 gi (t + ∆t n ) ( 17
)
where gi (t) corresponds to g i (t) convolved with d(t), i.e. to g i (t) filtered by a bandpass filter in the frequency band 12), has thus been filtered by d(t). Moreover, if the system under study has no significant behaviour outside of [f 1 , f 2 ], then gi (t) = g i (t).

[f 1 , f 2 ]. The constant K, present in Eq. (
In Fig. 2, y * s(t) is represented. Because ∆t n ∝ ln(n) and f 2 > f 1 , the higher the order of linearity n the more in advance the corresponding gn (t) will be. Thus, if T m is long enough, the different gn (t) will not overlap. They are then easy to separate by windowing in the time domain. The separation of the contribution of the different orders of non-linearity by using exponential sweeps, which is mathematically demonstrated here, is already experimentally well known in the audio community [START_REF] Craven | Practical adaptive room and loudspeaker equaliser for HI-FI use[END_REF][START_REF] Norcross | A survey on the effects on nonlinearity on various types of transfer function measurements[END_REF][START_REF] Griesinger | Beyond MLS -Occupied hall measurement with FFT techniques[END_REF].

Next, using Eq. ( 18), the family { hn (t)} n∈ [1,N ] of the Kernels of the cascade of Hammerstein models describing the behaviour of the system in the frequency band [f 1 , f 2 ] can be fully extracted.

   h1 (t)

. . . A c is the Chebyshev matrix defined in the appendix without its first column and its first row and (.) T stands for matrix transposition. The first column and the first row of matrix A have been removed here as there is no continuous components here.

hN (t)    = A T c    g1 (t) . . . gN (t)    (18)
As hn (t) = d * h n (t), if the system under study has no significant behaviour outside of [f 1 , f 2 ], then hn (t) = h n (t). In most vibration application, systems are designed for a given frequency band (typically [20Hz, 20kHz] for audio applications). The border frequencies f 1 and f 2 can thus generally be selected to identify the real Kernels h i (t). If it is not possible, Kernels are identified between f 1 and f 2 and thus are only a band-limited version of the real Kernels.

Practical implementation

In this section, the practical discrete-time implementation of the method presented in section 3 is described.

Overview of the method

In Fig. 3, a global overview of the procedure is given. It can be decomposed in the following steps :

1. Design of the input sweep e(t) using Eq. ( 8). The choices of f 1 , f 2 and T are discussed in section 4.2.

2. Playing e(t) and recording s(t). The sampling frequency f s must be chosen to avoid any aliasing effects caused by the digital to analog converter in the frequency range of interest

[f 1 , f 2 ].
3. Generation of the inverse filter y(t) according to Eq. [START_REF] Klippel | Tutorial: Loudspeaker nonlinearities -causes, parameters, symptoms[END_REF]. A convenient way to implement this filter is described in section 3.4.

4. Convolution of the output signal s(t) with the inverse filter y(t) as in Eq. ( 17). This can be done in the frequency domain with a sufficient number of points to avoid temporal aliasing.

5.

Windowing in the temporal domain to obtain { gk (t + ∆t k )} k∈[1,N ] (cf. Fig. 2). Rectangular windows can be chosen to separate the different orders of non-linearity. N is the highest desired order in the cascade Hammerstein model. Methods to choose N and its influence are shown in sections 5.3 and 5.4.

6. Temporal shift of the different orders of non-linearity to recover { gk (t)} k∈ [1,N ] . A shift of a non-integer number of samples can be performed with a phase shift in the frequency domain.

Multiplication with

A T c to access { hk (t)} k∈[1,N ]
, according to Eq. [START_REF] Colloms | Distortion mechanisms of distributed mode (DM) panel loudspeakers[END_REF]. A c is the Chebyshev matrix defined in the appendix without the first column and the first row. The matrix A T c of order 8, which is sufficient in practice, is explicitly given is the appendix.

Choice of the parameters (f 1 , f 2 , T and N )

For satisfactory measurements, the sweep parameters f 1 , f 2 , T and N must be well chosen. These choices must be made considering the following aspects:

• The frequencies f 1 and f 2 must be chosen such that the interesting behaviour of the system under study is in the frequency range

[f 1 , f 2 ].
• The influence of noise on the identification results should also be minimized [START_REF] Morise | Warped-TSP: An acoustic measurement signal robust to background noise and harmonic distortion[END_REF]. By itself, the exponential sweep rejects correctly uncorrelated noise in quiet environment [START_REF] Stan | Comparison of different impulse response measurement techniques[END_REF]. Moreover its energy repartition in frequency is often adapted to the ambient noise [START_REF] Muller | Transfer-function measurement with sweeps[END_REF][START_REF] Griesinger | Beyond MLS -Occupied hall measurement with FFT techniques[END_REF]. The choice of this signal is thus interesting from this point of view as will be seen in section 5.3. If an excellent signal to noise ratio (SNR) is needed, the longer T , the better the SNR after step 4 will be at a given amplitude of the input signal.

• The number of points to be convolved at step 4 is limited by the available computational power. Thus, as T increases, the calculation time will increase. T should not be too large in order to avoid long calculation times.

• N should not be underestimated in order to guarantee good accuracy in identification. The optimal N is reached when is it impossible to extract the corresponding N th impulse response from the background noise. This will be shown in sections 5.3 and 5.4.

• The different peaks g k (t) which appear in the temporal domain after the convolution with the inverse signal (step 4, see Fig. 2) must not overlap each other. The global decay time of the system τ global is an upper bound of the decay times of each order of non-linearity. Parameters f 1 , f 2 and T such that ∆t N -∆t N -1 > τ global will thus avoid overlapping of the different orders of non-linearity [START_REF] Muller | Transfer-function measurement with sweeps[END_REF][START_REF] Majdak | Multiple exponential sweep method for fast measurement of head-related transfer functions[END_REF]. Because ∆t N = T m ln N ln f 2 /f 1 , considering the chosen value for N , T must be chosen to be long enough and f 2 /f 1 not so large in order to respect the previous condition.

Validation of the method

In this section, the proposed method of identification is tested on a simulated cascade of Hammerstein models.

Design of the simulated system

A cascade of Hammerstein models of order N = 4 has been chosen for simulation purposes. This system is fully represented by its 4 Kernels h 1 (t), h 2 (t), h 3 (t) and h 4 (t). For these Kernels, which correspond to linear subsystems, impulse responses of low order ARMA filters (2 poles and 2 zeros) with a roughly 5 ms decay time have been chosen. Parameters of the simulated system are given in Table 1.

n

Identification without noise

The method presented in section 3 and implemented as described in section 4 has been applied here with the parameters given in Table 2 in order to identify the different Kernels of the system.

The magnitude and phase of the frequency responses of the original and estimated Kernels H For a finer analysis, the relative errors in dB, defined in Eq. ( 19), are given for the four estimated Kernels in Fig. 5. As magnitude estimation errors and phase estimation errors are included in this relative error, this is a more compact way to access estimation errors.

1 (f ), H 2 (f ), H 3 (f )

Parameter

Value Normalized value

f 1 20 Hz f 1 /f s = 1.04 × 10 -4 f 2 20 kHz f 2 /f s = 1.04 × 10 -1 f s 192 kHz T 15 s T × f s = 2.88 × 10 6 samples N 4
Table 2: Parameters chosen to identify the simulated system.

γ k (f ) = 20 log 10 H orig k (f ) -H meas k (f ) H orig k (f ) (19)
In Fig. 5, it can be seen that the relative estimation error is lower than -20 dB over a large portion of the frequency range [f 1 , f 2 ]. The consequences of the errors made in magnitude and phase near the border of [f 1 , f 2 ] are clearly visible. To avoid these side effects, the frequency range of the sweep can be chosen larger than the frequency range of interest, depending on the desired accuracy. Near the poles and zeros of the ARMA filters, estimation errors also increase slightly.

The method proposed here gives very good results for identification purposes over a given frequency range, without added noise.

Sensitivity to noise

In a second step, the influence of noise on the estimated Kernels has been studied on the simulated system. A white Gaussian noise at different levels has been added to the output s(t) of the system under study (see Fig. 3). Signal to noise ratios relatively to the input (SNR I ) and before convolution (SNR B ), i.e. before step 4 (Sec. 4 and Fig. 3) are presented in Table 3. Kernel to noise ratios (SNR A ) after step 7 are also given in Table 3. SNR is understood here as the ratio between the root-mean-square (RMS) level of the signal (or of the Kernel for SNR A ) in absence of noise and the RMS level of noise in absence of signal (respectively in absence of Kernel). Signals are recorded at the input of the system for SNR I and at the output of the system for SNR B . Kernels are taken after the complete identification procedure for SNR A . SNR A corresponding to each identified Kernel is given individually. The different Kernels of the system under study have been estimated using the parameters of Table 2 in the different noise conditions. In Table 3, it can be seen that SNR A for N = 1 is 13.6 dB higher than SNR B . This confirms the fact that exponential sine sweeps reject a great part of the uncorrelated noise. Table 3: Signal to noise ratio relatively to the input (SNR I ), before convolution (SNR B ) and after the complete identification procedure (SNR A ). SNR A is indicated relatively to the n th Kernel (n ∈ {1, 2, 3, 4}).

SNR I (dB) SNR B (dB) SNR A (dB) N = 1 N = 2 N = 3 N =
In Fig. 6, the relative estimation errors, according to Eq. ( 19), are given for the different Kernels and for the different noise levels. It can be seen that the noise level has an influence on the quality of the estimation. For each Kernel, when the noise level is decreasing, the relative estimation error is decreasing too.

For a SNR I of 37 dB (SNR B of 15.2 dB), only the estimation of the first Kernel is acceptable. The three other Kernels, which have gains at least 20 dB lower than the first one (see Fig. 4), cannot be estimated correctly in this case. Let's consider in Tab. 3 the SNR A for n ≥ 2. They are all around or below 0 dB. That means that after identification, resulting noises have a RMS level equal to, or higher than, that of identified Kernels. This explains why H 2 (f ), H 3 (f ) and H 4 (f ) cannot be estimated correctly in this case.

When SNR I is 57 dB (SNR B of 35.2 dB), the SNR A corresponding to the second and third Kernels are higher than 10 dB. The second and third Kernels, H 2 (f ) and H 3 (f ), are then correctly estimated. But the last Kernel, with SNR A = -18.6 dB is still not identified.

When the SNR I is 77 dB (SNR B of 55.2 dB), the SNR A corresponding to the fourth Kernel will be 1.4 dB, and the fourth Kernel H 4 (f ) is then finally correctly estimated.

So, with the parameters given in Table 2 applied to the chosen system, it seems that any Kernel can be correctly estimated until its SNR A reaches ≃ 0 dB. Otherwise, the corresponding Kernel is completely mixed with noise and no information can be extracted. In practice, this defines a measurability criterion (see [START_REF] Chatterjee | Non-linear parameter estimation with Volterra series using the method of recursive iteration through harmonic probing[END_REF]) that can be used to estimate the optimal value of N for a given amplitude of the input signal.

Influence of the assumed order of non-linearity

A study of the influence of the order of non-linearity N on estimation has been conducted on the simulated system. The simulated system of order 4 has been identified using the parameters of Table 2 but with different assumed orders of nonlinearity N ranging from 2 to 6. A white Gaussian noise with a SNR I of 80 dB (SNR B of 57.2 dB) has also been added to the output s(t) of the simulated system. In Fig. 7, the relative estimation errors made on the Kernels estimated with the different orders of non-linearity are shown. The case N = 4 will be the reference as it is the exact order of non-linearity of the simulated system.

If this order of non-linearity is underestimated (cases N = 2 and N = 3), the method gives inaccurate results. This is because of the link which exists between the different estimated Kernels {h n (t)} n∈{1...N } and the extracted impulse responses {g n (t)} n∈{1...N } . This link is the matrix A T c , as seen in Eq. ( 18). The first coefficients of the matrix A T c are given in Eq. ( 29) in the appendix. By viewing these coefficients, it's obvious that odd Kernels depend only on odd extracted impulse responses and that the same stands for even Kernels. If an impulse response of order n odd (or even) is not taken into consideration, it will have consequences on all the Kernels of order i odd (or even) lower than n. When N = 2, the non-linearities of order 3 and 4 are not taken into consideration and induce estimation errors on the Kernels of order 1 and 2. And when N = 3, the non-linearity of order 4 is not taken into consideration and induces estimation errors on the Kernel of order 2 only.

On the other hand, if the order of non-linearity N is over-estimated (case N = 6), some portion of noise will be interpreted as extracted impulse responses. As a consequence, estimations of the Kernels are slightly less precise than in the reference case. However, as can be seen in Fig. 7, the loss in accuracy is acceptable.

Thus, to ensure an estimation which is as close as possible to reality, it is better to choose the order of non-linearity N as large as possible. The upper limit of N is reached when it is impossible to extract the corresponding impulse response from the background noise. As has been shown in section 5.3, this case occurs when the SNR A corresponding to the N th Kernel reaches a certain level (≃ 0 dB for the chosen system and parameters). In practice, experimental SNR A can be calculated and used to determine the upper limit of N , as in Sec. 6.2.

Modelling acoustical transducers with cascade of Hammerstein models

In this section, acoustical transducers are represented by cascade of Hammerstein models and their Kernels are identified using the method presented in section 3.

Experimental setup

Experiments have been conducted on two acoustical transducers to identify their Kernels in a cascade of Hammerstein models representation. A standard electrodynamic loudspeaker and a prototype panel-type transducer have been tested. The panel-type transducer consists of a light and stiff plate of sandwich material (40 cm × 60 cm) on the back of which an exciter has been glued. The plate is suspended by elastics at the middle of two of its side to approximate free boundary conditions. All measurements have been made on axis at one meter from the motor of both transducers with a microphone. Measurements have been done in a semi-anechoic room.

Measured cascade of Hammerstein models Kernels

The Kernels corresponding to both systems have been measured using the previously described experimental setup. As the cascade of Hammerstein models is a non-linear model, its Kernels should be independent of the amplitude of measurement. To assess this, measurements of the Kernels corresponding to both transducers have been done using parameters given in Table 4 for 10 different amplitudes. Amplitudes were ranging from 74 to 94 dB in pressure at 1 kHz for the electro-dynamic loudspeaker. This corresponds to normal and relatively high listening levels. For the panel-based transducer, amplitude where higher, ranging from 90 dB to 110 dB in pressure at 1 kHz.

Parameter

Value Normalized value

f 1 20 Hz f 1 /f s = 1.04 × 10 -4 f 2 20 kHz f 2 /f s = 1.04 × 10 -1 f s 192 kHz T 15 s T × f s = 2.88 × 10 6 samples N 5
Table 4: Parameters chosen to identify the real system.

Measurability of each Kernels, using the criterion defined in Sec. 5.3, has been studied for the different amplitudes at which Kernels have been identified. As it is not possible to experimentally remove noise from measurements, an experimental Kernel to noise ratio (SNR XP ) is defined which corresponds to the ratio between the RMS level of the Kernel with noise and the RMS level of noise in absence of Kernel, after the complete identification procedure. This SNR XP can be computed for each Kernel separately and can be interpreted as SNR A (see Sec. 5.3).

For both transducers, SNR XP corresponding to each Kernel are presented in Fig. 8 as a function of the measurement amplitude. As expected, the SNR XP for N = 1, i.e. for the linear transfer function is linear with the measurement amplitude. One can notice that it is not the case for Kernels of order n ≥ 2. From Fig. 8, it can also be seen that for the lower amplitudes, Kernels of order 3, 4 and 5 have a SNR XP close to 0 dB. They are thus not measurable and pollute slightly the other Kernels. As has been seen in Sec. 5.4 a limited pollution is acceptable. Moreover, these Kernels become measurable as the measurement amplitude increases. On the other hand, Kernels of order n ≥ 6 are always hardly measurable for both transducers in the chosen range of amplitude with the present experimental setup. Consequently, the choice of N = 5 in the identification procedures (see Tab. 4) seems to be a reasonable compromise between pollution of the identified Kernels by noise and incomplete modelling of the system. Mean measured Kernels and their corresponding standard deviation across measurement amplitude are then given in Fig. 9 for the electro-dynamical loudspeaker and in Fig. 10 for the panel-based loudspeaker. Due to their different underlying physical principles, the linear and non-linear responses of the two transducers are quite different. The panel-based transducer has a modal behaviour and as a consequence exhibits a linear response with more dips than the electrodynamic loudspeaker. The amplitude of the different Kernels of order n ≥ 2 decreases with frequency for the electrodynamic loudspeaker, which is consistent with the physical analysis of Ref. [START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF]. For the panel-based loudspeaker, the amplitude stays globally constant with frequency. The major non-linear phenomena involved in these two transducers do not have the same physical origins.

The variability of the measured Kernels with the excitation amplitude is studied afterwards. The linear part, which is by definition independent of amplitude, exhibits no variations among the different measurements for both transducers. The non-linearities have thus been removed successfully from the linear part using the proposed method. The identified Kernels of order n ≥ 2 depends slightly on the amplitude at which they have been measured. As a consequence, the assumption that these two transducers can be represented by a cascade of Hammerstein models is a correct approximation in the chosen range of amplitude.

Prediction of the harmonic distortion generated by transducers

In this section, the previously identified cascade of Hammerstein models will be used to predict the harmonic distortion generated by both transducers.

Link between HD n , THD and cascade of Hammerstein models

To characterize distortion generated by an acoustic transducer, the following approach is classically adopted. The input of the system is assumed to be sinusoidal and non-linearities generate harmonic components at frequencies higher than the input fundamental frequency. The amplitudes of these harmonics compared to the amplitude of the fundamental are considered as representative of the non-linearity of the transducer. Total harmonic distortion (THD) and harmonic distortion of order n (HD n ) are common tools to quantify this [START_REF] Czerwinski | Multitone testing of sound system components some results and conclusions -part 2: Modeling and application[END_REF]. The THD is the square root of the ratio of the power contained in the harmonics to the power contained in the fundamental. The HD n is the same but for the n th harmonic only.

For a sinusoidal input signal x(t) = X cos(2πf t) which enters a cascade of Hammerstein models identified at the amplitude X 0 , the output signal z(t) can be written as in Eqs. (20) by using Eq. ( 2) and [START_REF] Stan | Comparison of different impulse response measurement techniques[END_REF].

z(t) = N n=1 |Γ n (X, nf )| cos [2πnf t + ∠(Γ n (X, nf ))] (20a) with Γ n (X, f ) = N k=1 X X 0 HD n (X, f ) = Γ n (X, nf ) Γ Tot (X, f ) (21a) THD(X, f ) = N n=2 Γ n (X, nf ) Γ Tot (X, f ) 2 (21b) with Γ Tot (X, f ) = N n=1 [Γ n (X, nf )] 2 (21c)
The knowledge of the Kernels in the frequency range [f 1 , f 2 ] allows for the direct computation of the THD and HD n using Eqs. [START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF]. This can be done for any value of input amplitude X and for any frequency

f in [f 1 , f 2 ].

Prediction of HD n and THD at a given amplitude

Using the different sets of Kernels measured in the previous section, the HD n and THD for the two transducers has been predicted using Eqs. [START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF]. To compare with predictions, traditional measurements using pure tones have been done using the experimental protocol depicted in section 6.1. HD n and THD have been measured this way for 50 frequencies between 50 Hz and 12 kHz. This has been done for 10 different amplitudes ranging from 74 to 94 dB in pressure for the electro-dynamic loudspeaker and from 90 to 110 dB for the panel-based one.

In Figs. 11 and 12 the predictions for the total harmonic distortion (THD), and for the harmonic distortion of order 2 and 3 (HD 2 and HD 3 ) made using equations ( 21) are shown for the electrodynamic loudspeaker and the panel-based loudspeaker. It can be seen that the agreement between measured and predicted data is satisfying over the entire frequency range for the electrodynamic loudspeaker (Fig. 11). For the panel-based loudspeakers, the agreement is also good (Fig. 12). Below 200 Hz the predictions sometime underestimate HD 2 , HD 3 . However, evaluation of harmonic distortion using series does not necessary converge to the desired result [START_REF] Chatterjee | Non-linear parameter estimation with Volterra series using the method of recursive iteration through harmonic probing[END_REF][START_REF] Chatterjee | Convergence analysis of Volterra series response of nonlinear systems subjected to harmonic excitation[END_REF]. Indeed, the number of terms N to be used in the series of Eq. ( 21) has to be carefully chosen for predictions to converge to measurements. To study that point, THD has been predicted using 2, 3, 4 or 5 terms in Eqs.

Prediction of HD n and THD for different amplitudes

To have an overview of the quality of the predictions depending on the amplitude at which Kernels have been measured (X m ) and on the amplitude at which prediction are made (X p ), a mean error has been introduced. This error is defined in Eq. 22. This error has been computed for the following frequency bands: [45, 355] Hz (octave bands 63 Hz to 250 Hz), [355, 2800] Hz (octave bands 500 Hz to 2 kHz) and [2800, 11200] Hz (octave bands 4 kHz to 8 kHz). The error in dB in each of these frequency bands is shown for the two transducers in Figs. 14 and15.

∆ [fA,fB ] (X m , X p ) = ... 20 log 10 mean [fA,fB ] THD meas (X p , f ) -THD pred (X m , X p , f ) THD meas (X p , f ) (22)
Figs. 14 (a), (b), and (c) give the resulting error for the electrodynamic loudspeaker. In Fig. 14 (a), the error in the frequency band [45,355] Hz is shown. It can be seen that this error is acceptable, around -6 dB. The minimum of -8 dB is reached when the THD is predicted for low values of X p . As the amplitude of prediction X p increases, the error increases too in this frequency band. In Figs. 14 (

Discussion

The presented method has been applied to two real acoustic transducers, a classic electro-dynamical loudspeaker and a panel-based transducer. As a complete model describing the non-linear behaviour of these transducers is accessible using the previous method (see section 6), an analysis of the two systems can be performed from this point of view. The first difference between them appears when viewing their linear response. The panel-based transducer exhibits a well known diffuse behaviour in the high frequency range, whereas the electrodynamic loudspeaker does not. This diffuse behaviour is also present in the Kernels of higher order. For the two transducers, the amplitudes of the higher order Kernels can also be compared. Kernels of the panel-based loudspeaker generally have a lower amplitude than Kernels of the electro-dynamic loudspeaker. For the electro-dynamic loudspeaker, the amplitude of Kernels of order n ≥ 2 have a tendency to decrease with the frequency, which is coherent with Refs. [START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a Volterra series expansion[END_REF][START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF], whereas for the panel-based loudspeaker there is no global variation with frequency. Panel-based transducers thus seems to generate less distortion than electrodynamic loudspeakers, and such distortion is almost constant with frequency.

In section 7, prediction of the harmonic distortion of order n (HD n ) and the total harmonic distortion (THD) from the identified Kernels have been performed. The originality of the present approach is that an analytical formula involving the identified Kernels and the amplitude at which they have been identified allows one to predict HD n and THD at different amplitudes. This is an advancement compared to current literature where HD n and THD are usually predicted only for a given amplitude [START_REF] Muller | Transfer-function measurement with sweeps[END_REF][START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a Volterra series expansion[END_REF]. The results obtained for the two systems when comparing this approach to the traditional one at different amplitudes for HD n and THD are satisfying. These results validate that cascade of Hammerstein models are a well adapted model.

Conclusion

In this paper a simple and rapid model based procedure to measure non-linearities of a vibrating system has been presented mathematically, validated by simulation and finally applied to two acoustical transducers. Cascade of Hammerstein models have been chosen here to model the non-linearities. It has been shown on simulations that the identification method is very accurate. Harmonic distortion generated by these devices is afterwards precisely predicted using this model. This method, coming from the audio community, can be of great interest in the more general field of vibrations. In modal analysis, for example, a common limitation to access high frequencies is the signal to noise ratio (SNR). Continuous acoustical excitations with high levels are thus of great interest to increase the SNR, but only if the resulting signals are not polluted by non-linearities. Using the presented method, structures can be acoustically excited at high amplitude levels with the assurance that the non-linear part of the excitation present in the measurements can be completely removed. This method can thus help to solve practical problems which are commonly encountered in experimentations involving vibrations.

This method can also be interesting for transducer quality assessment. It is now known that traditional non-linear measurements tools (such as HD n and THD) correlate poorly with subjective experiences [START_REF] Voishvillo | Assessment of nonlinearity in transducers and sound systems -from THD to perceptual models[END_REF]. However, the present approach gives a fine, input-independent representation of the linear and non-linear characteristics of real transducers. As a consequence, simulations of the non-linear responses of identified or calculated transducers can be easily performed. This can be useful when searching for new criteria to assess the decrease of quality caused by non-linearities in acoustical transducers.

Subsequently, it can easily be shown that they satisfy the recurrence relation given in Eq. [START_REF] Baumgartner | Complete identification of a class of nonlinear-systems from steady-state frequencyresponse[END_REF].

k = 0 T 0 (x) = 1 (24a) k = 1 T 1 (x) = x (24b) k > 1 T k+1 (x) = 2xT k (x) -T k-1 (x) (24c)
Then, by writing the polynomials as in Eq. 25 , one can obtain Eq. ( 26), using Eq. ( 24), and find the coefficients of the matrix A.

T k (x) = k i=0 A(i, k)x i (25)
i = 0 A(0, k + 1) = -A(0, k -1) (26a)

0 < i < k A(i, k + 1) = 2A(i -1, k) -A(i, k -1) (26b) i k A(i, k + 1) = 2A(i -1, k) (26c)
The linearisation of the polynomials can now be rewritten in a matrix form, as in Eq. [START_REF] Boyd | Fading memory and the problem of approximating nonlinear operators with Volterra series[END_REF].

    1 cos(x) ... cos(N x)     = A     1 cos(x) ... cos N (x)     (27)
By inverting Eq. ( 27), Eq. ( 7) is directly obtained and this results in Eq. ( 28) which gives explicitly the C matrix.

C = A -1 (28)
The matrix A T c , necessary to access to {h n (t)} n∈[1,N ] , is the matrix A without the first column and the first row, as seen in Eq. [START_REF] Colloms | Distortion mechanisms of distributed mode (DM) panel loudspeakers[END_REF]. To avoid the implementation of the recurrence, the transpose of the A c matrix of order 8, which is sufficient for practical use, is given in Eq. [START_REF] Bicken | Frequency domain Kernel estimation for 2nd-order Volterra models using random multi-tone excitation[END_REF].

A T c =             1

Figure 1 :

 1 Figure 1: Block diagram representation of a cascade of Hammerstein models.

4 Figure 2 :

 42 Figure 2: Separation of the different orders of non-linearity after convolution with y(t).

Figure 3 :

 3 Figure 3: Overview of the method used to identify a cascade of Hammerstein models.

  and H 4 (f ) are shown in Fig 4. The estimated Kernels are very close to the original ones over almost the entire frequency range [f 1 , f 2 ].For the frequency regions close to f 1 and f 2 , the estimated Kernels deviates slightly from the original ones, especially the highest orders. This illustrates the limits of the band-limited inverse filter defined by Eq. (15).

Figure 4 :

 4 Figure 4: Magnitude (a) and phase (b) of the frequency responses of each Kernel of the simulated system. Originals are shown in solid or dotted lines and estimations with o,+,* or x.

Figure 5 :

 5 Figure 5: Relative errors made during the estimation of the different Kernels. (a) : H 1 (f ), (b) : H 2 (f ), (c) : H 3 (f ), (d) : H 4 (f ).

Figure 6 :

 6 Figure 6: Relative errors made during the estimation of the different Kernels for different noise levels. (a) : H 1 (f ), (b) : H 2 (f ), (c) : H 3 (f ), (d) : H 4 (f ).

Figure 7 :

 7 Figure 7: Relative errors made during the estimation of the different Kernels for different chosen orders of non-linearity N . (a) : H 1 (f ) estimated with N = {2, 3, 4, 6}. (b) : H 2 (f ) estimated with N = {2, 3, 4, 6}. (c) : H 3 (f ) estimated with N = {3, 4, 6}. (d) : H 4 (f ) estimated with N = {4, 6}.

Figure 8 :

 8 Figure 8: Experimental Kernel to noise ratio after the complete identification procedure (SNR XP ) as a function of the amplitude of measurement for each Kernel. (a) Electro-dynamic loudspeaker. (b) Panel-based transducer.

Figure 9 :

 9 Figure 9: Mean measured Kernels of the cascade of Hammerstein models for the electrodynamic loudspeaker (solid line) and the corresponding standard deviations (dashed lines). (a) : H 1 (f ), (b) : H 2 (f ), (c) : H 3 (f ), (d) : H 4 (f ).

Figure 10 :

 10 Figure 10: Mean measured Kernels of the cascade of Hammerstein models for the panel-based loudspeaker (solid line) and the corresponding standard deviations (dashed lines). (a) : H 1 (f ), (b) : H 2 (f ), (c) : H 3 (f ), (d) : H 4 (f ).

Figure 11 :

 11 Figure 11: Comparisons between measurements (circles) and predictions (solid lines) at 85 dB for the HD 2 , HD 3 and THD of the electrodynamic loudspeaker. Kernels identified at 86 dB have been used for predictions with N = 5. Harmonic distortion of order 2, HD 2 . (b) Harmonic distortion of order 3, HD 3 . (c) Total harmonic distortion, THD.

  [START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF] and compared to measurements. The mean error in frequency between predictions and measurements for different values of N is presented in Fig.13for both transducers. It has been computed for the following frequency bands: [45, 180] Hz (octave bands 63 Hz and 125 Hz), [180, 710] Hz (octave bands 250 Hz and 500 Hz), [710, 2800] Hz (octave bands 1 kHz and 2 kHz) and [2800, 11200] Hz (octave bands 4 kHz and 8 kHz). For the electro-dynamic loudspeaker in the two upper frequency bands, as N increases the mean error become lower. For the lower frequency bands, the mean error increases with N until N = 4 where it starts decreasing. For the panel-based loudspeaker, mean error globally decreases with N . Thus, the choice of N = 5 (see Tab. 4) in the identification procedure and for the predictions leads to globally convergent results in the chosen ranges of amplitude and frequency.

Figure 12 :

 12 Figure 12: Comparisons between measurements (circles) and predictions (solid lines) at 101 dB for the HD 2 , HD 3 and THD of the panel-based transducer. Kernels identified at 105.5 dB have been used for predictions with N = 5. (a) Harmonic distortion of order 2, HD 2 . (b) Harmonic distortion of order 3, HD 3 . (c) Total harmonic distortion, THD.

Figure 13 :

 13 Figure 13: Convergence of the prediction error averaged over different frequency bands for the electro-dynamic and panelbased transducers. (a) Electro-dynamic loudspeaker: Kernels identified at 88 dB have been used for predictions at 85.1 dB. (b) Panel-based loudspeaker: Kernels identified at 99.5 dB have been used for predictions at 92.2 dB.

  Figs.14 (a), (b), and (c) give the resulting error for the electrodynamic loudspeaker. In Fig.14 (a), the error in the frequency band[45, 355] Hz is shown. It can be seen that this error is acceptable, around -6 dB. The minimum of -8 dB is reached when the THD is predicted for low values of X p . As the amplitude of prediction X p increases, the error increases too in this frequency band. In Figs.[START_REF] Griesinger | Beyond MLS -Occupied hall measurement with FFT techniques[END_REF] (b) and (c), errors for frequency bands [355, 2800] and [2800, 11200] Hz are shown. Error values in these frequency bands are significantly lower than in the previous one. The minimums of these errors, which are -12 dB and -10 dB, can be seen around the diagonals. Predictions are then precise in these frequency bands. Figs. 15 (a), (b), and (c) give the same errors for the panel-based transducer. In Figs. 15 (a) and (c), the errors in frequency bands [45, 355] Hz and [2800, 11200] Hz are shown. These results are acceptable though less satisfying than the results obtained for the electrodynamic loudspeaker. In Fig. 15 (b), errors for frequency band [355, 2800] are shown. Error values in these frequency bands for the panel-type loudspeakers are of the same magnitude as for the electrodynamic loudspeaker, remaining low, around -10 dB.

  ∀k ∈ N, cos(kφ) = T k [cos(φ)]

Figure 14 :

 14 Figure 14: Mean error in the frequency band [45, 355] Hz (a), [355, 2800] Hz (b) and [2800, 11200] Hz (c) for the electrodynamic loudspeaker. Amplitude of measurements of the Kernels are given on the x-axis.

Figure 15 :

 15 Figure 15: Mean error in the frequency band [45, 355] (a), [355, 2800] (b) and [2800, 11200] Hz (c) for the panel-based transducer. Amplitude of measurements of the Kernels are given on the x-axis.

Table 1 :

 1 f zeros (kHz) f zeros /f s |p zeros | f poles (kHz) f poles /f s Parameters used for the simulation of the cascade of Hammerstein models of order N = 4.

	|p poles | Gains

k C(k, n)H k (f )(20b)THD and HD n can thus be directly identified from Eq. (20) and expressed as :

Acknowledgements

The authors would like to thank Thomas Hélie for his helpful advice on this work.

Appendix: Computation of the matrix C

Chebyshev polynomials {T k [cos(φ)]} k∈N are defined by Eq. ( 23).