Antonin Novak  homepage
Postdoctoral Research Fellow
Laboratoire d'Acoustique de l'Université du Maine 
Index Terms: Generalized Hammerstein model, Higher Harmonic Frequency Response, identification, nonlinear system, nonlinear systems in series.
ALMOST all realworld systems exhibit nonlinear behavior to varying degrees and must be described using a nonlinear model.
For such nonlinear systems, the Volterra series representation is usually considered as an effective model and its use in nonlinear system identification and analysis has become widespread from the early 80’s [1], [2], [3]. Nevertheless, the Volterra model lays down
the calculation of multidimensional kernels and most applications are limited to the second or third order, because of analytical
difficulties and computational cost. Simplified Volterrabased models, namely Hammerstein model (NL), Wiener model (LN),
HammersteinWiener cascade (NLN) or WienerHammerstein cascade (LNL) are then often preferred in the case of openloop
systems, because of their simpler structure and lower computational cost. Furthermore, for a better accuracy of the estimation, these
simple models can be extended to socalled generalized models, such as the generalized Hammerstein model, as illustrated in Fig. 1, in
the case of a SingleInput SingleOutput (SISO) system. This generalized Hammerstein model is made up of N parallel
branches, with each branch consisting of a linear filter G_{n}(f) preceded by a nth order power static nonlinear function, for
n = 1,N.
The identification of a nonlinear system using a generalized Hammerstein model consists in estimating the unknown linear filters
G_{n}(f) from the known input and output signals x(t) and y(t) respectively. Most of identification methods of Hammerstein or Wiener
systems are based on parametric models [4], which means that each unknown linear filter G_{n}(f) is modeled by a
parametric frequency response function (FRF) model such as ARMA structure. However, a nonparametric method
based on exponential sweptsine input signal x(t) and nonlinear convolution has been proposed in [5] and [6] for
the estimation of the modulus G_{n}(f), and recently extended to the estimation of both the modulus and the phase
of G_{n}(f) through the synchronization of the excitation swept sine signal [7], [8]. This method, called Synchronized
SweptSine Method in this paper, has been successfully applied to the study of various nonlinear systems [9], [10],
[8].
However, for experimental situations, where a physical source (sound, light, vibration...) is driven by a sweptsine input signal for
exciting the system under test (SUT), the unwanted nonlinearities of the source will be mixed with the nonlinearities of the SUT to be
estimated. Then it might be of great interest to model the whole system (source + SUT) as two nonlinear dynamic systems in series, in
order to be able to estimate independently each physical subsystem.
The present work consequently aims at developing a method allowing the identification of the nonlinear SUT without taking into
account the effects of the nonlinearities of the first subsystem. Recently, Ege et al. [11] have proposed such a method for the estimation
of the nonlinearities of a piano soundboard, when the soundboard is acoustically excited with a loudspeaker. Since both the source and
the soundboard are nonlinear, the characterization of the piano soundboard requires measurements of both the distorted pressure
generated by the loudspeaker and the local acceleration of the soundboard. The method presented by Ege et al. then allows the
estimation of the distortion rate of the soundboard only, but does not allow to identify the different contributions of the nonlinearities of
the SUT.
In this paper, an offline method allowing the identification of the second subsystem of two nonlinear systems in
series is presented. This method is based on the estimation of the Higher Harmonic Frequency Responses (HHFRs)
[9]. In Section II, the principle of the method is detailed. In Section III, the method is experimentally validated in
the wellknown framework of nonlinear propagation of acoustic waves. Finally, conclusions are presented in section
IV.
Consider two SISO stable, timeinvariant nonlinear systems (NLS) connected in series in an open loop, as depicted in
Fig. 2.
The signal x(t) is the input of the first NLS. The signal u(t) is the output of the first NLS and the input of the second NLS. The signal y(t) is the output of the second NLS. Consequently, y(t) may be seen as the output of the whole system for the input signal x(t). Lastly, it is supposed that these three signals are available for the analysis.
The method presented here then allows the identification of the second NLS. Adding to the assumption that both subsystems are in
series without any feedback (openloop systems), as shown in Fig. 2, we suppose in the following that the first NLS belongs to the
periodinsameperiodout (PISPO) class of systems [4], and that the second NLS may be described by a Nthorder generalized
Hammerstein model, as shown in Fig. 3.
The identification of the system is then equivalent to estimating the linear filters G_{n}(f), n = 1,N, from the measured signals u(t) and y(t), the input signal x(t) being known. The identification process is based on the offline estimation of both the HHFRs _{l}^{(u,x)}(f) between x(t) and u(t), and the HHFRs _{l}^{(y,x)}(f) between x(t) and y(t), for l = 1,L.
We recall that, given an input signal a(t) and an output signal b(t) of a NLS, the HHFR _{l}^{(b,a)}(f) may be seen as the contribution, in both amplitude and phase, of the lth harmonic at the output, for a sine at frequency f at the input, as

(1) 
Therefore, the HHFRs can be interpreted as frequency responses of fundamental and higher harmonics [9]. Several methods have been developed to estimate the HHFRs. The most intuitive but time consuming one is based on an harmonic excitation of the SUT and the process is repeated by changing the input frequency step by step. In [12], two methods have been proposed, the first one being based on FFT techniques to estimate the autospectrum and phase information of HHFRs and the second one using IQ demodulation. In this paper, the socalled Synchronized SweptSine Method [9] is used to estimate the HHFRs, but any other technique allowing the estimation of _{l}^{(y,x)}(f) and _{l}^{(u,x)}(f) in both amplitude and phase can be used. The method we propose also involves estimating the HHFRs _{l}^{(un,x) }(f) of powers of the signal u(t), for n = 1,N, corresponding to each polynomial input of the generalized Hammerstein model describing the second NLS. The signal u(t), already distorted by the first nonlinear subsystem NL_{1}, is then taken to the power of n in order to estimate the HHFRs _{l}^{(un,x) }(f).
The HHFRs _{l}^{(y,x)}(f) of the output signal y(t) result in the combination of all HHFRs _{l}^{(un,x) }(f) after filtering by filters G_{n}(f). As detailed in Appendix A, the relation between the HHFRs _{l}^{(un,x) }(f), _{l}^{(y,x)}(f) and the linear filters G_{n}(f) can indeed be written in a matrix form as

(2) 
Equation (2) can then be solved for unknown G_{n}(f) by using a square matrix inversion in the case L = N, or by using a
pseudoinversion in the case L > N. The number of harmonics L chosen for the estimation of HHFRs must be equal or greater than
the number of the branches N of the generalized Hammerstein model. The matrix (pseudo)inversion must be computed for each
frequency f separately. Nevertheless, for frequencies corresponding to weak values of the output u(t) of the first nonlinear system, the
matrix L × N may be illconditioned. This can lead to meaningless solutions. Consequently, the values of _{l}^{(u,x)}(f) have to be
calculated first in order to achieve a good matrix conditioning. Appendix A details the method when exciting the two nonlinear systems
in series by a sine wave of frequency f_{0}.
III. Validation on an Acoustic Waveguide
A. Experimental set up description
To illustrate and validate the method described in the previous section, a particular case of two nonlinear systems
connected in series is considered here. It consists of a compression driver exciting an acoustic waveguide as shown in Fig.
4.
A high level excitation voltage is supplied to the driver so that the behavior of the driver is supposed to be nonlinear. The driver generates a high level acoustic signal at the input of the waveguide so that a nonlinear acoustic propagation occurs in the waveguide. Nonlinear acoustic propagation due to a high level source can lead to shock wave for long distance propagation [13]. Indeed, local small nonlinear perturbations of the acoustic wave due to high sound pressure level are cumulative along the propagation and may distort the waveform considerably for sufficiently large sourcetoreceiver distances. Nevertheless, interest is focused here on a shorter distance, leading to the observation of weak nonlinear distortion, mathematically described by Burgers’ equation [14], as summarized in Appendix B.
In addition to this weak nonlinear distortion, the acoustic pressure generated at the input of the waveguide is distorted, due to the
nonlinear behavior of the driver, in such a way that the acoustic pressure propagating in the waveguide contains additional
loudspeakerdriven nonlinear components. Identification of the nonlinearities solely due to nonlinear propagation in the waveguide is
then carried out, using the identification method presented in the previous section. The experimental setup is depicted in
Fig.5.
The setup is made up of a cylindrical airfilled tube (6 m long and 58 mm internal diameter) coupled to a compression driver (JBL model 2446H) at one of its ends. The other end of the tube is loaded by an absorbing termination in order to avoid sound reflections. The input signal x(t) is a swept sine satisfying the conditions required for the Synchronized SweptSine Method, as detailed in [7]. Two microphones M1 and M2 (acceleration compensated piezoelectrical gauges PCB M116B) are flush mounted on the pipe wall. The first one is set at 20 cm from the source, whereas the second one is set at the distance of 4.2 m from the source. In this experiment, the first NLS consists of the power amplifier, the compression driver and the part of the guide between the driver and the first microphone. The second NLS is the part of the guide between the two microphones.
Pressure level provided by the driver at the input of the tube is chosen to be high enough for exhibiting weakly nonlinear acoustic propagation. In this experiment, the RMS pressure is 650 Pa. Measurements are performed from 300 Hz to 5000 Hz and sampling frequency is set to 96 kHz. For this frequency range, the absorbing termination is efficient, so that the outgoing wave reflected from the end of the tube is negligible. Thus, the identification results should be compared with the nonlinear traveling plane wave theory, based on Burgers’ equation resolution and detailed in Appendix B. Note that the cutoff frequency, below which the plane wave approximation is valid, is around 3300 Hz.
B. Estimation of Higher Harmonic Frequency Responses
Both acoustic pressure signals u(t) and y(t) (from microphones M1 and M2) are independently analyzed thanks to the Synchronized SweptSine Method, by estimating HHFRs up to the 9thorder HHFR, with x(t) the input signal. In the following, only the first, second and third HHFRs are plotted to improve readability of the graphs. In Fig. 6, the moduli of the HHFRs _{l}^{(u,x)}(f) at the first microphone location are given. Similarly, Fig. 7 shows the moduli of the HHFRs _{l}^{(y,x)}(f) at the second microphone location.
Fig. 6 and 7 may be read as follows. For a sine wave excitation at a given frequency f, the acoustic pressure signal under analysis is made up of a fundamental frequency, a second and a thirdorder frequency components, which amplitude level is available at f, 2f and 3f frequencies, respectively.
Such representation of lthorder harmonics allows to read easily the levels of the HHFRs, relatively to a given value of the
fundamental and consequently to evaluate the distortion, for each input frequency.
In Fig. 6, fundamental level of first microphone acoustic pressure signal, _{1}^{(u,x)}(f), is quite constant, around 140150 dBSPL
over the whole frequency span. The level of second harmonic _{2}^{(u,x)}(f) is lower than the fundamental one and
increases slightly with frequency. _{3}^{(u,x)}(f) level is lower than _{2}^{(u,x)}(f) except below 800 Hz input frequency,
leading to a minimum dynamic value between _{1}^{(u,x)}(f) and _{3}^{(u,x)}(f) of 25 dB, for an input frequency of
450 Hz.
Comparing the results at both microphones locations, the level of _{1}^{(y,x)}(f) measured by the second microphone is 1 dB lower than the level of _{1}^{(u,x)}(f) measured by the first microphone, for low frequencies. Moreover, the difference between fundamental levels at both microphones locations increases slightly with frequency above 2000 Hz. Both these effects are due to viscothermal losses, which theoretically increase with frequency [15]. Furthermore, we notice in Fig. 7 that the acoustic pressure at the second microphone location is more distorted than at the first microphone location. This is due to the effects of nonlinear propagation: _{2}^{(y,x)}(f) and _{3}^{(y,x)}(f) globally increase with frequency, as predicted by the theory [16].
C. Generalized Hammerstein Model and Burgers’ Theory
From experimental signals x(t), u(t) and y(t), HHFRs _{l}^{(un,x)
}(f) and _{l}^{(y,x)}(f) are calculated, with l = 1,9. Then, the linear
filters G_{n}(f), for n = 1,3, are estimated from (2). They characterize the behavior of the second nonlinear system,
as far as a generalized Hammerstein model correctly fits the nonlinear propagation in the waveguide. The HHFRs
_{l}^{(y,u)}(f) for a given input u(t) may then be calculated from the identified generalized Hammerstein model as depicted in
Fig. 8.
Moreover, the nonlinear acoustic propagation in the waveguide is completely described by the theory of Burgers detailed in Appendix B. A theoretical expression of HHFRs for a given input pressure is then available, noted H_{l}(f).
As a validation, we aim at comparing the calculated HHFRs _{l}^{(y,u)}(f) from the estimated model and the theoretical HHFRs H_{l}(f). Besides, we also compare in the following the waveforms of the measured, theoretical and modelbased resynthesized acoustical pressures.
Two comparisons are presented here. First, in Figs. 9 and 10, the theoretical results calculated using Burger’s theory are compared with the measured data including the whole system (NL_{1} in series with NL_{2}). This comparison represents a case for which the nonlinearities of the first subsystem (NL_{1}) are mixed with the nonlinearities of the SUT (NL_{2}).
In Fig. 9, the theoretical HHFRs H_{l}(f) are compared with the HHFRs _{l}^{(y,x)}(f) of the whole system which contains
information of nonlinearities of both systems NL_{1} and NL_{2}. Discrepancies between theoretical H_{l}(f) and measured
_{l}^{(y,x)}(f) are then obvious, especially for second and third HHFRs. In particular, _{3}^{(y,x)}(f) overestimates H_{3}(f) in
the frequency span from 300 Hz up to 800 Hz, due to the compression driver distortion contribution, as noted in
Fig. 6.
Fig. 10 shows the comparison of the time waveforms of the acoustical pressure for input signal at 1400 Hz. The theoretical
waveform calculated using the Burger’s theory (dashed line) clearly differs from the experimental signal captured by the microphone 2
(solid line), the RMS error between both waveforms reaching 54.3 Pa due to the nonlinearities of the first subsystem (NL_{1}).
The second comparison (Figs. 11 and 12) shows the theoretical results calculated using Burger’s theory compared with the estimated
model. In other words, the identification of the second nonlinear subsystem (the acoustical waveguide) is achieved while getting rid of
the effects of the first subsystem (the compression driver), according to Equation (2).
In Fig. 11, the theoretical HHFRs H_{l}(f) are compared with the HHFRs _{l}^{(y,u)}(f) obtained from the estimated model of the SUT
(NL_{2}). The results show that the identification of the propagation in the waveguide (NL_{2}) is in agreement with Burgers’ theory, contrary
to the results proposed in Fig. 9.
The same conclusion can be drawn with regard to the comparison of time waveforms of the acoustical pressure (Fig. 12). The waveform (solid line) resynthesized from the estimated model matches well with the expected theoretical waveform calculated using the Burger’s theory (dashed line), the RMS error being 9.3 Pa, almost six times lower than in the previous case.
The accuracy of the proposed identification is also confirmed by calculating the root mean squared error E_{l} for each HHFR _{l}^{(y,u)}(f) as

for which values are given, for both cases, in Table I.
In this paper, a method for the identification of two NLS in series has been presented.
The cascade of two NLS means that the second NLS under test is excited by a distorted signal, which is the output of the first nonlinear subsystem. In practice, this corresponds to the classical identification case of a system under test excited by a nonlinear device, as a loudspeaker or a shaker.
The method proposed in this paper identifies the second system using a generalized Hammerstein representation. On one hand, one of the characteristics of the method is to properly estimate the HHFRs of input and output of nonlinear SUT. Among several available methods, the Synchronized SweptSine method has been chosen for both its rapidity and robustness [7].
On the other hand, it is also worth noting that the method is low time consuming. The proposed method is indeed easy to implement and no special algorithm is required. In addition, the method has no need for any knowledge of the first subsystem, as far as the generalized Hammerstein model correctly fits the nonlinear behavior of the (second) nonlinear system under test.
The method allows in particular to regenerate an output signal corresponding to any given input signal. Thanks to this property we validate the method for a theoretically very wellknown physical system, i.e. the weakly nonlinear acoustic propagation in a waveguide.
The results show good agreement between both modelbased and theoretical system outputs.
Appendix A: Detailed Procedure
The following example details the estimation procedure of the unknown filters G_{n}(f) when exciting the cascade of nonlinear systems by a sine wave signal. In this example, the nonlinear system under test (the second nonlinear system) is represented by a third order generalized Hammerstein model. The estimation procedure is detailed for a frequency f_{0}, leading to estimation of the three unknowns G_{1}(f_{0}), G_{2}(f_{0}) and G_{3}(f_{0}).
First, we excite the cascade of nonlinear systems with a sine wave signal x(t) with frequency f_{0}, as depicted in Fig. 13. For frequency f_{0}, the relation between the input u(t) and the output y(t) of the nonlinear system under test can be written as

(3) 
This equation describes the relation between the first harmonics of signals u(t), u^{2}(t), u^{3}(t) and y(t).
Next, we change the excitation frequency of the excitation signal x(t) to f_{0}/2 and we study what happens for the same frequency f_{0} as in the previous case. This situation is illustrated in Fig. 14 and leads to an equation that describes the relation between the second harmonics of signals u(t), u^{2}(t), u^{3}(t) and y(t), as

(4) 
Finally, we change the excitation frequency of the excitation signal x(t) to f_{0}/3 and we still study what happens for frequency f_{0}. This situation is illustrated in Fig. 15 and leads to an equation that describes the relation between the third harmonics of signals u(t), u^{2}(t), u^{3}(t) and y(t), as

(5) 
The systems of three equations Eqs.(35) in three unknowns G_{1}(f_{0}), G_{2}(f_{0}) and G_{3}(f_{0}) can be solved. Generalizing the number of
unknowns (number of branches of the generalized Hammerstein model) to N and the number of harmonics to L, and repeating the
procedure for all the desired frequencies f leads to the matrix expression (2).
Appendix B: Burgers’ Equation for Traveling Plane Waves
Generalized Burgers’ equation with thermoviscous losses in the boundary layer, for a pure traveling wave is [16] [17] [18]

(6) 
Equation (6) is built from an adapted method for cumulative phenomena, the Multiple Scale Method (MSM), which has been validated up to a distance close to the critical shock wave distance [14]. This equation has no analytical known solution.
Solution is therefore obtained numerically in the frequency domain using a Fourier series decomposition [17], expressed as

(7) 
Introducing (6) into (7) leads to the following expressions

(8) 

(9) 
In the frame of this work, these Fourier series are truncated to L harmonics, l = 1,3. Then a resolution on σ by small steps Δσ
(finite difference method) knowing the initial condition at σ = 0 is performed. Spatial step is typically here 1cm. The resolution method
used is a prediction correction method. At first order, the classical Euler method is used, and correction is performed using the Adams
Moulton second order method.
The solution of Burgers’ equation was experimentally verified by Menguy et al. [17]. This theory serves in this work as a reference for testing the proposed identification method, by deducing the theoretical moduli of HHFRs, according to

(10) 
Acknowledgment
The authors would like to thank Joel Gilbert and Michal Bednarik for their help in Burgers’ theory for the numerical simulation of the HHFRs.
© 2014 Antonin Novak  homepage (All Rights Reserved)