© T. Lebrun and T. Hélie, Published by EDP Sciences, 2020

1 Introduction

When a loudspeaker diaphragm vibrates, its displacement modulates the propagation time of the generated acoustic waves, between its surface and a listening point. This moving-boundary phenomenon, similar to the Doppler effect, is usually considered to be negligible in practice, based on distortion measurements [1–4]. This is especially relevant for narrow-band or harmonic excitations. However, distortion increases for more complex signals, due to a significant intermodulation between the low-frequency content (large displacements) and the higher-frequency content (large accelerations). It has been shown [5, 6] that this phenomenon can be a dominant source of intermodulation distortion (compared to other sources such as force factor, suspension, etc.) for full-range speakers at high frequencies.

Initial investigations were carried out by Beers and Belar [7], who measured sound distortion generated by a two-tone vibrating diaphragm and derived a criterion for the evaluation of the Doppler distortion magnitude. Various models have been proposed to describe this phenomenon: (i) van Wulfften Palthe [1] adopted a pulsating sphere as a loudspeaker radiation approximation and performed calculation of intermodulation distortion, considering both the moving-boundary effect (Doppler) and the nonlinear propagation; (ii) Braun [8] established a time-domain formulation of the problem, followed by Butterweck [9] who derived a series solution for plane waves propagation generated by a vibrating piston in an infinite baffle; (iii) Zóltogórski [10] added nonlinear propagation effects to Butterweck’s model, and presented a preliminary design of an anti-Doppler filter. Another correction system has been proposed by Klippel [11], in which the electrical input is delayed to oppose the time-shift introduced by the moving boundary.

The present paper models this distortion effect in a similar way as in [9, 10], using a regular perturbation method. It examines its significance and addresses its correction based an algorithm built on inverse Volterra kernels.

The content is organised as follows. First, the considered model and both the direct and inverse problems are described in Section 2. Then the direct problem is solved in Section 3 by: (a) converting of the (Lagrangian) position of the moving source into its equivalent (Eulerian) velocity field at a fixed position; (b) deriving the acoustic pressure radiated from this velocity field. Finally, the inverse (anti-Doppler) problem is examined in Section 4, that is, the derivation of a piston displacement that generates a targeted Eulerian velocity field. Conclusions are drawn in Section 5 on the relevance of Doppler distortion and anti-Doppler efficiency in the case of realistic loudspeaker diaphragm motions.

2 Problem statement

Consider a baffled circular piston of radius R ₀, localised at position z = 0 (see Fig. 1a). For an Eulerian velocity field excitation V ₀(t), the corresponding far-field pressure p(z, t) on the symmetry axis at position z > 0 is given in the time domain by [12]

$$p(z,t)=\frac{{\rho }_0 R_0^2}{2z} \frac{\mathrm{d}}{\mathrm{d}t}{V}_0\left(t-\frac{z}{c_0}\right),$$(1) where ρ ₀ is the air density and c ₀ is the speed of sound.

Figure 1 Description of the physical model composed of (a) on-axis baffled piston radiation for the evaluation of the sound pressure and (b) conversion of the piston displacement into an equivalent Eulerian velocity source $V_0\left(t\right)$.

The core of this study focuses on the (nonlinear) mapping between the rigid piston displacement and the equivalent Eulerian velocity field V ₀(t), as described in Figure 1b. This issue is addressed under the following assumptions:

• H1 (geometry): The piston displacement ξ(t) is small compared to its radius (|ξ(t)| < ξ _max ≪ R ₀)

• H2 (fields equivalence approximation): In the domain Ω = [−ξ _max, ξ _max], the particle velocity ξ′(t) at the piston position ξ(t) results from the conservative plane wave propagation of an Eulerian velocity field V ₀(t) at z = 0

• H3 (waveform): the piston displacement ξ is a smooth function and no shockwave propagates in Ω

• H4 (radiation): the piston radiates into a semi-infinite space (no backward wave)

Following (H1–H2–H4), the particle velocity $v$ is described by a forward wave V ₀

$$v(z,t)=V_0(t-z/c_0).$$(2) and satisfies a moving-boundary condition at the piston position that reads

$$V_0(t-\xi (t)/c_0)={\xi }^{\mathrm{\prime }}\left(t\right).$$(3)

Moreover, condition H3 implies that the Eulerian velocity field V ₀ must be a regular function.

Let us denote $\mathcal{D}$ the operator that converts the piston displacement into its equivalent Eulerian velocity (see Fig. 2), such that $V_0\left(t\right)=\mathcal{D}\left[\xi \right]\left(t\right)$. This study examines the following problems:

Figure 2 Block diagram of the system under consideration, with input $\xi \left(t\right)$ and output $p(z,t)$.

• P1 (Direct problem): Is-it possible to derive a solver for $\mathcal{D}$ that satisfies (2)–(3)?

• P2 (Inverse problem): What is the corrected displacement ξ _c(t) to provide to $\mathcal{D}$ so that its output is the target $V_0\left(t\right)={\xi }^{\prime }\left(t\right)$?

The direct problem P1 is tackled in Section 3 and the inverse problem P2 is examined in Section 4.

3 Direct problem: simulation of the Doppler distortion

This section addresses the direct problem (Doppler effect) for smooth excitations: in Section 3.1, an exact solution of the Eulerian velocity V ₀ is determined for a linear propagation, based on the method of characteristics; in Section 3.2, a truncated series expansion that approximates the Eulerian velocity V ₀ by a closed-form solution is derived for simulation purposes, based on a perturbation method; the distortion is evaluated on the radiated pressure in Section 3.3; finally Section 3.4 includes the influence of nonlinear acoustic wave propagation.

3.1 Method of characteristics

This part addresses the existence and uniqueness of regular solutions of P1, based on the method of characteristics. It provides a regular expression of the velocity waveform V ₀ under a necessary and sufficient condition on ${\mathcal{C}}^1$-regular functions ξ. In the following,

• “a” refers to quantities related to the arrival time of the acoustic wave at the observing point,

• “d” denotes quantities related to the departure time of the wave from the piston position.

Definition 3.1 (Characteristic)

Let us define the regular functions

$$\begin{array}{c}{\tau }_{\xi }^a: \mathbb{R}\mapsto \mathbb{R}\\ t\mapsto t-\frac{\xi \left(t\right)}{c_0}\end{array}$$(4) and

$$\begin{array}{ll}{K}_{\xi }: {\mathbb{K}}_{\xi }^d& \mapsto {\mathbb{K}}_{\xi }^a\\ (z,t)& \mapsto (z,{\tau }_{\xi }^a(t)+z/c_0),\end{array}$$(5) where the domain ${\mathbb{K}}_{\xi }^d$ and codomain ${\mathbb{K}}_{\xi }^a$ are

$$\begin{array}{lll}{\mathbb{K}}_{\xi }^d& =& \left\{\right(z,t)\in {\mathbb{R}}^2 s.t.\hspace{0.5em}z\ge \xi (t\left)\right\},\\ {\mathbb{K}}_{\xi }^a& =& \left\{K_{\xi }\right(z,t) \mathrm{for} (z,t)\in {\mathbb{K}}_{\xi }^d\}.\end{array}$$(6)

The characteristic lines K _ξ defined above are depicted in Figure 3. They represent the linear propagation of the acoustic wave over space and time. A departure time t ^d is mapped to an arrival time t ^a at a given position Z through K _ξ . One should note that

1. the propagation is conservative so that the particle velocity is constant on the characteristic line described by K _ξ (r, t ^d) when z varies, for given a departure time t ^d. For the boundary condition (3), this observation yields

2. the departure and arrival times are related by $t^{\mathrm{a}}=t^{\mathrm{d}}+\frac{Z-\xi \left(t^{\mathrm{d}}\right)}{c_0}$. This can be written using the characteristic equation (6),

Figure 3 Illustration of the mapping from departure coordinates $\left(\xi \right(t^{\mathrm{d}}),t^{\mathrm{d}})$ to arrival coordinates $(Z,t^{\mathrm{a}})$ through the characteristic line (described by $K_{\xi }(z,t^{\mathrm{d}})$ when $z$ varies), for a sinusoidal piston motion. Note the variation in the propagation time $t^{\mathrm{a}}-t^{\mathrm{d}}=\frac{Z-\xi \left(t^{\mathrm{d}}\right)}{c_0}$.

$$v(Z,t^{\mathrm{a}})={\xi }^{\prime }\left(t^{\mathrm{d}}\right),$$(7)

$$(Z,t^{\mathrm{a}})=K_{\xi }(Z,t^{\mathrm{d}}).$$(8)

It is clear from these observations that the particle velocity v at any space-time coordinates $(Z,t^{\mathrm{a}})\in {\mathbb{K}}_{\xi }^{\mathrm{a}}$ can be computed from (7), if the departure time t ^d is known. However, equation (8) shows that the existence of t ^d is conditionned by the existence of $K_{\xi }^{-1}$, the reciprocal of K _ξ . This point is examined in the property detailed in Appendix A, which yields a condition on the Mach number of the piston velocity,

$$\left|{\xi }^{\mathrm{\prime }}\right(t\left)\right|<c_0.$$(9)

Finally, the solution of P1 is given in Theorem 3.2. It provides the regular waveform V ₀, solution of the propagation equation (2) and the boundary condition (3), under the Mach number condition (9).

Theorem 3.2 (Regular solution)

Let $\xi$ be a ${\mathcal{C}}^{n+1}$ -regular function and suppose that the Mach number condition (9) is fulfilled. Then the regular solution of (2)–(3) is given by the ${\mathcal{C}}^n$ -regular function

$$\begin{array}{ll}{V}_0:{\mathbb{R}}^{+}& \mapsto \mathbb{R}\\ t& \mapsto {\xi }^{\prime }\circ {\tau }_{\xi }^d\left(t\right)\end{array}$$(10) where ${\tau }_{\xi }^d$ is the reciprocal of ${\tau }_{\xi }^a$ .

Proof. From the property detailed in Appendix A, $K_{\xi }$ is a ${\mathcal{C}}^n$-diffeormorphism, so that the reciprocal $K_{\xi }^{-1}$ exists and is ${\mathcal{C}}^n$-regular. Moreover,

$$K_{\xi }^{-1}\left(\xi \right(t),t)=\left(\xi \right(t),t),$$(11) since (ξ(t), t) is a fixed point of K _ξ . The reciprocal function $K_{\xi }^{-1}$ can be developped

$$\left(\xi \left(t\right),{\tau }_{\xi }^{\mathrm{d}}\left(t-\frac{\xi \left(t\right)}{c_0}\right)\right)=\left(\xi \left(t\right),t\right),$$(12) where the second component of (12) reads

$${\tau }_{\xi }^{\mathrm{d}}\left(t-\frac{\xi \left(t\right)}{c_0}\right)=t.$$(13)

Finally,

$$V_0(t-\xi (t)/c_0)={\xi }^{\prime }\circ {\tau }_{\xi }^{\mathrm{d}}\left(t-\frac{\xi \left(t\right)}{c_0}\right)={\xi }^{\prime }\left(t\right),$$(14) therefore the regular solution (10) verifies the boundary condition (3). Moreover, ${\xi }^{\prime }\circ {\tau }_{\xi }^{\mathrm{d}}$ is ${\mathcal{C}}^n$-regular because ${\xi }^{\prime }$ and ${\tau }_{\xi }^{\mathrm{d}}$ are ${\mathcal{C}}^n$-regular, which concludes the proof. □

This theorem can be summarized in equations (15) and (16), with the expression of the Eulerian source

$$V_0\left(t\right)={\xi }^{\prime }\circ {\tau }_{\xi }^{\mathrm{d}}\left(t\right),$$(15) where ${\tau }_{\xi }^{\mathrm{d}}\left(t\right)$ maps the current time t at z = 0 to the departure time of the wave at the piston position, satisfying the implicit equation

$${\tau }_{\xi }^{\mathrm{d}}\left(t\right)=t+\frac{\xi \circ {\tau }_{\xi }^{\mathrm{d}}\left(t\right)}{c_0}.$$(16)

The velocity source V ₀(t) cannot be directly computed from the piston displacement ξ(t) because of the implicit equation (16). The derivation of a direct solver for (15) and (16) is addressed in the next part.

3.2 Regular perturbation method

A direct solver is established in this subsection, assuming a regular function $\xi \in {\mathcal{C}}^{\mathrm{\infty }}$ and $\left|{\xi }^{\mathrm{\prime }}\right(t\left)\right|<c_0$, so that ${\tau }_{\xi }^{\mathrm{d}}\in {\mathcal{C}}^{\mathrm{\infty }}$ from Theorem 3.2. First, function ${\tau }_{\xi }^{\mathrm{d}}$ is substituted for the ${\mathcal{C}}^{\mathrm{\infty }}$-regular function $ϵ:t\mapsto {\tau }_{\xi }^{\mathrm{d}}\left(t\right)-t$. Then equations (15) and (16) can be reformulated

$$V_0\left(t\right)={\xi }^{\prime }(t+ϵ(t\left)\right),$$(17)

$$ϵ\left(t\right)=\frac{\xi (t+ϵ(t\left)\right)}{c_0}.$$(18)

Equation (17) highlights that $ϵ\left(t\right)$ acts as a perturbation of the standard (linear) solution V ₀(t) = ξ′(t).

Now the implicit equation (18) is solved using the regular perturbation method. The piston displacement function ξ is marked with an auxiliary amplitude α > 0, so that we define ξ(t) = αu(t). Then we seek a power series solution of $ϵ\left(t\right)$ of the form

$$ϵ\left(t\right)={\displaystyle \sum _{n=0}^{\mathrm{\infty }}}\mathrm{}{ϵ}_n\left(t\right){\alpha }^n.$$(19)

Calculations, developed in Appendix B, yield the exact series solution of the form

$$V_0\left(t\right)={\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathrm{}{v}_n^D\left(t\right),$$(20)in which each term $v_n^D\left(t\right)$ are defined as follows,

$$v_n^D\left(t\right)=\{\begin{array}{ll}{\xi }^{\mathrm{\prime }}\left(t\right)& \mathrm{if} n=1,\\ {\displaystyle \sum _{k=1}^{n-1}}\mathrm{}\frac{{\xi }^{(k+1)}\left(t\right)}{k!}\sum _{i\in {\mathcal{C}}_{n-1,k}}\mathrm{}{ϵ}_{i_1}\left(t\right)...{ϵ}_{i_k}\left(t\right)& \mathrm{otherwise},\end{array}$$(21)where ξ ⁽ⁱ⁾ stands for the ith derivative of ξ and the ϵ _i is defined recursively in Appendix B.

In the sequel, computation of V ₀(t) is carried out by truncating the series at N = 3, that appears sufficient to capture the Doppler distortion effect for typical loudspeaker diaphragm motions. The expression of the truncated series is

$$\begin{array}{ll}{V_0\left(t\right)|}_{N=3}& =\underset{n=1}{\underbrace{{\xi }^{\prime }\left(t\right)}}+\underset{n=2}{\underbrace{\frac{1}{c_0}\xi \left(t\right){\xi }^{″}\left(t\right)}}\\ & +\underset{n=3}{\underbrace{\frac{1}{c_0^2}\left(\frac{1}{2}{\xi }^{\mathrm{″\prime }}\left(t\right)\xi (t{)}^2+{\xi }^{″}\left(t\right){\xi }^{\prime }(t\left)\xi \right(t)\right)}}.\end{array}$$(22)

The term n = 1 corresponds to the linear solution without taking into account the Doppler effect: the Eulerian field V ₀(t) equals the piston velocity. Considering the term n = 2, the product ξ(t)ξ″(t) will be of high amplitude if the piston velocity signal is composed of low frequency components (amplified by ξ) together with high frequencies (amplified by ξ″), therefore generating intermodulation distortion in V ₀(t).

A realisation structure of $\mathcal{D}$ is drawn from this truncated solution, with input ξ(t) and output V ₀(t), depicted in Figure 4. Time-domain simulations are handled for any input signal ξ(t) by choosing appropriate discrete time approximation of the operator $\frac{\mathrm{d}}{\mathrm{d}t}$.

Figure 4 Continuous time realisation structure of $\mathcal{D}$ truncated at N = 3, composed of sums and products of the diaphragm displacement and its time derivatives.

3.3 Evaluation of the Doppler distortion effect

The distortion due to the Doppler effect is evaluated through the computation of the (truncated) Eulerian source (22) and the on-axis pressure field (1) at z = 1 m. First, a piston velocity signal composed of a low-frequency tone at f ₀ and a high-frequency tone at f _a is chosen, expressed as

$${\xi }^{\mathrm{\prime }}\left(t\right)=A\left(\sin \left(2\pi f_{0}t\right)+\sin \left(2\pi f_{a}t\right)\right),$$(23)where A, f ₀ and f _a correspond to a realistic high-amplitude loudspeaker motion parameters (see Tab. 1).

The amplitude spectrum of the acoustic pressure P(f) for the piston velocity signal (23) is depicted in Figure 5, where P(z = 1, f) is the discrete Fourier transform or p(z = 1, t). One should note that the level of the low-frequency f ₀ is lower than that of the high-frequency tone f _a because of the time derivation of the velocity signal in (1). Two kinds of signal distortion are observed:

• Harmonic distortion (HD) at frequency 2f ₀ = 80 Hz. The level of the harmonic is −56 dB below the fundamental f ₀. In the case of loudspeaker diaphragm motion, this effect is often considered as neglectable [6].

• Intermodulation distortion (IMD) at frequencies f _a  ± f ₀, f _a  ± 2f ₀ and f _a  ± 3f ₀. The low-frequency tone f ₀ modulates the high-frequency tone, generating sidelobes around f _a  = 1 kHz. The first-order intermodulation peaks f _a  ± f ₀ reach −35 dB below the fundamental f _a .

Figure 5 Amplitude spectrum of the acoustic pressure for a piston velocity signal composed of a low-frequency tone at $f_0=40$ Hz and a high-frequency tone at $f_a=1$ kHz. Intermodulation distortion is visible around $f_a$, together with a slight harmonic distortion at $2f_0$.

The above-mentionned intermodulation distortion effect can be seen in the time domain in Figure 6, where the signal p(z = 1, t) (solid line) is represented for a few periods 1/f _a . A time shift of the high-frequency waveform is observed, compared to the acoustic pressure without Doppler effect (dashed line). This corresponds to the time-domain version of the frequency modulation caused by the Doppler effect.

Figure 6 Time domain signal $p(z=1,t)$. A time shift of the acoustic pressure (solid line) is observed, compared to the case without Doppler (dashed line).

In order to evaluate the influence of the high-frequency tone on the Doppler distortion, a second input velocity signal is examined, composed of a low-frequency tone at f ₀ and a logarithmic chirp whose instantaneous frequency varies between f _a and f _b . The piston velocity expression reads

$${\xi }^{\mathrm{\prime }}\left(t\right)=A\left(\sin \left(2\pi f_{0}t\right)+\cos \left(2\pi {\int }_0^t\mathrm{}{f}_b{\left(\frac{f_b}{f_a}\right)}^{t/T}\mathrm{d}t\right)\right),$$(24)where f _b and T are defined in Table 1. The spectrograms of the piston velocity signal (24) and the acoustic pressure are depicted respectively in Figures 7 and 8.

Figure 7 Spectrogram of the input piston velocity signal ξ′(t), consisting of a chirp in the range [$f_a$, $f_b$] and a low-frequency tone $f_0$ (not visible here).

Figure 8 Spectrogram of the output acoustic pressure $p(z=1,t)$. Intermodulation distortion is observed, increasing with the chirp frequency $f\in [f_a,f_b]$.

Only the chirp part of the signal in the range [f _a , f _b ] is visible on these graphs. Sidelobes are still observed along the rising frequency. Second-order intermodulation products of frequency f ± f ₀ are noticed with amplitudes in the range [−35 dB, −25 dB]. Weak third order intermodulation (f ± 2f ₀) is also noted. This highlights once again the intermodulation distortion created by the low-frequency tone $f_0$. Moreover the sidelobes level clearly increase with the chirp frequency: the first sidelobe f ± f ₀ starts at −35 dB at f _a  = 1000 Hz to reach −25 dB at f _b  = 2000 Hz.

This numerical evaluation states that the Doppler effect of a vibrating piston causes mostly intermodulation distortion in the radiated pressure, in agreement with previous studies. This distortion increases with the ratio f _max/f _min, where f _min and f _max are respectively the minimum and maximum frequency of excitation. For a realistic loudspeaker diaphragm motion, it is shown that this effect can potentially produce sound corruption at −25 dB below the unaltered pressure signal.

3.4 Nonlinear wave propagation effect

Several authors [1, 10] have pointed out that the moving-boundary (Doppler) effect should be tackled together with the nonlinear propagation phenomenon, since both appear at large diaphragm displacements. This section briefly investigates this coupling by substituting the linear propagation in H2 for the inviscid Burger’s equation, describing the nonlinear propagation of progressive plane waves,

$${\mathrm{\partial }}_{t}v(z,t)+(c_0+\beta v(z,t\left)\right){\mathrm{\partial }}_{x}v(z,t)=0,$$(25)where β is a coefficient derived from the nonlinear relationship between pressure and air density. This equation can be solved by using the method of characteristics as in Section 3.1. In that case, the characteristic reads

$$K_{\xi }:\hspace{1em}(z,t)\mapsto \left(z,t+\frac{z-\xi \left(t\right)}{c_0+\beta {\xi }^{\mathrm{\prime }}\left(t\right)}\right),$$(26)where the term 1/(c ₀ + βξ′(t)) stands for the nonlinear propagation effect: the propagation speed differs from c ₀ and is higher for high-amplitude waves. The moving-boundary is still taken into account by the term z − ξ(t). The existence of $K_{\xi }^{-1}$ is limited (i) in space by the nonlinear propagation (z should be small enough so that no shockwave appears) and (ii) in amplitude by the Mach number condition as in (9).

A formal series solution for the particle velocity v(z, t) is derived by applying a perturbation method similar to Section 3.2. The three first orders are listed below:

$$v_1(z,\tau )={\xi }^{\prime }\left(\tau \right),$$(27)

$$v_2(z,\tau )=v_2^{\mathrm{D}}\left(\tau \right)+v_2^{\mathrm{C}}(z,\tau ),$$(28)

$$v_3(z,\tau )=v_3^{\mathrm{D}}\left(\tau \right)+v_3^{\mathrm{C}}(z,\tau ),$$(29)where τ = t − z/c ₀, $v_2^{\mathrm{D}}$ and $v_3^{\mathrm{D}}$ correspond to the terms solely due to the Doppler effect (defined in (22)) and the convective terms $v_2^{\mathrm{C}}$ and $v_3^{\mathrm{C}}$ are

$$v_2^{\mathrm{C}}(z,\tau )=\frac{\beta z}{c_0^2}{\xi }^{\prime }\left(\tau \right){\xi }^{″}\left(\tau \right),$$(30)

$$\begin{array}{c}{v}_3^{\text{C}}(z,\tau )=\frac{-\beta }{c_0^2}{\xi }^{\prime }\left(\tau \right){\xi }^{\prime }(\tau )\xi (\tau )c+\frac{\beta z}{c_0^3}[{\xi }^{\prime \prime \prime }\left(\tau \right){\xi }^{\prime }(\tau )\xi (\tau )\\ +{({\xi }^{\prime \prime }(\tau ))}^2\xi (\tau )+\left(1-\beta \right){\xi }^{\prime \prime }(\tau ){({\xi }^{\prime \prime }(\tau ))}^2]\\ +\frac{{\beta }^2{z}^2}{c_0^4}\left[{({\xi }^{\prime \prime }\left(\tau \right))}^2{\xi }^{\prime }(\tau )+\frac{{\xi }^{\prime \prime \prime }(\tau ){({\xi }^{\prime }(\tau ))}^2}{2}\right],\end{array}$$(31)

The contribution of the convection can be separated from the Doppler effect up to the order 3, as described in equations (27)–(29). One should note that this separation is not possible for higher order terms since both effects will be coupled. Moreover the convective terms at z = 0 yield

$$v_2^{\mathrm{C}}(r=0,t)=0,$$(32)

$$v_3^{\mathrm{C}}(r=0,t)=\frac{-\beta }{c_0^2}{\xi }^{″}\left(t\right){\xi }^{\prime }\left(t\right)\xi \left(t\right).$$(33)

It appears that orders 1 and 2 does not influence the equivalent Eulerian source $V_0\left(t\right)$ at z = 0. However the third order term $v_3^{\mathrm{C}}$ has a nonzero value at z = 0. This corresponds to the influence of nonlinear propagation in the domain Ω = [−ξ _max, ξ _max], that alternates the equivalent Eulerian source at z = 0 calculated in (22). Therefore the Doppler effect can be treated independently from the nonlinear propagation at orders 1 and 2, but higher orders require consideration of both effects coupled. In the following part, the anti-Doppler system is derived from third order expansion of V ₀(t) presented in Section 3.2, thus neglecting the term (33) introduced by the nonlinear propagation.

4 Inverse problem: anti-Doppler system

This section presents the derivation of a corrector ${\mathcal{D}}^{\star }$ that provides a diaphragm displacement ξ _c (t) compensating for the Doppler effect. This feed-forward controller relies on the design of a pre-inverse system, as described in Figure 9. First, the operator $\mathcal{D}$ is recasted into the Volterra series formalism in Section 4.1, for which pre-inversion methods exist and are easily tractable [13]. Then a realisation structure of ${\mathcal{D}}^{\star }$ is built from the calculated pre-inverse kernels in Section 4.2.

Figure 9 Description of the pre-inverse system ${\mathcal{D}}^{\star }$. The tandem system $\mathcal{D}$ $\circ {\mathcal{D}}^{\star }$ generates the identity.

4.1 Volterra series formulation

Proposition 4.1

The system $\mathcal{D}$ admits a representation into Volterra series and its transfer kernels are given by

$$D_n\left(s\right)=\{\begin{array}{ll}s& \mathrm{if} n=1,\\ {\displaystyle \sum _{k=1}^{n-1}}\mathrm{}\sum _{i\in {\mathcal{C}}_{n-1,k}}\mathrm{}{s}_n^{k+1}{D}_{i_1}^{ϵ}\left({\sigma }^1\right(s\left)\right)...D_{i_k}^{ϵ}\left({\sigma }^k\right(s\left)\right)& \mathrm{otherwise},\end{array}$$(34) where

$$D_n^{ϵ}\left(s\right)=\{\begin{array}{ll}1/c_0& \mathrm{if} n=1,\\ {\displaystyle \sum _{k=1}^{n-1}}\mathrm{}\sum _{i\in {\mathcal{C}}_{n-1,k}}\mathrm{}\frac{s_n^k}{c_0}{D}_{i_1}^{ϵ}\left({\sigma }^1\right(s\left)\right)...D_{i_k}^{ϵ}\left({\sigma }^k\right(s\left)\right)& \mathrm{otherwise},\end{array}$$(35) and $s\in {\mathbb{C}}^{1\times n}\left|\mathfrak{R}\right(s)>0$ , ${\sigma }^a=(s_{i_1+...+i_{a-1}+1},...,s_{i_1+...+i_a})$ .

The proof is given in Appendix C. The three first kernels calculated from this proposition are listed below:

$$D_1\left(s_1\right)=s_1,$$(36)

$$D_2(s_1,s_2)=\frac{s_1^2}{c_0},$$(37)

$$D_3(s_1,s_2,s_3)=\frac{1}{c_0^2}\left(\frac{s_1^3}{2}+s_1{s}_2^2\right).$$(38)

4.2 Third-order corrector calculation

The purpose of the anti-Doppler filter is to reshape the diaphragm displacement waveform ξ(t) into another displacement ξ _c (t) in order to reach the following target Eulerian velocity field,

$$V_0^{\star }\left(t\right)={\xi }^{\prime }\left(t\right).$$(39)

The pre-inverse system presented in Figure 9 is calculated so that the tandem system $\mathcal{D}\circ {\mathcal{D}}^{\star }$ of input $V_0^{\star }\left(t\right)={\xi }^{\mathrm{\prime }}\left(t\right)$ must be equal to the identity. This property is translated into the Volterra series formalism by

$$1=D_1^{\star }\left(s_1\right)D_1\left(s_1\right).$$(40)

Moreover, higher order terms of the series composition must vanish, yielding one equation per order. Solving recursively for n = 1, 2, 3 yields

$$D_1^{\star }\left(s_1\right)=1/s_1,$$(41)

$$D_2^{\star }(s_1,s_2)=-\frac{s_1}{c_0{s}_2(s_1+s_2)},$$(42)

$$D_3^{\star }(s_1,s_2,s_3)=\frac{\frac{s_1{s}_2}{s_3(s_2+s_3)}+\frac{s_1^2}{2s_2{s}_3}}{c_0^2(s_1+s_2+s_3)}.$$(43)

The realisation structure of (41)–(43) is presented in Figure 10. A prior integration of V ₀(t) = ξ′(t) is assumed, so that ξ(t) is the input of ${\mathcal{D}}^{\star }$. This correction algorithm is also decomposed into orders of homogeneous degree, where for instance ξ _c1(t) = ξ(t) is the contribution from the first order.

Figure 10 Continuous time realisation structure of ${\mathcal{D}}^{\star }$ truncated at N = 3, built from the transfer kernels $D_1^{\star }$ $D_2^{\star }$, $D_3^{\star }$. This structure is composed of products and sums of differentiators and integrators.

Note that the second-order term

$$\begin{array}{ll}{\xi }_{c2}\left(t\right)& =-\frac{1}{c_0}{\int }_0^t\mathrm{}{\xi }^{\mathrm{″}}\left(t\right)\xi \left(t\right)\mathrm{d}t\\ & =\frac{1}{c_0}\left({\int }_0^t\mathrm{}\left[{\xi }^{\mathrm{\prime }}\right(t){]}^2\mathrm{d}t-\xi (t\left){\xi }^{\mathrm{\prime }}\right(t)\right)\end{array}$$(44)is not centered on zero and constantly increases with time due to the integral of [ξ′(t)]². A similar divergence can be observed for the third-order term. This problem, illustrated in Figure 11, can be interpreted as follows: the particle velocity generated by the moving piston is slightly more compressed for positive values and relaxed for negative values (see top-right picture in Fig. 11), leading to a signal assymetry and thus a non-zero (negative) mean value. Therefore the required diaphragm displacement ξ _c (t) to perfectly compensate for the Doppler effect must constantly move toward the listener (see down-left picture in Fig. 11). To get around this problem, another corrector ${\mathcal{D}}_{\mathbf{c}}^{\star }$ is proposed, that substitutes the perfect integrator $1/s$ in the Laplace domain for 1/(s + 2πf _c ), where f _c is below the lowest frequency of interest. The impact of this substitution is evaluated on the orders 1 and 2 by connecting in tandem ${\mathcal{D}}_{\mathit{c}}^{\star }$ and $\mathcal{D}$ in equations (45) and (46).

• The identity is retrieved ($1$ if n = 1 and $0$ for n = 2) by setting a = 0.

• The composition $D_1\circ D_{c,1}^{\star }$ behaves like a gain of $1$ for frequencies greater than f _c .

• The composition $D_2\circ D_{c,2}^{\star }$ is equivalent to the product of two filters s/(s + 2πf _c ), connected in tandem with an integrator 1/(s + 2πf _c ) and a gain of −2πf _c . The result of this composition is very close to $0$ for frequencies greater than f _c .

$$D_1(s_1)\circ D_{c,1}^{\star }=\frac{s_1}{(s_1+a)}$$(45)

$$D_2(s_1,s_2)\circ D_{c,2}^{\star }=\frac{-s_2{s}_{1}a}{(s_1+a)(s_2+a)(s_1+s_2+a)}$$(46)

Figure 11 Illustration of the (slow) divergence of the displacement signal for a centered target Eulerian field. The centered displacement signal $\xi \left(t\right)$ (top-left) generates an Eulerian field $V_0\left(t\right)$ (top-right) with a negative mean value. Therefore, a centered target Eulerian field (down-right) requires a corrected piston velocity ${\xi }_c^{\mathrm{\prime }}\left(t\right)$ thas has a positive mean value (down-middle) and the corresponding displacement ${\xi }_c\left(t\right)$ constantly increases with time (down-left). For the sake of illustration, calculations have been carried out with exaggerated displacement amplitude of 1 m at 10 Hz.

Similar observations can be made by evaluating the third order. It can be concluded that ${\mathcal{D}}_{\mathbf{c}}^{\star }$ tends to equal ${\mathcal{D}}^{\star }$ for frequencies higher than $f_c$, therefore restoring the exact correction for sufficiently high frequencies.

4.3 Numerical evaluation of the corrector

The evaluation of the corrector is carried out with the same piston velocity signal as in the previous Section 3, composed of a low-frequency tone and a rising chirp. The cutoff frequency of ${\mathcal{D}}_{\mathbf{c}}^{\star }$ is set to $f_c=0.1$ Hz and the numerical integration is achieved by the bilinear transform of $1/(s+2\pi f_c)$. A complete block diagram of the correction evaluation is presented in Figure 12: the targeted Eulerian field is sent to the corrector ${\mathcal{D}}_{\mathbf{c}}^{\star }$ that generates a corrected piston displacement ${\xi }_c\left(t\right)$. Then the Eulerian source $V_0\left(t\right)$ is calculated through $\mathcal{D}$ and the acoustic pressure is evaluated at z = 1 m by (1).

Figure 12 Block diagram of the corrector evaluation, with the target Eulerian field as input $V_0^{\star }\left(t\right)={\xi }^{\mathrm{\prime }}\left(t\right)$ and the acoustic pressure as output $p(z,t)$.

Spectrogram of the corrected acoustic pressure $p(z=1,t)$, presented in Figure 13, should be compared to the case without correction in Figure 8. Second order intermodulation is now below −36 dB (compared to −25 dB without correction) and third order is no longer visible.

Figure 13 Spectrogram of the acoustic pressure $p(z=1,t)$ for a corrected diaphragm displacement ${\xi }_c\left(t\right)$. Intermodulation distortion is reduced compared to Figure 8.

Finally, an intermodulation factor is evaluated at the swept sine frequency $f$ from the following formula,

$${\mathrm{IMD}}_P(f)=100\frac{\sqrt{\frac{1}{4}{\displaystyle \sum _{n=1}^2}\mathrm{}\left[P(f+n{f}_0{)}^2+P(f-n{f}_0{)}^2\right]}}{\left|P\right(f\left)\right|},$$(47)where $P\left(f\right)$ is the discrete Fourier transform of the acoustic pressure $p(z=1,t)$. This factor is calculated at various ratio $f/f_0$, with and without correction of the piston displacement.

The results, presented in Figure 14, confirms the increase of intermodulation distortion with frequency: the case without correction reaches 15% of IMD for $f/f_0=100$. For $f/f_0=50$, which corresponds to a realistic case $f_0=40$ Hz and f = 2 kHz, the IMD is about 7.5%. The correction algorithm ${\mathcal{D}}_{\mathbf{c}}^{\star }$ maintains the intermodulation level below 4%. The slight loss of efficiency of ${\mathcal{D}}^{\star }$ with frequency is due to (i) the truncation at order $3$ and (ii) the phase and amplitude limitations of the implemented differentiator filter.

Figure 14 Evolution of the intermodulation distortion $\mathrm{IM}{D}_P\left(f\right)$ with $f/f_0$. The intermodulation distortion due to Doppler effect (in blue) reaches 15% at $f/f_0=100$. The correction algorithm ${\mathcal{D}}^{\star }$ (in red) maintains the distortion level below 4%.

5 Conclusion

The influence of the Doppler effect induced by a vibrating source has been investigated in terms of intermodulation distortion, based on the truncation of an exact series expansion that generates an equivalent Eulerian velocity field. The sound artefact generated by this nonlinear phenomenon increases with the distance between the highest and the lowest frequencies of the source signal.

Distortions can be expected up to 7.5% for diaphragm displacements with a wide frequency content (typically [40 Hz, 2 kHz]) and large amplitudes (a few millimeters). In this situation, truncating the series expansion at order 3 is sufficient to represent this distortion effect with a good accuracy. Also, the nonlinear propagation (due to convection) has been examined: this effect only appears from order 3. More precisely, compared to the Doppler alone, the convection effect only modifies one coefficient of the third order contribution.

In addition, an algorithm to compensate for this effect while preserving a centered version of the membrane displacement has been proposed. The solution is derived by inverting the series of the direct problem and replacing each time-integrator by a “band-limited integrator”. Numerical tests show that this solution yields a satisfying distortion reduction.

A future work is concerned with the experimental validation of the proposed Doppler system on loudspeakers. Another one will be devoted to the validation of the anti-Doppler system, through its combination with correction techniques to reject the driver electromechanical nonlinearities.

Conflict of interest

Authors declared no conflict of interests.

References

Appendix A

Property 1: $K_{\xi }$ is a diffeomorphism

Property A.1 (${\mathit{K}}_{\mathit{\xi }}$ is a ${\mathcal{C}}^n$-diffeomorphism)

Let $\xi$ be a ${\mathcal{C}}^1$ -regular function that satisfies the condition 9 on the Mach number. If function $\xi$ is also ${\mathcal{C}}^{n+1}$ -regular with $n\in \mathbb{N}$ , then

• $K_{\xi }$ is a ${\mathcal{C}}^n$ -regular diffeomorphism. Consequently the functions $K_{\xi }^{-1}$ and ${\tau }_{\xi }^d=[{\tau }_{\xi }^a{]}^{-1}$ exist and are ${\mathcal{C}}^n$ -regular,

• ${\mathbb{K}}_{\xi }^a={\mathbb{K}}_{\xi }^d={\mathbb{K}}_{\xi }$ .

Proof. (i) Bijection. By construction of the image set ${\mathbb{K}}_{\xi }^{\mathrm{a}}$ (see (6)), $K_{\xi }$ is a surjective function. Moreover $\xi$ is ${\mathcal{C}}^1$-regular, so that ${\tau }_{\xi }^{\mathrm{a}}$ and then $K_{\xi }$ are also ${\mathcal{C}}^1$-regular functions. For all $(z,t)\in {\mathbb{K}}_{\xi }^{\mathrm{d}}$, the Jacobian of $K_{\xi }$ is given by

$$J_{K_{\xi }}(z,t)=\left(\begin{array}{ll}1& 0\\ 1/c_0& 1-{\xi }^{\mathrm{\prime }}\left(t\right)/c_0\end{array}\right),$$in which $1-{\xi }^{\mathrm{\prime }}\left(t\right)/c_0$ is strictly positive from (9), therefore $K_{\xi }$ is bijective.

(ii) Diffeomorphism. Now let $n\in \mathbb{N}$ and consider $\xi \in {\mathcal{C}}^{n+1}$-regular. Then, $K_{\xi }^{-1}$ exists and can be expressed as

$$\begin{array}{ll}{K}_{\xi }^{-1}: {\mathbb{K}}_{\xi }^{\mathrm{a}}& \mapsto {\mathbb{K}}_{\xi }^{\mathrm{d}}\\ (z,t)& \mapsto (z,{\tau }_{\xi }^{\mathrm{d}}(t-z/c_0\left)\right),\end{array}$$(A1)where ${\tau }_{\xi }^{\mathrm{d}}$ is the reciprocal of ${\tau }_{\xi }^{\mathrm{a}}$. The Jacobian of $K^{-1}$ is given by

$$\begin{array}{ll}{J}_{K^{-1}}: \hspace{1em}{\mathbb{K}}_{\xi }^{\mathrm{a}}& \to {\mathcal{M}}_{\mathrm{2,2}}\left(\mathbb{R}\right)\\ & (z,t)\mapsto \left(\begin{array}{ll}1& 0\\ \frac{-1}{c_0-{\xi }^{\mathrm{\prime }}\left({\tau }_{\xi }^{\mathrm{d}}\right(t-z/c_0\left)\right)}& \frac{1}{1-{\xi }^{\mathrm{\prime }}\left({\tau }_{\xi }^{\mathrm{d}}\right(t-z/c_0\left)\right)/c_0}\end{array}\right)\end{array}.$$(A2)

Now, we prove by induction that ${\tau }_{\xi }^{\mathrm{d}}$ is ${\mathcal{C}}^p$-regular for $1\le p\le n$.

• Case $p=1$ : ${\tau }_{\xi }^d\in {\mathcal{C}}^1$ .

• From (4), ${\tau }_{\xi }^{\mathrm{a}}$ is ${\mathcal{C}}^1$-regular, so that ${\tau }_{\xi }^{\mathrm{d}}$ exists and is continuous. Then the Jacobian $J_{K^{-1}}$ defined in (A2) is a continuous function, so that $K_{\xi }^{-1}$ and then ${\tau }_{\xi }^{\mathrm{d}}$ are ${\mathcal{C}}^1$-regular.

• Case $p\ge 2$ : If ${\tau }_{\xi }^d$ is ${\mathcal{C}}^p$ -regular, then ${\tau }_{\xi }^d$ is ${\mathcal{C}}^{p\pm 1}$ -regular.

• From (A2), the Jacobian $J_{K^{-1}}$ is ${\mathcal{C}}^p$-regular. Therefore and $K^{-1}$ ${\tau }_{\xi }^{\mathrm{d}}$ are ${\mathcal{C}}^{p+1}$-regular, which proves that $K_{\xi }$ is a diffeomorphism.

Moreover, $K_{\xi }\left(\xi \right(t),t)=\left(\xi \right(t),t)$ proves that the codomain ${\mathbb{K}}_{\xi }^{\mathrm{a}}$ is also bounded by the fixed point $\left(\xi \right(t),t)$, leading to ${\mathbb{K}}_{\xi }^{\mathrm{a}}={\mathbb{K}}_{\xi }^{\mathrm{d}}={\mathbb{K}}_{\xi }$. □

Appendix B

Derivation of the perturbation method

Injecting (19) into (18) yields

$${\displaystyle \sum _{n=0}^{\mathrm{\infty }}}\mathrm{}{ϵ}_n(t){\alpha }^n=\alpha u(t+{\displaystyle \sum _{n=0}^{\mathrm{\infty }}}\mathrm{}{ϵ}_n(t\left){\alpha }^n\right)/c_0,$$(B1)where function u can be expanded into Taylor series at point t, leading to

$${\displaystyle \sum _{n=0}^{\mathrm{\infty }}}\mathrm{}{ϵ}_n\left(t\right){\alpha }^n=\frac{1}{c_0}\alpha {\displaystyle \sum _{m=0}^{\mathrm{\infty }}}\mathrm{}\frac{u^{\left(m\right)}\left(t\right)}{m!}{\left({\displaystyle \sum _{n=0}^{\mathrm{\infty }}}\mathrm{}{ϵ}_n\left(t\right){\alpha }^n\right)}^m.$$(B2)

The power series composition at the right-hand side is simplified as follows,

$${\displaystyle \sum _{n=1}^{\mathrm{\infty }}}\mathrm{}{ϵ}_n\left(t\right){\alpha }^n=\frac{\alpha u\left(t\right)}{c_0}+{\displaystyle \sum _{q=1}^{\mathrm{\infty }}}\mathrm{}{\alpha }^{q+1}{\displaystyle \sum _{k=1}^q}\mathrm{}\frac{u^{\left(k\right)}\left(t\right)}{k!c_0}\sum _{i\in {\mathcal{C}}_{q,k}}\mathrm{}{ϵ}_{i_1}\left(t\right)...{ϵ}_{i_k}\left(t\right)$$(B3)where ${\mathcal{C}}_{n,k}=\left\{\right(i_1,...,i_k\left)\right|i_1+...+i_k=n\}$ is the set of compositions of $n$ into $k$ parts. Finally, equating (B3) for each power of $\alpha$ yields the following recursive relation, for all $n\in {\mathbb{N}}^{\mathrm{*}}$,

$${ϵ}_n\left(t\right)=\{\begin{array}{ll}\xi \left(t\right)/c_0& \mathrm{if} n=1,\\ {\displaystyle \sum _{k=1}^{n-1}}\mathrm{}\frac{{\xi }^{\left(k\right)}\left(t\right)}{k!c_0}\sum _{i\in {\mathcal{C}}_{n-1,k}}\mathrm{}{ϵ}_{i_1}\left(t\right)...{ϵ}_{i_k}\left(t\right)& \mathrm{otherwise},\end{array}$$(B4)where the amplitude $\alpha$ is set to $1$, so that $u=\xi$.

Now that (18) is solved, (17) can be written

$$V_0\left(t\right)={\xi }^{\mathrm{\prime }}\left(t+{\displaystyle \sum _{q=1}^{\mathrm{\infty }}}\mathrm{}{ϵ}_q\left(t\right)\right),$$(B5)where ξ′ can be expanded into Taylor series at point t

$$V_0\left(t\right)={\displaystyle \sum _{n=0}^{\mathrm{\infty }}}\mathrm{}\frac{{\xi }^{(n+1)}\left(t\right)}{n!}{\left({\displaystyle \sum _{q=1}^{\mathrm{\infty }}}\mathrm{}{ϵ}_q\left(t\right)\right)}^n.$$(B6)

Appendix C

Proof of Proposition 4.1

Proof. The formula (35) is proven by induction.

• Case $p=1$: $D_1^{ϵ}$.

The proof is straightforward since ${ϵ}_1\left(t\right)=\xi \left(t\right)/c_0$, so that its associated transfer kernel is $D_1^{ϵ}\left(s_1\right)=1/c_0$.

• Case $p\ge 2$: if $D_p^{ϵ}$ is true, then $H_{p\pm 1}^{ϵ}$ is true.

From its construction using the regular perturbation method, ${ϵ}_p\left(t\right)$ is a homogeneous multivariate polynomial of degree $p$, so that its associated transfer kernel only consists of order $p$. Thus, the transfer kernel of following product

$$p_{i,k}\left(t\right)={\xi }^{\left(k\right)}\left(t\right){ϵ}_{i_1}\left(t\right)...{ϵ}_{i_k}\left(t\right)$$(C1)is given by

$$P_{i,k}\left(s\right)=s_{i_1+...+i_k+1}^k{h}_{i_1}^{ϵ}\left({\sigma }^1\right(s\left)\right)...h_{i_k}^{ϵ}\left({\sigma }^k\right(s\left)\right).$$(C2)

Therefore the transfer kernel associated with

$${ϵ}_{p+1}\left(t\right)={\displaystyle \sum _{k=1}^p}\mathrm{}\sum _{i\in {\mathcal{C}}_{p,k}}\mathrm{}{p}_{i,k}\left(t\right)/c_0$$(C3)is

$$D_{p+1}^{ϵ}\left(s\right)={\displaystyle \sum _{k=1}^p}\mathrm{}\sum _{i\in {\mathcal{C}}_{p,k}}\mathrm{}{P}_{i,k}\left(s\right)/c_0,$$(C4)that concludes the proof by induction.

The transfer kernels $H_n\left(s\right)$ are derived similarly. □

All Tables

All Figures

	Figure 1 Description of the physical model composed of (a) on-axis baffled piston radiation for the evaluation of the sound pressure and (b) conversion of the piston displacement into an equivalent Eulerian velocity source $V_0\left(t\right)$.
In the text

	Figure 2 Block diagram of the system under consideration, with input $\xi \left(t\right)$ and output $p(z,t)$.
In the text

Figure 3 Illustration of the mapping from departure coordinates $\left(\xi \right(t^{\mathrm{d}}),t^{\mathrm{d}})$ to arrival coordinates $(Z,t^{\mathrm{a}})$ through the characteristic line (described by $K_{\xi }(z,t^{\mathrm{d}})$ when $z$ varies), for a sinusoidal piston motion. Note the variation in the propagation time $t^{\mathrm{a}}-t^{\mathrm{d}}=\frac{Z-\xi \left(t^{\mathrm{d}}\right)}{c_0}$.

In the text