PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique Applied to Nonlinear Audio Devices Emulation* Thomas Schmitz, AND Jean-Jacques Embrechts, AES Member (
[email protected])
(
[email protected])
Department of Electrical Engineering and Computer Science, University of Li`ege, Belgium
Numerous audio systems are nonlinear. It is thus of a great importance to study them and understand how they work. Volterra series model and its subclass (cascade Hammerstein model) are usual ways to model nonlinear systems. In this paper we first describe some difficulties that are encountered during the implementation of these identification methods. Solutions are proposed to solve each of these problems, some are originals, others are inspired by the literature. The different corrections provided in this paper lead to the Hammerstein Kernels Identifications by Sine Sweep (HKISS) method. To evaluate its efficiency, we emulate a distortion guitar effect and two tubes guitar amplifiers.
0 INTRODUCTION The modeling of nonlinear systems has been a central topic in many engineering areas, as most real-world devices exhibit nonlinear behaviors. The study of these topics is of great interest as they include a large range of applications such as industrial processes, control systems, economic data, life sciences, social systems, physics and many more. This project focuses on one particularly interesting application: the modeling of a guitar player’s audio chain. Such a chain is usually composed of compressors, distortion effects, preamplifiers, amplifiers, loudspeakers and microphones. Most of these elements exhibit nonlinear characteristics [1, 2, 3]. The idea is to replace these devices by a computer emulation that could be suitable for guitar players in many situations [1, 4]. Most weakly nonlinear systems can be approximated by linearization. In that case, the modeling can be achieved by the convolution technique [5], as for classical Linear Time Invariant (LTI) systems. However, stronger nonlinearities need more specific models in order to keep an accurate approximation of the real system. Such models include the Volterra [6], the multiple-input-single-output (MISO) [7], the single-input-single-output (SISO) [8], the multiple-input-multiple-output (MIMO) [9], the nonlinear auto-regressive moving average with exogenous input (NARMAX) [10], the neural network [11], the hybrid ge* To whom correspondence should be addressed Tel: +32-4366-
2706; e-mail:
[email protected] J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
netic algorithm [12], the extended Kalman filtering [13] and the particle filtering models [14].
x[n]
(.) 1
h1
(.) 2
h2
(.) M
hM
y[n]
Fig. 1. The Global Polynomial Hammerstein model
In our research, we focus on a technique based on the Volterra series model [6, 15], which recently received an increased interest with the improvement of computer’s performance. Considering only static (memoryless) nonlinearities, this model can be simplified into a Hammerstein model (see Fig. 1) which is a cascade of M parallel (the system’s order) branches where each branch m is composed by a static monomial nonlinearity (putting the input signal to the power m) followed by a linear filter hm [n] 1
Th.Schmitz AND J.J.Embrechts
which is the mth Hammerstein kernel (also called diagonal Volterra kernel) [7, 15, 16, 17]. In the last few years, a nonlinear convolution method [18] has been elaborated with the aim of emulating nonlinear audio systems. Farina et al. [18, 19] indeed proposed to use the properties of the Exponential Sine Sweep (ESS) signal in order to identify the Hammerstein kernels. With this signal at the input of a nonlinear system, they showed that it is possible to separate from each other the harmonic responses created by the harmonic distortions. This method has been quite successfully applied to model nonlinearities in loudspeakers [1, 20], guitar distortion generators [21, 22] (the famous Tube-Screamer is often taken as an example) and compressors [23]. Recently, Novak et al. [24] have brought important improvements to the ESS method with their use of a synchronized swept-sine (SSS) signal. This new formulation of the sine sweep signal corrects some inaccuracies that were inherent to the original implementation, in particular the problem of the phase residue that will be described in the next section. In the same paper, these authors have also used an efficient analytical formula for the deconvolution of the sweep signal. Also recently, Vanderkooy et al. [25] have proposed a different method to correct the problem of the phase residue. In this paper, we first highlight some difficulties that must be considered in the implementation of the ESS method, in order to perform an accurate identification of nonlinear (NL) systems. The nature of these problems is described and solutions (either original or inspired by the literature) are proposed for each of them. In a second part of the paper, these solutions are validated by the emulations of real nonlinear audio devices. The paper is organized as follows : the basic identification technique is recalled in section 1. This, theory is illustrated with a simple third-order model and then extended to a given order M. In section 2 , five problems occurring in the implementation of the identification method are described. We propose an algorithm for testing the reconstruction of the Hammerstein kernels and use it to illustrate the methods that can solve these problems. This leads to the definition of an Hammerstein Kernels Identification by Sine Sweep (HKISS) method. In section 3, we apply this method and prove its efficiency on three real nonlinear audio devices: two tubes amplifiers and a distortion guitar effect. The section 4 offers an overview on the limitations of the HKISS method. Finally section 5 presents a Matlab toolbox helping to the implementation of the HKISS method.
PAPERS
1.1 Volterra Series Considering time invariant nonlinear systems, the Volterra series expresses the relationship between an input signal x(t) and its corresponding output signal y(t) as an infinite sum of multiple convolutions [6, 15]: ∞
y(t) =
Z +∞
m=1 −∞
2
−∞
vm (τ1 , . . . , τm ) (1)
m
× ∏ x(t − τi ) dτ1 . . . τm i=1
Where {vm (τ1 , . . . , τm )}∀m∈N+ are the Volterra kernels characterizing the system. 1.2 Cascade of Hammerstein Models It can be shown that any continuous nonlinear system can be represented by a set of M parallel branches composed by a static non linearity Pm (.) surrounded by two linear filters [26, 27]. The cascade Hammerstein model is a subclass of this general model where we consider a static nonlinear monomial function followed by a linear filter hm (t) (the Hammerstein kernels) as shown in Fig.1. Considering a nonlinear discrete time causal system with a finite order M, we can express the input/output relationship as: M
y[n] =
∑ hm [n] ~ xm [n]
(2)
m=1
Where ~ is the symbol of the convolution operator. We can see from Eqs. (1) and (2) that the Hammerstein model corresponds to the Volterra series model if the Volterra kernels are diagonal (i.e. the Volterra kernels differ from zero only for the values τ1 = τ2 = ... = τM ). That is why the cascade of Hammerstein models is also called a diagonal Volterra model. 1.3 Hammerstein Kernels Identification with the ESS method In order to extract the Hammerstein kernels, we need to excite the nonlinear system with a suitable input signal. Such a signal should allow us to separate the different Hammerstein kernels from each other after the deconvolution by its inverse filter. The next two subsections explain which suitable input signal is chosen and how it can be used to compute the Hammerstein kernels. 1.3.1 Exponential Sine Sweep Phase Properties In the following, our development starts from the original formulation of Farina et al. [18, 19]. The equations are expressed in discrete time for a better correspondence with the practical implementation. Let ss[n] be the ESS signal:
1 BASIC IDENTIFICATION METHOD This section briefly presents the principles of the Volterra series and the Hammerstein model. In order to ensure a good understanding of the method, we illustrate it with a third-order model, before extending it to a general order M.
Z +∞
...
∑
ss[n] = A sin(φ [n])
(3)
The phase φ [n] grows exponentially such as : φ [n] = ω1
R n .(e R − 1) fs
∀n ∈ [0, N − 1]
(4)
Where: J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
r R = (N − 1).(log ω2 )−1 is the inverse frequency changω1 ing rate. r N is the length of the sweep in samples. r ω1 , ω2 are the initial and final angular frequencies respectively. r fs is the sample rate.
The goal is to obtain M equations giving the M Hammerstein kernels hm [n] in Eq. (2). The method consists in exciting the nonlinear system by the ESS given in Eq. (3) and then to deconvolve the result by its inverse filter. Including the ESS signal (3) as an input of the model from Eq. (2) gives:
An interesting property of the ESS [1] is that the mth harmonic of the ESS can be derived from the ESS itself delayed by ∆m , as:
y[n] = A sin(φ [n]) ~ h1 [n] + A2 sin2 (φ [n]) ~ h2 [n] 3
+ A sin (φ [n]) ~ h3 [n]
m.φ [n] = φ [n + ∆m ] − B(m − 1) ∀n ≤ N − 1 − ∆m (5) Where: B = ω1 .R fs ∆m = R. log(m)
Using the trigonometric formulas [34] to decompose in order-one formulation, we have: (6)
Choosing B ∈ 2kπ, with k ∈ N allows to neglect the B term in Eq. (5) as it just adds a multiple of 2π to the sine phase. However, this constraint is not mandatory as it will be seen later. Note also that the values ∆m are real and non-integer (except for m=1). This will necessitate a phase correction during the extraction of the harmonic impulse responses. This will be explained in subsection 2.4. The property of the ESS signal in Eq. (5) is interesting because the contribution of an harmonic m in the impulse response will be delayed by −∆m from the first order impulse response, after its deconvolution by the ESS inverse filter ss[n]. In the particular case of B ∈ 2kπ, we can write: N−1 ss(φ [n]) ~ ss(φ [n]) = ss[k].ss[n − k] = C.δ [n − n ]
∑
0
k=0 ss(m.φ [n]) ~ ss(φ [n]) = C.δ [n − (n0 − ∆m )]
(7) Where C is an amplitude constant. n0 is the position along the time axis of the first order (linear) impulse response which depends on the ESS length (n0 = N − 1 in [28]). The inverse ESS can be calculated in several ways [17, 24, 28, 29, 30, 31]. According to [24], we can compute an analytical version of the Fourier transform of ss[n]: r h i ω − j L ω 1−ln( ωω ) −ω1 − π4 1 SS(ω) = 2 e (8) 2πL The advantages of the analytical inverse sine sweep computation are that C=1 in the equation Eq. 7 and the bandwidth of the inverse filter can be extended until the Nyquist frequency if needed. We must note at this point that if m 6= 1, the second convolution in Eq. (7) is not an exact delayed version of the Dirac pulse with the same amplitude as in the case m = 1. This problem will be addressed in section 2.6. 1.3.2 A Simple Third Order Model In order to illustrate the origin of the problems mentioned in the introduction, we propose to establish the Hammerstein kernels equations for a simple Hammerstein model developed up to the third order. An extension of the model up to any order M is presented in subsection 1.3.3. J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
(9)
3
y[n] = A sin(φ [n]) ~ h1 [n] A2 [cos(2φ [n]) − 1] ~ h2 [n] 2 A3 + [3 sin(φ [n]) − sin(3φ [n])] ~ h3 [n] 4
−
(10)
In order to deconvolve the signal y[n] properly, we want to change the cosine term in a sine expression. This can be done using the Hilbert Transform [32],[33] : H {sin(u(t))} = sin(u(t)) ~ h¯ (t) = − cos(u(t)) Where H {.} is the Hilbert Transform and h¯ (t) = troducing (11) in (10) and using (5) gives:
(11) 1 πt .
In-
y[n] = A sin(φ [n]) ~ h1 [n] A2 [sin(φ [n + ∆2 ] − B) ~ h¯ [n] + 1] ~ h2 [n] 2 A3 + [3 sin(φ [n]) − sin(φ [n + ∆3 ] − 2B)] ~ h3 [n] 4 +
= ss[n] ~ h1 [n] " # A ss[n + ∆2 ] cos(B) ~ h¯ [n] + ~ h2 [n] 2 − ss[n + ∆2 ] sin(B) 3ss[n] 2 A + ~ h3 [n] − ss[n + ∆3 ] cos(2B) 4 − ss[n + ∆3 ] sin(2B) ~ h¯ [n] (12) The DC offset of the even powers has been neglected as it is only related to the mean value of the ESS. Moreover, in practice, most audio devices stop the DC component. In the Fourier domain, using the trigonometric exponential property, the Euler formula and the Hilbert Transform property [32], we obtain for positive frequencies: Y ( jω) = SS( jω).H1 ( jω) jω∆2 A − jB f s − SS( jω). j.e .e .H2 ( jω) 2 jw∆3 A2 − j2B f 3.SS( jω) − SS( jω).e e s .H3 ( jω) + 4 (13) 3
Th.Schmitz AND J.J.Embrechts
PAPERS
With the following Hilbert Transform property: F {H (s(t))} = − j.sign(ω)F {s(t)}
(14)
Note that since the Hammerstein kernels are real, we can use the Hermitian property in the Fourier domain [5] to find the Hammerstein kernels for negative frequencies: Hm (− jω) = Hm∗ ( jω). By deconvolving the signal Y ( jω) with the inverse filter SS( jω) and by grouping in the Gm ( jω) terms having the ω same phase delay e j fs ∆m , we obtain: 3A2 1 0 G1 ( jω) H1 ( jω) 4 − jB G2 ( jω) = 0 − jAe2 0 H2 ( jω) (15) 2 − j2B −A e G3 ( jω) H3 ( jω) 0 0 4
By inverting this relation, we finally obtain the Hammerstein kernels: 1 0 3e j2B G1 ( jω) H1 ( jω) jB H2 ( jω) = 0 2 je 0 G2 ( jω) (16) A j2B −4e G ( jω) H3 ( jω) 3 0 0 2 A
1.3.3 Extension to the Order M The generalization up to the order M can be done by means of the trigonometric power formulas of sine [34]: −1n n 2n+1 sin (t) = (−1)k Ck2n+1 sin([2n + 1 − 2k]t) n ∑ 4 k=0 ∀n ∈ N Cn2n −1n n−1 sin2n (t) = 2n−1 ∑ (−1)k Ck2n cos([2n − 2k]t) + 2n 2 2 k=0 ∀n ∈ N0 (17) By choosing B ∈ 2kπ in Eq. (5) and using it in Eq. (17), we derive power formulas of the ESS signal for the positive frequency domain defining SSn as the Fourier transform of ssn : −1n n SS2n+1 =SS n ∑ (−1)k Ck2n+1 e jω∆2n+1−2k / fs 4 k=0 =an .SS ∀n ∈ N (18) −1n n−1 SS2n = j.SS. 2n−1 ∑ (−1)k Ck2n e jω∆2n−2k / f s 2 k=0 =bn .SS ∀n ∈ N0 The term glected. 4
Cn2n 22n
is the DC component and has been ne-
1.3.4 Deconvolution Process Following the same steps as in section 1.3.2 but applying (18) for sine powers, the deconvolution of the Hammerstein model (Eq. (2)) in the frequency domain leads to: jωn0 / fs Z( jω) = Y ( jω)SS( jω)e MM
= H1 ( jω) +
∑ [am H2m+1 ( jω) + bm H2m ( jω)]
m=1
(19) Where MM = (M − 1)/2. In the time domain, Eqs. (18,19) show that z[n] is a superposition of several versions of each Hammerstein kernel hm [n] delayed by ∆k (e jω∆k / fs in the Fourier domain). We name gm [n] the superposition of the hm [n] having the same delay ∆m . Furthermore, we consider that all these harmonic impulse responses ’start’ at n=0, which necessitates a delay of ∆m applied to all of them: G1 ( jω) H1 ( jω) .. .. (20) = T . . GM ( jω)
HM ( jω)
We can develop the elements Tu,v with u, v ≤ M using Eq. (19) in the matrix form and grouping in Gm ( jω) the Hm ( jω) having the same phase delay. 1−u v−u u+v (−1) 2 Cv 2 )=0 ∀v ≥ u & mod( v−1 Tu,v = 2 2 0 else (21) Finally the Hammerstein kernels in the positive frequency domain are given by the following relationship: H1 ( jω) G1 ( jω) .. .. −1 (22) = (T) . . HM ( jω)
GM ( jω)
1.4 Hammerstein Model And Input Level Dependency The cascade of Hammerstein models takes the amplitude of the input signal into account, for example, Fig. 2 compares the result of a power series when the input signal is a sinus of amplitude 0.5 at 1000Hz with its emulation using the Hammerstein kernels calculated with an ESS of amplitude 1. The signal is reconstructed correctly. However, for a real device, the model of the cascade of Hammerstein could no longer match perfectly and therefore the Hammerstein kernels can be slightly different, according to the amplitude of the input signal [22, 35]. For example, Fig. 3 presents the comparison between an ESS (focus around the frequency f=1000Hz) with α = 0.5 passing through the distortion part of the amplifier TSA15h [36] and its emulation using Hammerstein kernels calculated with an ESS of amplitude α = 1 ( Hm1 ( jω) ). We can see that the signal coming from the distortion effect (y1 ) is different from the emulated signal (y2 ). From now on Hmα ( jω) are the Hammerstein kernels obtained using an ESS with an amplitude α : ss[n] = α sin(φ [n])
(23)
J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
For the rest of the paper, the symbol α over the Hammerstein kernels will be omitted if α = 1. 1 y1
0.8
y2
Amplitude
0.6 0.4 0.2
the Nyquist band in such a way that the created model is equivalent to a power series. Using the ESS as input signal, the output y1 is the result of the power series (see section 2.7 for aliasing consideration). Deconvolving y1 with the ESS inverse filter ss[n] allows to calculate the Hammerstein kernels of this nonlinear system. Then, using Eq. (2) we can compute a signal y2 based on the Hammerstein kernels and the ESS powers. If the Hammerstein kernels are correctly reconstructed, y1 and y2 have to be identical.
0 ss[n]
(.) M
−0.2 ssm [n]
−0.4 0
20
40 60 Sample [n]
80
y1
100
ss[n]
Fig. 2. Comparison between a sine of amplitude 0.5 passing through a series of power (y1 ) and its emulation (y2 ) with Hammerstein kernels measured for an ESS of amplitude 1.
ssm [n] z[n]
hm [n]
Comput e Volt erra kernels
1
y1
Nonlinear Convolut ion
y2
Fig. 4. Test algorithm: y1 and y2 have to be identical
y2
Amplitude
0.5
0
−0.5
−1 0
50
100 Sample [n]
150
200
Fig. 3. Comparison between an ESS of amplitude 0.5 passing through the distortion effect TSA15h (y1 ) and its simulation (y2 ) by means of Hammerstein kernels measured for an ESS of amplitude 1.
2.2 The Best Case Following our top down approach, we first present our best match between the power series y1 (up to the order M=6) and its emulation by means of the Hammerstein kernels y2 (see Fig. 5). For clarity, we choose to zoom on a part of the sweep signals, (arbitrarily) around the frequency f = 1000Hz. The power series is correctly reconstructed, as it can be seen in Fig. 5. This best match has been obtained by following all the precautions and corrections described in the next subsections. 1.8
y1
1.6
y2
1.4
2 OPTIMIZATION OF THE HAMMERSTEIN KERNELS IDENTIFICATION Eq. (22) can be used to calculate the Hammerstein kernels. However several difficulties must be taken into account to make an accurate emulation of a nonlinear system. In order to illustrate them, we choose a top down approach, starting with the optimized case (subsection 2.2) where all the precautions are taken and then removing them one at a time to prove their necessities (subsections 2.3, 2.4, 2.5, 2.6). Note that other problems such as artifacts in the impulse responses have been identified and analyzed in [37, 38, 39]. 2.1 Test Algorithm We propose a test algorithm (see Fig. 4) based on the Hammerstein model described by Eq. (2). The kernels hm [n] are chosen to be simple all-pass filters restricted to J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
Amplitude
1.2 1 0.8 0.6 0.4 0.2 0 0
50
100 Sample [n]
150
Fig. 5. Comparison between y1 [n] and y2 [n] at 1000Hz (best case)
2.3 Hammerstein Kernels Phase Errors The basic model (16) shows that the Hammerstein kernels depend on a constant factor B introduced in the phase of the ESS. During the general formulation given at Eq. (22) we have made the simplifying assumption that 5
Th.Schmitz AND J.J.Embrechts
PAPERS
y
1.2
1
y
2
1 0.8 Amplitude
B ∈ 2kπ, with k ∈ N (see section 1.3.1 ). Forcing B close enough from 2kπ would imply to tweak the ESS parameters. This constraint and other conditions discussed in the next subsection would to a difficult optimization problem that could have no exact solution, leading to compromises on the ESS parameters values.As mentioned in the introduction, Novak et al. [24] have overcome this difficulty by their use of a synchronized swept-sine signal. In this paper based on the original formulation of the sweep signal [18, 19], we rather propose to directly include a correction by multiplying term by term the (T)−1 matrix with CorrB , a corrective M by M matrix, taking into account the B factor: j0 e . . . e j(v−1)B .. CorrB (u, v) = ... . . . (24) .
0.6 0.4 0.2 0 −0.2 −0.4 −0.6 0
50
100 Sample [n]
150
200
Fig. 7. Comparison between y1 [n] and y2 [n] at 1000Hz when the B correction is not applied
e j0 . . . e j(v−1)B This correction is equivalent to the one given in [25].In this last paper, the authors have introduced a generalized form of the Hilbert Transform which results in shifting all the phases by the same amount, for all frequencies. Indeed, after application of the CorrB matrix in Eq. 24, all the phases of the Hammerstein kernels are equal to zero (see Fig.6), which is the condition imposed in [25]. Fig. 7 shows the comparison between y1 and y2 when the matrix CorrB is not applied (with mod(B, 2π) = 0.7). The emulated signal is very sensitive to the Hammerstein kernels B factor as it no longer corresponds to the power series signal. Hv1 Hv2 Hv3 Hv4 Hv5 Hv6
150 100
Phi °
50 0 −50 −100
phase mismatch between the kernels gm [n] that impacts the computation of the Hammerstein kernels (see Eq. (22)).
Ceil(∆4 )
Ceil(∆3 )
Ceil(∆2 )
g1 [n] g2 [n] g3 [n] g4 [n]
∆2 ∆3
∆4
4
[n]
3
2
N-1
Fig. 8. Deconvolution of an ESS passing through a power series, phase mismatch εm
Choosing sweep parameters in order to minimize the εm could be very restrictive and difficult when the order M increases. Novak proposes in [24] to overcome this problem by shifting, the signal z[n] by (∆m − Ceil(∆m )) before extracting each gm [n]: ω
j fs (Ceil(∆m )−∆m ) 1 Zm ( jω) = Z( jω)e
(25)
−150 0.5
1
1.5 Hz
2 4 x 10
Fig. 6. Hammerstein kernels phase diagrams of a power series after application of the matrix CorrB
2.4 Phase Mismatch The harmonic impulse responses (gm [n]) appear in the deconvolved signal z[n] at specific positions on the time axis. They can be extracted, starting at n0 − ∆m for the first sample. However, no matter how we adjust the ESS parameters, it is impossible to satisfy the condition ∆m ∈ N, ∀m ∈ [1, M] (except for the trivial case M = 1). Therefore, if we extract gm [n] from n0 − Ceil(∆m ) an error εm is created on the starting position as it could be seen at Fig. 8. As this error εm depends on the index ’m’, this leads to a 6
This phase mismatch mostly introduces some reconstruction errors at high frequencies. Fig. 9 shows the comparison between y1 and y2 if we don’t apply the phase mismatch correction. A significant difference between the two signals starts at 3000Hz and increases with frequency. 2.5 ESS Fade In/Out In order to avoid spectral leakage due to abrupt termination of the ESS, appropriate fade-in and fade-out procedures have been applied, according to [24]. If the loss of bandwith is expressed as a fraction of f1 and f2 , for the fadeIn and fadeout respectively, then: length f adeIn = R ∗ log(1 + loss) 1 − loss )+1−N length f adeOut = R ∗ log( f2 ∗ f1
(26)
J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
y1
1.2
y2 1
0.6 0.4 0.2 0
0
5
10
15 Sample [n]
20
25
Fig. 9. Comparison between y1 [n] and y2 [n] at 6000Hz when the phase correction from Eq. (25) is not applied (y2 [n] perfectly matches y1 [n] if the correction is applied)
2.6 Spreading of the Dirac Impulse We have assumed in Eq. (7) that the deconvolution of ss(m.φ [n]) is a delayed Dirac impulse δ [n − (n0 − ∆m )]. As the signal is band-limited, the Dirac impulse is in fact a Sync function δm0 . It follows that each measured kernel gm [n] is in fact: gm,measured [n] = gm [n] ~
δm0
1.2
y’1 y2
1
(27) 0.8
Moreover, as the signal ss(m.φ [n]) starts at frequency m. f1 instead of f1 , δm0 is more and more band limited as the order m increases. These limitations spread the kernels gm [n] around the positions n0 − ∆m .
Amplitude
Amplitude
0.8
this offset. Thus,we propose to cut the kernels gm [n] at n0 − ∆m − 1000, compute Hammerstein kernels and then delete the first thousand samples of the Hammerstein kernels. Fig. 11 shows that this method gives excellent results. In the case where the bandwidth is short (measurement of a tweeter for example) the Dirac impulses become larger and larger. Simply cutting the first thousand samples after the computation of the Hammerstein kernels could lead to a wrong representation. In this case, a compromise has to be found between the accuracy of the Hammerstein kernels and their possible use in real-time emulations. In the Toolbox described in section 5, the choice of the offset before the beginning of the Hammerstein kernels is therefore left open to the user, depending on the bandwidth and the real-time constraint of the user’s application. One can notice that the sampling rate used here is always 44100Hz as we try to ensure compatibility with basic sound cards. Nevertheless, the use of higher sampling rate allows a better compactness of the impulse response as the harmonic m extends its bandwidth to the frequency m. f2 or until the Nyquist frequency of the sampling rate is reached.
0.6 0.4 0.2
Table 1. Mean error around 1000Hz between y01 and y2 for different offsets cutting offset
-1000
-500
-100
-10
0
0 −0.2 0
50
100 Sample [n]
150
200
Fig. 10. Comparison between emulated signals cutting gm [n] at n0 − ∆m − 1000 (see y01 ) and at n0 − ∆m (see y2 ). 1.8
y1
1.6
y2
1.4 1.2 Amplitude
Extracting the kernels gm [n] from z[n], by starting at the exact positions ∆m can therefore be problematic. We prefer cutting the kernels gm [n] a few samples before that position to have a better reconstruction of the Hammerstein kernels. In practice, 1000 samples have been chosen and there is no need to take a larger number of samples (at fs = 44100Hz). For example, Fig. 10 shows the comparison between the emulated signal using kernels gm [n] cut at n0 − ∆m − 1000 (signal y01 ) and the emulated signal using kernels gm [n] cut at n0 − ∆m (signal y2 ). The signal y1 is better aligned to the signal y1 from Fig. 5. Table 1 shows the mean error between y1 and y2 over 2 periods around the frequency f=1000Hz for different cutting offsets.
1 0.8 0.6 0.4 0.2
Mean error
0.005
0.007
0.014
0.062
0.155
Note. Bandwidth = [5,22000]Hz with fs = 44100Hz
For real time emulations, introducing such a delay in the impulse responses could affect the output signal in the same way since the impulse responses will be convolved with the signal to emulate, this signal will be delayed by J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
0 0
50
100 Sample [n]
150
Fig. 11. Comparison between the output power series y1 and the emulated signal y2 calculated by cutting gm [n] at n0 − ∆m − 1000 and deleting the first thousand samples of hm [n]
7
Th.Schmitz AND J.J.Embrechts
PAPERS
2.7 Computing the Powers of the Input Signal Eq. (2) shows that implementing the Hammerstein model requires the powers of the input signal. Some precautions must be taken to avoid aliasing. In order to raise the input signal x[n] to the power m, the following procedure is used: r Resample the input signal at m. fs by zeros padding where fs is the original sample frequency. r Filter the signal in the initial Nyquist band fny = fs . 2 r Raise the obtained signal to the power m. r Filter the result in the Nyquist band fny . r Downsample the filtered signal at the fs frequency.
compares the signal y1 coming from the distortion device with the signal y2 computed with the HKISS method using 6 Hammerstein kernels. The two curves are very close. The mean error is : N−1
∑
y1 [n] | max(|y |) − 1
y2 [n] max(|y2 |) |
N
n=0
< 1%
(28)
The spectral contents of a pure tone (500Hz) passing through the distortion device or through the simulator (using 9 Hammerstein kernels) are compared in Table 2. A deviation smaller than 1.5db up to the 9th harmonic can be noticed. 0.4
y
3 MEASUREMENTS AND EMULATIONS
1
y2
0.3
In the precedent section, a set of precautions necessary to a better computation of the Hammerstein kernels has been presented. Taking these precautions (on phase error, phase mismatch, last sample, spreading of the Dirac impulse) leads to a Hammerstein Kernels Identification by Sine Sweep (HKISS) method. This section presents several tests made on real nonlinear audio devices in order to prove the efficiency of this method.
0.2 0.1 0 −0.1 −0.2 −0.3
3.1 Evaluation Algorithm We propose two methods to measure the accuracy of the emulated signal. 1. The first method, presented in Fig. 12, compares the output signal y1 coming from a nonlinear device driven by an ESS with its emulated signal y2 obtained using the Hammerstein kernels calculated from y1 . 2. The second method compares the spectral content of the signal coming from the nonlinear device and from the simulator when they are driven by a pure tone.
NL device
ss[n]
y1
under t est
ss[n] Comput e Volt erra
(.) M
kernels
−0.4 0
20
40
60
80
100
120
140
160
Fig. 13. Comparison between the sweep signals y1 and y2 around 1000Hz for the preamplifier distortion TSA15h Table 2. Harmonic amplitudes at the output of the TSA15h distortion for a pure tone input signal at 500Hz Harmonic order
Distortion effect(db)
Simulation(db)
1
92.2
92.2
3
65.3
64.9
5
56.4
56.1
7
51
50.8
9
46.5
45
Note. Only the odd harmonics are presented here as the distortion effect clips the signal in a symmetric way.
z[n]
hm [n] Nonlinear
ssm [n]
Convolut ion
y2
Fig. 12. Comparison between the signal coming from the nonlinear device under test y1 and its emulation y2
3.2 A Guitar Distortion Effect The first nonlinear audio device under test is the distortion effect coming from the guitar Tube Screamer Amplifier TSA15h [36] with its knob Gain at its maximum. Fig 13 8
3.3 Tubes Amplifier TSA15h One of the most popular application in the field of audio nonlinear modeling is the emulation of tubes amplifier so dear to guitarists. The emulation of the tubes amplifier TSA15h with all the knobs Gain at their maximum is presented in Figs. 14 ,15 & 16 around the frequencies 200, 1000 and 5000Hz respectively. The simulations result y2 (using 6 Hammerstein kernels) is very close to the measured signal y1 , especially at middle and high frequencies The comparison of the spectral content for a pure tone of 500Hz passing through the amplifier or through the simuJ. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
lator is slightly worse than for the distortion effect but still acceptable (see Table 3). Several reasons can explain these differences :
0.5
0
−0.5
−1 0
Finally, we have compared the signals coming from the amplifier and from the simulator using different numbers of Hammerstein kernels during the emulation (see Fig. 17). A good agreement can already be obtained with only few Hammerstein kernels (around five or six) which is interesting for real time implementations. 0.8
0
20
25
y1 y2(order = 1)
1
y2(order = 5) y2(order = 14)
Amplitude
Amplitude
0.2
10 15 Sample [n]
1.5
y2
0.4
5
Fig. 16. Comparison between the sweep signals y1 and y2 for the tubes amplifier TSA15h around 5000Hz
y1
0.6
y1 y2
Amplitude
r the nonlinearities created by the tube amplifier are more pronounced than those generated by the distortion effect (see Tables 2&3); r the tube amplifier has an output transformer whose hysteresis behavior can discard from a simple Hammerstein system; r the tubes power compression (the sagging effect [41]) is not a time invariant effect.
1
0.5
0
−0.2 −0.5 −0.4 −0.6 −0.8 0
−1 0 100
200 300 Sample [n]
400
Fig. 14. Comparison between the sweep signals y1 and y2 for the tubes amplifier TSA15h around 200Hz
0.8
y1 y2
0.6 0.4 Amplitude
20
40 60 Sample [n]
80
100
500
0.2 0
Fig. 17. Comparison between sinusoidal signals y1 and y2 for the tube amplifier TSA15h at 500Hz with different order of simulation
between y1 and y2 around 200Hz and 1000Hz respectively. Note that the system is strongly nonlinear. In this case, the emulation has been made with 14 Hammerstein kernels. A little overshoot can be observed on the emulation signal at 1000Hz. The mean error as defined by (28) is less than 2.3%.
−0.2 −0.4
Table 3. Harmonic amplitude at the output of the amplifier TSA15h for a pure tone input signal at 500Hz
−0.6 −0.8 0
20
40
60 Sample [n]
80
100
120
Fig. 15. Comparison between the sweep signals y1 and y2 for the tubes amplifier TSA15h around 1000Hz
3.4 Tubes Amplifier: Engl Retro Tube 50 E762 Another tubes amplifier has been tested, the Engl retro tube 50 E762 [42]. Figs. 18 & 19 present the comparison J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
Harmonic order
Tube amplifier(db)
Simulation(db
1
99.5
99.5
2
89
88.9
3
86.1
84.5
4
79.3
78.3
5
75.4
76.7
9
Th.Schmitz AND J.J.Embrechts
PAPERS
0.8
1.4
y
y
1
1
y2
0.6
2
1.2
0.4
1
0.2
Amplitude
Amplitude
y
0 −0.2
0.8 0.6 0.4
−0.4
0.2
−0.6 −0.8 0
100
200
300 Sample [n]
400
500
Fig. 18. Comparison between y1 and y2 for the tubes amplifier Engl at 200Hz. 0.6
0 0
600
y2
y1 y2
0.4
Amplitude
Amplitude
200
0.3
0 −0.2 −0.4
0.2 0.1 0
−0.6
−0.1 50
100 Sample [n]
150
200
Fig. 19. Comparison between y1 and y2 for the tubes amplifier Engl at 1000Hz.
4 DISCUSSION AROUND THE LIMITATIONS OF THE METHOD The HKISS method is pretty stable. Applying again the test algorithm of section 2 (power series model) up to the order 20 causes no problem, as it can be seen in Fig. 20. The model is also accurate for high frequencies (see Fig. 21). The remaining limitations are discussed in this section. 4.1 Order Limitation As the amplitude of the kernels gm [n] decreases when the order m increases, it finally reaches the background noise. This creates a maximum order limit. Improving the signal to noise ratio (SNR) is thus important if high order emulation is needed. It could be convenient to take several measurements of a device and averaging the results in order to improve the SNR. However, we must be sure that the system is perfectly time invariant. The measurement of a loudspeaker at high power for example could be difficult as the heating of the coil changes its characteristics. 4.2 Low Frequency Limitation As the ESS is band limited in the frequency interval [ω1 , ω2 ], its mth harmonic is band limited in [m.ω1 , m.ω2 ] 10
150
0.5
0.2
−0.8 0
100 Sample [n]
Fig. 20. Comparison between the sweep y1 and y2 for a power series up to the 20th order around 1000Hz
y1
0.4
50
−0.2 0
5
10
15
20 25 Sample [n]
30
35
40
Fig. 21. Comparison between the sweep y1 and y2 for a power series up to the 20th order around 16000Hz
(the upper bound could be ω2 if we initially chose ω2 equal to the Nyquist frequency). It means that the deconvolution by the inverse filter is also band limited in [m.ω1 , m.ω2 ] which gives a band limited Dirac impulse δ˜m . Eq. (7) can be rewritten as: ss (m.φ [n]) ~ ss(φ [n]) = D.δ˜m [n − (n0 − ∆m )]
(29)
In the frequency domain, the impulse δ˜m delayed by ∆m is a step function I( jω) on the domain [m.ω1 , mω2 ] (see Fig. 22). Therefore : SSm ( jω).SS( jω) = I[m.w1 ,m.w2 ] ( jω)
(30)
In order to understand the impact on the Hammerstein kernels computations, let’s take Eq. (13) again and convolve this output signal y(t) with the ESS inverse filter using Eq. (30). We obtain: 3A2 G1 ( jω) = I H3 ( jω).I[w1 ,w2 ] [w1 ,w2 ] .H1 ( jω) + 4 A G2 ( jω) = j.e− jB .H2 ( jω).I[2.w1 ,2.w2 ] 2 2 −A G3 ( jω) = e− j2B .H3 ( jω).I[3.w1 ,3.w2 ] 4 (31) J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
100
δ1 δ3
Amplitude db
95
δ6
90 f1
3.f
1
6.f
1
85
80
75
70
1
2
10
10 Frequency in Hz
Fig. 22. FFT of δ˜m with m ∈ {1, 3, 6}
Fig. 23. Comparison between the first Hammerstein kernel coming from a power series up to the order 6 and its reconstruction by the HKISS method.
By inverting the relation to have the Hammerstein kernels, we have: −4 H3 ( jω).I[3.w1 ,3.w2 ] = 2 e j2B .G3 ( jω) A −2 jB H2 ( jω).I[2.w1 ,2.w2 ] = j.e .G2 ( jω) A 3A2 H ( jω).I H3 ( jω).I[w1 ,w2 ] 1 [w1 ,w2 ] = G1 ( jω) − 4 (32) This last expression indicates that to calculate H1 ( jω).I[w1 ,w2 ] , H3 ( jω).I[w1 ,w2 ] must be known but we only have H3 ( jω).I[3.w1 ,3.w2 ] . However the formula is still correct if only the domain [3ω1 , 3ω2 ] is considered. Therefore, we can calculate the Hammerstein kernels for a nonlinear system of order M only if: ωvalid ≥ M.ω1
(33)
To illustrate the problem, the Fig. 23 presents the comparison in frequency domain between the first Hammerstein kernel of a power series (up to the 6th order) and the first Hammerstein kernels reconstructed with the HKISS method. The reconstructed one does not correspond to the initial one until it reaches the frequency f = 6 ∗ f1 .A way to get around this problem is to choose the ESS parameter ω1 according the desired bandwidth. For example, if we would like to emulate a 10 order nonlinear system on the bandwidth [100, 20000]Hz, we must take f1 = 10Hz. 4.3 Hammerstein Kernels Input Level Dependency As discussed in section 1.4, the accuracy of the emulated signal depends on the match between the input signal amplitude and the amplitude α of the sine sweep used to measure and compute the Hammerstein kernels. During the emulation test of the amplifier TSA15h, we have noticed that the Hammerstein kernels Hmα ( jω) have the same shape for a given order m but have different relative amplitudes (especially for higher-order kernels). For example, the first and the fifth Hammerstein kernels of the device TSA15h have been measured using an ESS with full amplitude (α = 1) and half amplitude (α = 0.5). Their Fourier J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
transforms are compared in Fig. 24. It involves that we can calculate the Hammerstein kernels Hmα ( jω) from Hm1 ( jω) by changing the amplitude factor A appearing in Eq. (16). The Fig. 25 shows an example of an adapted Hm1 ( jω) obtained by changing the A factor, in order to fit the Hm0.5 ( jω) kernels. With these adapted Hammerstein kernels based on the measure at full amplitude, we can emulate the signal at half amplitude (see Fig. 26) with less than 4.7% of mean error. The parameter A has been manually adjusted for our example. To avoid having a large set of impulse responses Hmα ( jω) corresponding to each level in the simulator (space memory) and to avoid some difficulties with a variable convolution kernels, the Hammerstein model equation can be weighted as follows: M
y[n] =
∑
m=1
1 [A(α)]m−1
h1m [n] ~ xm [n]
(34)
A way to determine the function A(α) is to find for each level α the A factor that minimize the areas between the Fourier Transform of the Hammerstein kernels measured at the level α and the A corrected Hammerstein kernels measured at the level α = 1. Unfortunately, some devices could have their Hammerstein kernels Hmα too different from each other to use this method (in general strong nonlinear devices). The Hammerstein kernels can then be approximated for a continuous range of levels by means of an interpolation procedure as proposed by Tronchin and Coli in [35] where each Hammerstein kernel hαm [n] needs its own Am (α) function. 5 A MATLAB TOOLBOX The practical implementation of Hammerstein identification technique is complex and may be subject to a lot of possible errors. Hence, we have proposed a Matlab Toolbox (see Fig. 27) for the Hammerstein identification base on the methods presented in the previous sections. The toolbox is composed of three parts: the ESS generation, where the signal to send to the nonlinear device under test 11
Th.Schmitz AND J.J.Embrechts
PAPERS
20
h1
Full
h5Full
15
h10.5 h50.5
10 Amplitude db
DUT and the emulation of the nonlinear DUT. The toolbox can be downloaded in [43] and more information about it could be found in [44].
5 0 −5 −10 −15 2 10
3
4
10 Frequency in Hz
10
Fig. 24. FFT of the first and fifth Hammerstein kernels of the distortion TSA15h measured with an ESS with full and half amplitudes 20
h1Full h5Full
15
h10.5 h50.5
Amplitude db
10
Fig. 27. The Hammerstein Identification Matlab Toolbox 5 0
5.1 Power series example A short example using the toolbox to compute the Hammerstein kernels of a power series up to the order six is given here.
−5 −10 −15 2 10
3
4
10 Frequency in Hz
10
Fig. 25. FFT of the first and fifth Hammerstein kernels of the distortion TSA15h measured with an ESS with full amplitude and adapted to match the kernels of half amplitude changing the A factor. 1
y1 y2
Amplitude
0.5
0
−0.5
−1 0
50
100 Sample [n]
150
200
Fig. 26. Comparison between an ESS of amplitude 0.5 passing through the distortion effect TSA15h and it simulation (y2 ) by means of Hammerstein kernels measured for an ESS of amplitude 1 but corrected with the A factor (see text).
(DUT) is generated, the calculation of the Hammerstein kernels through the signal recovered at the output of the 12
5.1.1 Sine sweep generation The tab Sweep generation is the first step to measure a nonlinear device. It allows the generation of the ESS where the parameters are : f1 the initial frequency of the ESS, f2 the final frequency of the ESS, N the length of the ESS in samples, fade in (or out) is the percentage of the value of f1 (or f2 ) until a fade in (or out) is allowed. For example if f1 = 10Hz and fade in = 0.1, it means that a fade in will be apply between the frequency f1 = 10Hz and ten percents more so 10+10*0.1 = 11Hz. For the fade out it will be apply between the frequency f2 and f2 minus ten percents. These values are shown in the frame Bandwidth. The generate button exports the ESS in a .wav format. 5.1.2 Hammerstein kernels calculation Once the ESS signal has been sent to the DUT, the y signal at the output of the DUT has to be convolved with the inverse filter of the ESS. In this example, the output signal of a power series until the 6th order named yPowerSeries is provided. It can be loaded in the toolbox by pushing the Load y response button. After having pressed the Deconvolve button, the z impulse is obtained as has depicted at Fig.28. Before cutting z into gm kernels, we must specify the following parameters : Nb Kernels is the required number of Hammerstein kernels for the emulation and the Length Kernels is the desired size of the impulse responses. To compute the Hammerstein kernels, we also J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
have to specify the Decaycut which is the number of samples kept before the theoretical beginning of the Hammerstein kernels. This number has to be comprised in the interval [-1000 0]. 5.1.3 Emulation of the DUT The third tab allows the emulation of the DUT by convolving the Hammerstein kernels with any chosen signal in a .wav format. For this example, we can load the sweep.wav signal and convolve it with the Hammerstein kernels to emulate the power series giving back a signal similar to yPowerSeries.
Fig. 28. The Hammerstein kernels identification tabs
6 CONCLUSION Many researches have been done on Volterra kernels identification in the last few years. As the impulse response characterizes a linear system very precisely, we expect to extend the method to nonlinear systems with an equivalent success. As we have shown, the cascade of Hammerstein models can lead to a great sensitivity to implementation inaccuracies, especially when the order M increases. In this paper, we have highlighted different kinds of errors and we have shown how they can affect the emulation of a nonlinear system. Moreover we have proposed a way to get around each problem leading us to an accurate emulation of three nonlinear systems: the Tubes-Screamer like overdrive effect, the TSA15h tube amplifier and the Engl Retro Tube 50 E762 amplifier. We have also described the limitations of the model (i.e. maximum order, low frequency limitation, input level dependency). The method seems well suited for real time emulation of audio nonlinear systems. Finally, we have proposed a solution in Eq (34) for the input amplitude level dependency which seems to be valid for some nonlinear devices. However, for strongly nonlinear systems, further research is still needed to solve this practical problem. J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
7 REFERENCES [1] T. Schmitz, J.J. Embrechts, ”Nonlinear Guitar Loudspeaker Simulation,” presented at the 134th Convention of the Audio Engineering Society (2013 May),convention eBrief 96. [2] W. Klippel, ”Loudspeaker NonlinearitiesCauses, Parameters, Symptoms,” presented at the 119th Convention of the Audio Engineering Society (2005 October). [3] I. Cohen, T. Helie, ”Simulation of a Guitar Amplifier Stage for Several Triode Models: Examination of Some Relevant Phenomena and Choice of Adapted Numerical Schemes,” presented at the 127th Convention of the Audio Engineering Society (2009 October) [4] M.J. Kemp, ”Analysis and Simulation of NonLinear Audio Processes Using Finite Impulse Responses Derived at Multiple Impulse Amplitudes,” presented at the 106th Convention of the Audio Engineering Society (1999 May) [5] A.V. Oppenheim, A.S. Willsky, Signals and Systems (Prentice Hall, 1996). [6] M.Schetzen, The Volterra & Wiener Theory of Nonlinear Systems (John Wiley & Sons, 1980). [7] J.S. Bendat, Nonlinear System Techniques and Applications, New York: Wiley, 1998. [8] M. Schoukens, K. Tiels, M. Ishteva, & J. Schoukens, ”Identification of Parallel Wiener-Hammerstein Systems with a Decoupled Static Nonlinearity,” presented at the 19th World Congress of the International Federation of Automatic Control, Cape Town, South Africa, pp. 24– 29, (2014 August). http://dx.doi.org/10.3182/ 20140824-6-ZA-1003.00496. [9] M. Schoukens, G. Vandersteen, & Y. Rolain, ”An Identification Algorithm for Parallel Wiener-Hammerstein Systems,” presented at the 52th Annual Conference on Decision and Control (CDC), Firenze,(2013 December). http://dx.doi.org/10.1109/CDC.2013. 6760659. [10] F. Thouverez, & L. Jezequel, ”Identification of NARMAX Models on a Modal Base,” Journal of Sound and Vibration, Vol 189, No. 2, pp. 193–213, (1996). http://dx.doi.org/10.1006/jsvi.1996. 0015 [11] O. Nelles, Nonlinear System Identification: From Classical Approaches to Neural Networks and Fuzzy Models. Berlin, Germany: Springer-Verlag, (2001). http:// dx.doi.org/10.1007/978-3-662-04323-3. [12] Y.W. Chen, S. Narieda, & K. Yamashita”Blind Nonlinear System Identification Based on a Constrained Hybrid Genetic Algorithm,” IEEE Transaction on Instrumentation and Measurement, Vol 52, No. 3, pp. 898–902, (2003). [13] H. Sorenson, Kalman Filtering: Theory and Application, Montvale, NJ: IEEE Press, (1985). [14] O. Cappe, S.J. Godsill, & E. Moulines, ”An Overview of Existing Methods and Recent Advances in Sequential Monte Carlo,” Proceeding of IEEE, Vol 95, No. 5, pp 899–924, (2007).
13
Th.Schmitz AND J.J.Embrechts
[15] T.Ogunfunmi, Adaptive Nonlinear System Identification: The Volterra and Wiener Model Approaches Springer, (2007). http://dx.doi.org/10.1007/ 978-0-387-68630-1. [16] A. Novak,L. Simon, P. Lotton, & F. Kadlec, ”A New Method for Identification of Nonlinear Systems Using Miso Model with Swept-Sine Technique: Application to Loudspeaker Analysis,” presented at the 124th Convention of the Audio Engineering Society (2008 May), convention paper 7441. [17] M. R´ebillat, R. Hennequin, E. Corteel, & B. Katz, ”Identification of Cascade of Hammerstein Models for the Description of Nonlinearities in Vibrating Devices,” Journal of Sound and Vibration, Vol. 330, Issue 5, pp. 1018–1038 (2011 February). http://dx.doi.org/ 10.1016/j.jsv.2010.09.012. [18] A. Farina, A. Bellini, & E. Armelloni, ”Nonlinear Convolution: A New Approach for the Auralization of Distorting Systems,” presented at the 110th Convention of the Audio Engineering Society (2001 May), convention paper 5359. [19] A. Farina, ”Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique,” presented at the 108th Convention of the Audio Engineering Society (2000 February), convention paper 5093. [20] T. Schmitz, J.J. Embrechts, ”Improvement in Nonlinear Guitar Loudspeaker Sound Reproduction,” presented at the 21th International Conference on Systems, Signals and Image Processing (2014 May), international conference on IEEE. [21] L. Tronchin, ”The Emulation of Nonlinear TimeInvariant Audio Systems with Memory by Means of Volterra Series,” J. Audio Eng. Soc., Vol. 60, No. 12, (2012 December) [22] A. Novak, L. Simon, & P. Lotton, ”Analysis, Synthesis, and Classification of Nonlinear Systems Using Synchronized Swept-Sine Method for Audio Effects,” Eurasip Journal on Advances in Signal Processing, Vol. 2010:793816, (2010 July). http://dx.doi.org/10. 1155/2010/793816. [23] A. Novak, L. Simon, ”Nonlinear System Identification Using Exponential Swept-Sine Signal,” presented at the IEEE Transactions on Instrumentation and Measurement, Vol. 59, No 8, pp. 2220–2229, (2009 August). http://dx.doi.org/10.1109/TIM. 2009.2031836. [24] A. Novak, P. Lotton, & L. Simon, ”Synchronized Swept-Sine: Theory, Application and Implementation,” J. Audio Eng. Soc., Vol. 63, No. 10, (2015 October). http: //dx.doi.org/10.17743/jaes.2015.0071 [25] J. Vanderkooy, & S. P. I. Thomson, ”Harmonic Distortion Measurement for Nonlinear System Identification,” presented at the 140th Convention of Audio Engineering Society (2016 June), convention paper 9497. [26] G. Palm, ”On Representation and Approximation of Nonlinear Systems,” Journal Biological Cybernetics, Vol. 34, Issue 1, pp. 49–52 (1979 September). http: //dx.doi.org/10.1007/BF00336857.
14
PAPERS
[27] M. R´ebillat, R. Hennequin, E. Corteel, B. Katz, ”Identification of Cascade of Hammerstein Models for the Description of Nonlinearities in Vibrating Devices,” Journal of Sound and Vibration, Vol. 330, no. 5 (2011). http: //dx.doi.org/10.1016/j.jsv.2010.09.012. [28] M. Holters, T. Corbach, & U. Zlzer, ”Impulse Response Measurement Techniques and their Applicability in the Real World,” presented at the 12th Int. Conference on Digital Audio Effects (2009 September). [29] Q. Meng, D. Sen, S. Wang, & L. Hayes, ”Impulse Response Measurement with Sine Sweeps and Amplitude Modulation Schemes,” 2nd International Conference on Signal Processing and Communication Systems, (2008 December). http://dx.doi.org/10.1109/ ICSPCS.2008.4813749. [30] S. Norcross, G. Soulodre, & M. Lavoie, ”Evaluation of Inverse Filtering Techniques for Room/Speaker Equalization,” presented at the 113th Convention of the Audio Engineering Society (2002 October),convention paper 5662. [31] O. Kirkeby, & P. Nelson, ”Digital Filter Design for Inversion Problems in Sound Reproduction,” J. Audio Eng. Soc., Vol. 47, No. 7/8, (1999 July/August). [32] D. Zwillinger, Standard Mathematical Tables and Formulae, 32nd ed. CRC Press, (2011). [33] S. QIN, Y. QIN, Y. MAO, ”Accurate Instantaneous Frequency Estimation with Iterated Hilbert Transform and Its Application,” presented at the Proceedings of the 7th WSEAS International Conference on SIGNAL PROCESSING, ROBOTICS and AUTOMATION (ISPRA’08),(2008 February) [34] M.R. Spiegel, S. Lipschutz, & J. Liu, Mathematical Handbook of Formulas and Tables, 3rd ed. Schaum’s Outlines Series, McGraw-Hill eBook, (2009). [35] L. Tronchin, V.L. Coli, ”Further Investigations in the Emulation of Nonlinear Systems with Volterra Series,” J. Audio Eng. Soc., Vol. 63, No. 9, (2015 September). [36] Tube Screamer Amplifier TSA15H: http://www.ibanez.com/products/amp page15.php?year =2015&area id=3&cat id=5&series id=141&l id=4 [37] D. Ciric, M. Markovic, M Mijic, & S. SumaracPavlovic, ”On the Effects of Nonlinearities in Room Impulse Response Measurements with Exponential Sweeps,” Journal Applied Acoustics, Vol. 74, pp. 375–382, (2013). [38] A. Farina, ”Advancements in Impulse Response Measurements by Sine Sweeps,” presented at the 122th Convention of Audio Engineering Society (2007 May), convention paper 7121. [39] A. Torras-Rossel, F. Jacobsen, ”A New Interpretation of Distortion Artifacts in Sweep Measurements,” J. Audio Eng. Soc., Vol. 59, No. 5, (2011 May). [40] Scip Solver: http://scip.zib.de/ [41] J. Macak, & J. Schimmel, ”Real-Time Guitar Tube Amplifier Simulation Using an Approximation of Differential Equations,” presented at the proceeding of the 13th Int. Conference on Digital Audio Effects, Graz, Austria, (2010 September), DAFx-10.
J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
PAPERS
Hammerstein Kernels Identification by Means of a Sine Sweep Technique
[42] Engl Retro Tube E762 Amplifier: http://www.englamps.de/index.php?id=36&tx ddfproducts pi1[uid]=119&iPhoneBypass=1 [43] Matlab Toolbox for Identification of the Hammerstein Kernels: https://github.com/TSchmitzULG/HKISS
[44] T. Schmitz, & J.-J. Embrechts, ”A new toolbox for the identification of diagonal Volterra kernels allowing the emulation of nonlinear audio devices,” presented at the proceeding of the 22nd International Congress on Acoustics, Buenos-Aires, Argentina, (2016 September), ICA.
THE AUTHORS
Thomas Schmitz
Thomas Schmitz received the degree in Electrical Engineering (2012) from from the University of Liege (ULg). His final project focused on the emulation of a electrodynamics loudspeaker including its nonlinear behavior. He is presently a Ph.D. student in Laboratory for Signal and Image Exploitation (INTELSIG) research unit of the Electrical Engineering and Computer Science (EECS) department , University of Liege, Belgium. His research interests are on signal processing, nonlinear modeling and real time emulation of guitar audio systems.
r
J. Audio Eng. Sco., Vol. 1, No. 1, 2015 October
J-J. Embrechts Pr. J-J. Embrechts received the degree in Electrical Engineering (1981) and the Ph.D. degree (1987) from the University of Liege (ULg). Since 1999, he is a professor at the University of Liege, in the Department of Electrical Engineering and Computer Science, teaching acoustics, electroacoustics, audio and video engineering and lighting techniques. He is the author and co-author of more than 120 journal and conference papers in the field of audio signal processing, room acoustics, auralization and lighting techniques.
15