Correction of the Doppler distortion generated by a vibrating baffled piston

The Doppler effect is a phenomenon inherent to source motion, which introduces a variable propagation time between the source and a listening point. In the case of a vibrating piston, this is responsible for distortion of the radiated sound pressure. This moving-boundary phenomenon is part of the nonlinear effects involved in loudspeaker radiation. The present paper investigates the significance of this distortion, usually considered as neglectible, and addresses its correction. First, the direct problem is solved by: (a) converting 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. A series solution of (a) is derived and timedomain simulations of (b) are built from the truncated series combined with a baffled piston radiation model. Simulations show that Doppler distortion can be significant for realistic loudspeaker diaphragm motion with a wide spectral content. Second, the inverse (anti-Doppler) problem is examined, that is, the derivation of a piston displacement that generates a targeted Eulerian velocity field. The corrected piston velocity solution proves to be an uncentered signal, leading to a diverging displacement. In order to remove this practical problem, a centered approximation is preferred, based on modified inverse Volterra kernels. The anti-Doppler algorithm is reliable in the audio range.


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][2][3][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.

Problem statement
Consider a baffled circular piston of radius R 0 , localised at position z = 0 (see Fig. 1a). For an Eulerian velocity field excitation V 0 (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Þ where q 0 is the air density and c 0 is the speed of sound. The core of this study focuses on the (nonlinear) mapping between the rigid piston displacement and the equivalent Eulerian velocity field V 0 (t), as described in Figure 1b. This issue is addressed under the following assumptions: H1 (geometry): The piston displacement n(t) is small compared to its radius (|n(t)| < n max ( R 0 ) H2 (fields equivalence approximation): In the domain X = [Àn max , n max ], the particle velocity n 0 (t) at the piston position n(t) results from the conservative plane wave propagation of an Eulerian velocity field V 0 (t) at z = 0 H3 (waveform): the piston displacement n is a smooth function and no shockwave propagates in X 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 0 vðz; tÞ ¼ V 0 ðt À z=c 0 Þ: ð2Þ and satisfies a moving-boundary condition at the piston position that reads Moreover, condition H3 implies that the Eulerian velocity field V 0 must be a regular function.
Let us denote D the operator that converts the piston displacement into its equivalent Eulerian velocity (see Fig. 2), such that V 0 ðtÞ ¼ D½nðtÞ. This study examines the following problems: P1 (Direct problem): Is-it possible to derive a solver for D that satisfies (2)-(3)?
P2 (Inverse problem): What is the corrected displacement n c (t) to provide to D so that its output is the target V 0 ðtÞ ¼ n 0 ðtÞ?
The direct problem P1 is tackled in Section 3 and the inverse problem P2 is examined in Section 4.

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

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 0 under a necessary and sufficient condition on C 1 -regular functions n. 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.
where the domain K d n and codomain K a n are K d n ¼ fðz; tÞ 2 R 2 s:t: z ! nðtÞg; K a n ¼ fK n ðz; tÞ for ðz; tÞ 2 K d n g: The characteristic lines K n 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 n . One should note that 1. the propagation is conservative so that the particle velocity is constant on the characteristic line described by K n (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 This can be written using the characteristic equation (6), It is clear from these observations that the particle velocity v at any space-time coordinates ðZ; t a Þ 2 K a n 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 À1 n , the reciprocal of K n . This point is examined in the property detailed in Appendix A, which yields a condition on the Mach number of the piston velocity, Finally, the solution of P1 is given in Theorem 3.2. It provides the regular waveform V 0 , solution of the propagation equation (2) and the boundary condition (3), under the Mach number condition (9).
Theorem 3.2 (Regular solution). Let n be a 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 C n -regular function where s d n is the reciprocal of s a n .
Proof. From the property detailed in Appendix A, K n is a C n -diffeormorphism, so that the reciprocal K À1 n exists and is C n -regular. Moreover, since (n(t), t) is a fixed point of K n . The reciprocal function K À1 n can be developped where the second component of (12) reads Finally, therefore the regular solution (10) verifies the boundary condition (3). Moreover, n 0 s d n is C n -regular because n 0 and s d n are C n -regular, which concludes the proof. h This theorem can be summarized in equations (15) and (16), with the expression of the Eulerian source where s d n ðtÞ maps the current time t at z = 0 to the departure time of the wave at the piston position, satisfying the implicit equation The velocity source V 0 (t) cannot be directly computed from the piston displacement n(t) because of the implicit Figure 3. Illustration of the mapping from departure coordinates ðnðt d Þ; t d Þ to arrival coordinates ðZ ; t a Þ through the characteristic line (described by K n ðz; t d Þ when z varies), for a sinusoidal piston motion. Note the variation in the propagation equation (16). The derivation of a direct solver for (15) and (16) is addressed in the next part.

Regular perturbation method
A direct solver is established in this subsection, assuming a regular function n 2 C 1 and jn 0 ðtÞj < c 0 , so that s d n 2 C 1 from Theorem 3.2. First, function s d n is substituted for the C 1 -regular function : t7 !s d n ðtÞ À t. Then equations (15) and (16) can be reformulated Equation (17) highlights that (t) acts as a perturbation of the standard (linear) solution V 0 (t) = n'(t). Now the implicit equation (18) is solved using the regular perturbation method. The piston displacement function n is marked with an auxiliary amplitude a > 0, so that we define n(t) = au(t). Then we seek a power series solution of (t) of the form Calculations, developed in Appendix B, yield the exact series solution of the form in which each term v D n ðtÞ are defined as follows, v D n ðtÞ ¼ n 0 ðtÞ if n ¼ 1; P nÀ1 k¼1 n ðkþ1Þ ðtÞ k! P i2C nÀ1;k i 1 ðtÞ::: i k ðtÞ otherwise; where n (i) stands for the ith derivative of n and the i is defined recursively in Appendix B.
In the sequel, computation of V 0 (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 The term n = 1 corresponds to the linear solution without taking into account the Doppler effect: the Eulerian field V 0 (t) equals the piston velocity. Considering the term n = 2, the product n(t)n 00 (t) will be of high amplitude if the piston velocity signal is composed of low frequency components (amplified by n) together with high frequencies (amplified by n 00 ), therefore generating intermodulation distortion in V 0 (t).
A realisation structure of D is drawn from this truncated solution, with input n(t) and output V 0 (t), depicted in Figure 4. Time-domain simulations are handled for any input signal n(t) by choosing appropriate discrete time approximation of the operator d dt .

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 0 and a high-frequency tone at f a is chosen, expressed as where A, f 0 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 lowfrequency f 0 is lower than that of the high-frequency tone   (1). Two kinds of signal distortion are observed: Harmonic distortion (HD) at frequency 2f 0 = 80 Hz. The level of the harmonic is À56 dB below the fundamental f 0 . In the case of loudspeaker diaphragm motion, this effect is often considered as neglectable [6].
Intermodulation distortion (IMD) at frequencies f a ± f 0 , f a ± 2f 0 and f a ± 3f 0 . The low-frequency tone f 0 modulates the high-frequency tone, generating sidelobes around f a = 1 kHz. The first-order intermodulation peaks f a ± f 0 reach À35 dB below the fundamental f a .
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.
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 0 and a logarithmic chirp whose instantaneous frequency varies between f a and f b . The piston velocity expression reads 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.
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 0 are noticed with amplitudes in the range [À35 dB, À25 dB]. Weak third order intermodulation (f ± 2f 0 ) 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 0 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.

Nonlinear wave propagation effect
Several authors [1,10] have pointed out that the movingboundary (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, o t vðz; tÞ þ ðc 0 þ bvðz; tÞÞo x vðz; tÞ ¼ 0; where b 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 where the term 1/(c 0 + bn 0 (t)) stands for the nonlinear propagation effect: the propagation speed differs from c 0 and is higher for high-amplitude waves. The movingboundary is still taken into account by the term z À n(t).
The existence of K À1 n 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; sÞ ¼ n 0 ðsÞ; ð27Þ where s = t À z/c 0 , v D 2 and v D 3 correspond to the terms solely due to the Doppler effect (defined in (22)) and the convective terms v C 2 and v C 3 are 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 n 00 ðtÞn 0 ðtÞnðtÞ: ð33Þ It appears that orders 1 and 2 does not influence the equivalent Eulerian source V 0 ðtÞ at z = 0. However the third order term v C 3 has a nonzero value at z = 0. This corresponds to the influence of nonlinear propagation in the domain X = [Àn max , n 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 0 (t) presented in Section 3.2, thus neglecting the term (33) introduced by the nonlinear propagation.

Inverse problem: anti-Doppler system
This section presents the derivation of a corrector D H that provides a diaphragm displacement n 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 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 D H is built from the calculated pre-inverse kernels in Section 4.2.
The proof is given in Appendix C. The three first kernels calculated from this proposition are listed below:

Third-order corrector calculation
The purpose of the anti-Doppler filter is to reshape the diaphragm displacement waveform n(t) into another displacement n c (t) in order to reach the following target Eulerian velocity field, The pre-inverse system presented in Figure 9 is calculated so that the tandem system D D H of input V H 0 ðtÞ ¼ n 0 ðtÞ must be equal to the identity. This property is translated into the Volterra series formalism by Moreover, higher order terms of the series composition must vanish, yielding one equation per order. Solving recursively for n = 1, 2, 3 yields The realisation structure of (41)-(43) is presented in Figure 10. A prior integration of V 0 (t) = n 0 (t) is assumed, so that n(t) is the input of D H . This correction algorithm is also decomposed into orders of homogeneous degree, where for instance n c1 (t) = n(t) is the contribution from the first order. Note that the second-order term is not centered on zero and constantly increases with time due to the integral of [n 0 (t)] 2 . 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 n 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 D H c is proposed, that substitutes the perfect integrator 1=s in the Laplace domain for 1/(s + 2pf 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 D H c and 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 D H c;1 behaves like a gain of 1 for frequencies greater than f c . The composition D 2 D H c;2 is equivalent to the product of two filters s/(s + 2pf c ), connected in tandem with an integrator 1/(s + 2pf c ) and a gain of À2pf c . The result of this composition is very close to 0 for frequencies greater than f c .
Similar observations can be made by evaluating the third order. It can be concluded that D H c tends to equal D H for frequencies higher than f c , therefore restoring the exact correction for sufficiently high frequencies.

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 D H c is set to f c ¼ 0:1 Hz and the numerical integration is achieved by the bilinear transform of 1=ðs þ 2pf c Þ. A complete block diagram of the correction evaluation is presented in Figure 12: the targeted Eulerian field is sent to the corrector D H c that generates a corrected piston displacement n c ðtÞ. Then the Eulerian source V 0 ðtÞ is calculated through D and the acoustic pressure is evaluated at z = 1 m by (1).
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.
Finally, an intermodulation factor is evaluated at the swept sine frequency f from the following formula, where Pðf Þ 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 D H c maintains the intermodulation level below 4%. The slight loss of efficiency of D H with frequency is due to (i) the truncation at order 3 and (ii) the phase and amplitude limitations of the implemented differentiator filter.

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.

K À1
n : K a n 7 ! K d n ðz; tÞ 7 ! ðz; s d n ðt À z=c 0 ÞÞ; where s d n is the reciprocal of s a n . The Jacobian of K À1 is given by J K À1 : K a n ! M 2;2 ðRÞ ðz; tÞ 7 ! 1 0 À1 c 0 À n 0 ðs d n ðt À z=c 0 ÞÞ 1 1 À n 0 ðs d n ðt À z=c 0 ÞÞ=c 0 0 @ 1 A : Now, we prove by induction that s d n is C p -regular for 1 p n.
Case p ¼ 1: s d n 2 C 1 . From (4), s a n is C 1 -regular, so that s d n exists and is continuous. Then the Jacobian J K À1 defined in (A2) is a continuous function, so that K À1 n and then s d n are C 1 -regular. Case p ! 2: If s d n is C p -regular, then s d n is C pAE1 -regular. From (A2), the Jacobian J K À1 is C p -regular. Therefore and K À1 s d n are C pþ1 -regular, which proves that K n is a diffeomorphism.
Moreover, K n ðnðtÞ; tÞ ¼ ðnðtÞ; tÞ proves that the codomain K a n is also bounded by the fixed point ðnðtÞ; tÞ, leading to K a n ¼ K d n ¼ K n . h