Spherical shock waveform reconstruction by heterodyne interferometry

Hart, Carl R.; Lyons, Gregory W.; White, Michael J.

doi:10.1121/10.0024520

The indirect measurement of shock waveforms by acousto-optic sensing requires a method to reconstruct the field from the projected data. Under the assumption of spherical symmetry, one approach is to reconstruct the field by the Abel inversion integral transform. When the acousto-optic sensing modality measures the change in optical phase difference time derivative, as for a heterodyne Mach–Zehnder interferometer, e.g., a laser Doppler vibrometer, the reconstructed field is the fluctuating refractive index time derivative. A technique is derived that reconstructs the fluctuating index directly by assuming plane wave propagation local to a probe beam. With synthetic data, this approach is compared to the Abel inversion integral transform and then applied to experimental data of laser-induced shockwaves. Time waveforms are reconstructed with greater accuracy except for the tail of the waveform that maps spatially to positions near a virtual origin. Furthermore, direct reconstruction of the fluctuating index field eliminates the required time integration and results in more accurate shock waveform peak values, rise times, and positive phase duration.

I. INTRODUCTION

The visualization of shockwaves by acousto-optic sensing is an important diagnostic in compressible hydrodynamic studies.^1,2 Under special circumstances, such as cylindrical or spherical symmetry in the field, a single visualization or point measurement can provide quantitative data regarding shockwave characteristics. This follows from the projection-slice theorem, a fundamental result in the tomographic reconstruction of fields.³ Shock waveform reconstruction, by the inverse Abel transform, was demonstrated earlier for point homodyne interferometry.⁴ A reconstruction technique is proposed here for point heterodyne interferometry.

General classes of acousto-optic sensing modalities include shadowgraphy, schlieren, and interferometry.^5,6 Both shadowgraphy⁷ and schlieren⁸ generate projections of the field according to second or first spatial derivatives of the refractive index, respectively. In contrast, interferometry generates projections of the refractive index field that are measured as changes in the optical phase difference between a reference and probe beam. Thus, reconstruction of the index field proceeds from changes in optical phase difference. By the Gladstone–Dale relationship⁵ the density field can be computed from the index field, and then the pressure field from an appropriate equation of state. While interferometer noise levels are generally not of concern for shock wave measurements, low-noise interferometers can be suitable for mid- to high-frequency acoustic measurements.⁹

Little prior work exists on reconstructing shock waveforms from interferometer measurements. Notably, Smeets¹⁰ was the first to reconstruct pressure waveforms from spark-generated shockwaves. Two assumptions were made by Smeets:¹⁰ (1) the fluctuating index of refraction is constant over the in-air length of the probe beam (10 mm in length), and (2) fluctuating pressure is linearly related to fluctuating density. Time series of shock waveforms were estimated by point interferometer measurements. The plasma physics community has made a couple of studies. A refraction index profile, showing the leading shock front induced by a laser-induced plasma, was reconstructed from an interferogram.¹¹ Air density profiles, with time delays of 50–2650 ns after initiation of an in-air plasma, were reconstructed from interferograms as well.¹² Returning to research in the acoustics community, free-field measurements of spark-generated shockwaves were used to reconstruct pressure waveforms from 100 to 1000 mm in range from the origin.⁴ Similarly, measurements over solid and irregular surfaces were used to reconstruct pressure waveforms for spark-generated shockwaves at a range of 130 and 330 mm, respectively.^13,14 The fluctuating index of refraction was reconstructed from the measured temporal change in optical phase differences, giving consideration to spatial variations in the fluctuating index over the probe beam.^4,13,14

A spatial waveform is derived from time-series data by assuming the waveform changes little over the duration of the pulse and that the waveform is effectively scanned by a probe beam propagating at the ambient speed of sound toward the shockwave origin.⁴ This approximation permits use of the Abel integral transform and inverse integral transform for the processing of measurement data. Consequently, the spatial range coordinate is virtual, since the actual spatial waveform differs due to absorption, relaxation, and convective propagation effects. Considering simply convective propagation effects, i.e., higher overpressures propagate faster than lower overpressures,¹⁵ a linear time-to-space mapping has the effect of shifting in range both the shock front and back shock closer to the origin.

Synthetic waveform reconstructions, based on propagation predictions of a blast waveform¹⁶ with a generalized Burgers equation, show at most a relative error of 6% in reconstructing the peak negative pressure, 1% in the peak positive pressure, and 1.2% in the positive phase duration for a propagation distance of 100 mm.⁴ The greater relative error in the peak negative pressure prediction is to be expected since the Abel inverse integral transform weights portions of the waveform more heavily near the origin rather than farther from the origin,¹⁷ i.e., the back shock virtually shifted closer to the origin would have less reconstruction accuracy than the shock front. Despite the errors introduced by this linear time-to-space mapping, it is a practical and efficient way to reconstruct the shock pressure waveform.

Another optical system is a heterodyne Mach–Zehnder interferometer, which is typically the core of a laser Doppler vibrometer (LDV).¹⁸ Instead of measuring the change in optical phase difference, output is proportional to its time derivative, which is the velocity component parallel to the probe beam. This type of interferometer contains an optical diode (composed of a polarizing beam splitter and quarter waveplate) that permits the probe beam to travel back along its in-air path and then be combined with a modulated reference beam. The measured velocity is a combination of vibrations at a reflecting surface and a projection of the fluctuating index of refraction time derivative along the in-air probe beam path. Therefore, reconstruction of the fluctuating index by the Abel inverse transform integral requires time integration. By assuming the velocity data obeys a plane wave propagation equation, for short distances, it is possible to directly reconstruct the fluctuating index.

In Sec. II, the theory for estimating the shock pressure waveforms, by way of reconstructing the fluctuating index directly from velocity data, is described. This is followed by synthetic experiments in Sec. III. Waveform reconstruction is then applied to experimental measurements of laser-induced shockwaves in Sec. IV. Last, Sec. V provides concluding remarks.

II. ACOUSTO-OPTIC RECONSTRUCTION THEORY

The theory of reconstructing a shock waveform directly from projected data, i.e., LDV velocity, is given. A numerical approximation to the reconstruction procedure, applicable for discrete-time signals, is given utilizing two- and three-point interpolation for an inverse integral transform.

A. Inversion for fluctuating index of refraction field

Presume that a LDV probe beam is offset from the origin of a spherical shockwave as shown in Fig. 1. The shockwave origin is at point O, and the closest point of approach for the probe beam is at point A:

(x_{p}, 0, z_{p})

⁠. The probe beam is parallel to the y axis of the Cartesian coordinate system, and the offset distance is

r_{p} = \sqrt{x_{p}^{2} + z_{p}^{2}}

⁠. Positions along the probe beam have Cartesian coordinates

(x_{p}, y, z_{p})

⁠, and the distance to each position is

r = \sqrt{r_{p}^{2} + y^{2}}

⁠. As a shockwave propagates across the probe beam, the fluctuating index of refraction field causes a change in the optical path length of the probe beam. This in turn results in changes to the optical phase difference between probe and reference beams as¹⁸^,

ϕ (r_{p}, t) - ϕ (r_{p}, 0) = 2 (\frac{2 π}{λ}) \int_{- \infty}^{\infty} n (r, t) d y,

(1)

where the factor of 2 in front of

2 π / λ

accounts for a probe beam passing through the field twice, e.g., being reflected from a mirror; n is the fluctuating index of refraction field; and

ϕ

is the optical phase difference between the reference and probe beams. The shockwave does not exist for time t < 0, i.e., the optical phase difference prior to the generation of the fluctuating index field is due to ambient vibrations, thermal variations in air, and system optics. Taking the time derivative of Eq. gives

ϕ_{t} (r_{p}, t) = 2 (\frac{2 π}{λ}) \int_{- \infty}^{\infty} n_{t} (r, t) d y,

(2)

where

n_{t} (r, t) = \partial n (r, t) / \partial t

⁠, and

ϕ_{t} (r_{p}, t) = \partial ϕ (r_{p}, t) / \partial t

⁠. By definition,¹⁸^,

ϕ_{t} (r_{p}, t) = 2 (\frac{2 π}{λ}) v (r_{p}, t)

(3)

states that the time-varying optical path length, i.e., the velocity signal, is proportional to a time-varying change in the optical phase difference. While this is true for a vibrating surface that reflects a probe beam, it also holds for time variations in the optical path length caused by a fluctuating index of refraction field. Substituting Eq. into Eq. gives

v (r_{p}, t) = \int_{- \infty}^{\infty} n_{t} (r, t) d y .

(4)

Importantly, the quantity

v (r_{p}, t)

is the time-varying output signal from the LDV at range r_p, which is a projection of the fluctuating index time derivative. As opposed to a vibration velocity, this velocity signal is an apparent vibration velocity that arises from accumulated phase shifts along the probe beam, which are induced by a time-varying refractive index field. Given spherical symmetry in the fluctuating index field, the Abel forward integral transform of Eq. results in¹⁹

v (r_{p}, t) = 2 \int_{r_{p}}^{\infty} \frac{n_{t} (r, t) r d r}{\sqrt{r^{2} - r_{p}^{2}}} .

(5)

As Yuldashev et al.⁴ remark, there is no analytical solution for

n_{t} (r, t)

⁠, since it varies in space and time.

FIG. 1.

View large Download slide

(Color online) The coordinate system for the expanding shock front and probe beam. The shock origin is located at point O, and the closest point of approach for the probe beam is at point A: $(x_{p}, 0, z_{p})$ ⁠. The distance from O to A is r_p. The radial coordinate constrained along the probe beam is $(x_{p}, y, z_{p})$ with distance r.

Consider a linear space-time transformation that maps the fluctuating index field to space at an instant time. Assuming the propagation speed is equal to the ambient speed of sound c₀ results in⁴

t = t' - (r - r_{p}) / c_{0},

(6)

where

t'

is a travel duration to the probe beam position, say,

t' = r_{p} / c_{0}

⁠. Applying this space-time transformation to Eq. gives

\tilde{v} (r_{p}) = 2 \int_{r_{p}}^{\infty} \frac{{\tilde{n}}_{t} (r) r d r}{\sqrt{r^{2} - r_{p}^{2}}},

(7)

where

{\tilde{n}}_{t} (r) = n_{t} (r, t = t' - (r - r_{p}) / c_{0})

and

\tilde{v} (r_{p}) = v (r_{p}, t = t' - (r - r_{p}) / c_{0})

⁠.

By the Abel inversion integral transform,¹⁹

{\tilde{n}}_{t} (r) = - \frac{1}{π} \int_{r}^{\infty} \frac{{\tilde{v}}_{r_{p}} (r_{p}) d r_{p}}{\sqrt{r_{p}^{2} - r^{2}}},

(8)

where

{\tilde{v}}_{r_{p}} (r_{p}) = \partial \tilde{v} (r_{p}) / \partial r_{p}

⁠. Now, assume the fluctuating index waveform changes little over the duration of the pulse and obeys, local to the probe beam, a plane wave propagation equation for outgoing waves. In a similar manner, the projected field (LDV velocity signal) will obey an outgoing plane wave equation,¹⁵

{\tilde{v}}_{r_{p}} (r_{p}) = - \frac{1}{c_{0}} {\tilde{v}}_{t} (r_{p}) .

(9)

This approximation follows from Eq. ,

v_{r_{p}} (r_{p}, t) + \frac{1}{c_{0}} v_{t} (r_{p}, t) = \int_{- \infty}^{\infty} [n_{t r_{p}} (r, t) + \frac{1}{c_{0}} n_{t t} (r, t)] d y,

(10)

where, for plane wave propagation, the bracketed integrand on the right-hand side will be zero for all y. In other words, under the assumption of plane wave propagation for the fluctuating index field, Eq. shows that the projected field propagates as a plane wave, which upon application of the time-to-space mapping yields Eq. . While the plane wave approximation for the index field decreases in accuracy closer to the shockwave origin, with respect to the field geometry, as the curvature of the shock front decreases, the approximation improves. Substituting Eq. into Eq. and integrating in time yields

\tilde{n} (r) = \frac{1}{π c_{0}} \int_{r}^{\infty} \frac{\tilde{v} (r_{p}) d r_{p}}{\sqrt{r_{p}^{2} - r^{2}}} .

(11)

A time waveform for the reconstructed field is obtained from Eq. by inverting the space-time transformation,

\tilde{n} (t) = \tilde{n} (r = r_{p} - c_{0} (t - t'))

⁠, which is an approximate solution to

n (r_{p}, t)

⁠. To briefly summarize, the steps for reconstructing a shock waveform are (1) a time series is obtained for the projected field at the probe beam offset distance r_p, as in Eq. ; (2) a space-time transformation, Eq. , is applied to the projected field to place the waveform in virtual range; (3) an inversion integral transform, Eq. , obtains the fluctuating refraction field in virtual range coordinates; and (4) the space-time transformation is inverted to obtain the fluctuating refraction field in time, at the probe beam offset distance.

Equation (11) bears some resemblance to a typical Abel inversion integral transform, such as Eq. (8), but is different in two aspects. First, within the integrand is the projected field (LDV velocity signal) rather than its spatial derivative. Second, the sign is positive. Directly operating on the projected field is advantageous for noisy data, since numerical differentiation in space will typically increase the error in a reconstruction. The (non-dimensional) right-hand side of Eq. (11) contains a weighted ratio of velocities, which is consistent with the units expected for the index field. If spherical wave propagation is presumed for the velocity signal, then an additional term (velocity signal divided by probe beam offset distance) must be introduced on the right-hand side of Eq. (9). Reconstruction for the fluctuating index must then integrate this additional term in time, which Eq. (11) avoids.

B. Inverse transform by two-point and three-point interpolation

To solve Eq. with discrete-time LDV velocity data, a numerical approach is required. By subdividing the integral over each discrete interval, interpolating the velocity data, and analytically solving each subdivided integral, the inverse transform can be represented by a deconvolution operation,¹⁷

\tilde{n} (r_{i}) = \frac{1}{c_{0}} \sum_{j = 0}^{N - 1} D_{i j} \tilde{v} (r_{j}),

(12)

where D_ij contains weighting coefficients appropriate to the approximation made for the inverse transform integral,

r_{i} = i Δ r, r_{j} = j Δ r, i = 0, 1, \dots, N - 1

⁠, N is the sample count of the LDV velocity-time signal, and

Δ r = c_{0} Δ t

is the spatial resolution. Notably, Eq. does not explicitly depend on the spatial resolution

Δ r

⁠, but rather is a weighted sum of the velocity data divided by the ambient speed of sound. While Eq. can be computed by numerical quadrature, given discrete-time data and some appropriate handling of the singularity when r_p = r, the advantage of a deconvolution approach via Eq. is computational speed (approximately 100 times faster in comparison to adaptive quadrature) without significant losses in accuracy. In what follows, linear and quadratic interpolations are considered, while it is possible to consider cubic interpolation as well.⁸ It can be expected that higher orders of interpolation will lead to higher orders of numerical convergence.

First, consider a linear interpolation for the LDV velocity signal. Subdividing the inverse transform integral, Eq. , into discrete portions, analytically solving for each integral, and rearranging summation terms yields for the deconvolution operator

D_{i j}^{(2)} = {\begin{matrix} 0, & j < i, \\ J_{i j} (0) - J_{i j} (1), & j = i, \\ J_{i j} (0) - J_{i j} (1) + J_{i, j - 1} (1), & j > i . \end{matrix}

(13)

The term

J_{i j} (0)

is¹⁷

J_{i j} (0) = {\begin{matrix} 0, & j < i, \\ \frac{2}{π}, & j = i = 0, \\ \frac{1}{π} ln [\frac{1 + j + \sqrt{{(j + 1)}^{2} - i^{2}}}{j + \sqrt{j^{2} - i^{2}}}], & j \geq i, \end{matrix}

(14)

and

J_{i j} (1)

is defined as

J_{i j} (1) = {\begin{matrix} 0, & j < i, \\ \frac{1}{π}, & j = i = 0, \\ \frac{1}{π} [\sqrt{{(j + 1)}^{2} - i^{2}} - \sqrt{j^{2} - i^{2}}] - j J_{i j} (0), & j \geq i . \end{matrix}

(15)

Detailed derivation is given in Appendix A. The order of numerical convergence for

D_{i j}^{(2)}

⁠, independent of the time-to-space mapping, is

O (Δ r)

⁠.

Now, consider a quadratic interpolation for the LDV velocity signal. Subdividing the inverse transform integral, Eq. , into discrete portions, analytically solving for each integral, and rearranging summation terms yields for the deconvolution operator

D_{i j}^{(3)} = {\begin{matrix} 0, & j < i - 1, \\ \frac{1}{8} I_{i, j + 1} (2) - \frac{1}{4} I_{i, j + 1} (1), & j = i - 1, \\ \frac{1}{8} I_{i, j + 1} (2) - \frac{1}{4} [I_{i, j + 1} (1) + I_{i j} (2)] + I_{i j} (0), & j = i, \\ \frac{1}{8} [I_{i, j + 1} (2) + I_{i, j - 1} (2)] + \frac{1}{4} [I_{i, j - 1} (1) - I_{i, j + 1} (1) - I_{i j} (2)] + I_{i j} (0), & j \geq i + 1, \\ \frac{1}{8} I_{i, j + 1} (2) + \frac{1}{4} [I_{i, j - 1} (2) - I_{i, j + 1} (1) - I_{i j} (2)] + I_{i j} (0), & i = 0, j = 1, \end{matrix}

(16)

where in the last condition it is presumed that

v (r_{- 1}) = v (r_{1})

⁠, following Dasch.¹⁷ The terms

I_{i j} (0)

and

I_{i j} (1)

are¹⁷

I_{i j} (0) = {\begin{matrix} 0, & j = i = 0 or j < i, \\ \frac{1}{π} ln [\frac{1 + 2 j + \sqrt{{(2 j + 1)}^{2} - {(2 i)}^{2}}}{2 j}], & j = i \neq 0, \\ \frac{1}{π} ln [\frac{1 + 2 j + \sqrt{{(2 j + 1)}^{2} - {(2 i)}^{2}}}{- 1 + 2 j + \sqrt{{(2 j - 1)}^{2} - {(2 i)}^{2}}}], & j > i, \end{matrix}

(17)

I_{i j} (1) = {\begin{matrix} 0, & j < i, \\ \frac{1}{π} \sqrt{{(2 j + 1)}^{2} - {(2 i)}^{2}} - 2 j I_{i j} (0), & j = i, \\ \frac{1}{π} [\sqrt{{(2 j + 1)}^{2} - {(2 i)}^{2}} - \sqrt{{(2 j - 1)}^{2} - {(2 i)}^{2}}] - 2 j I_{i j} (0), & j > i . \end{matrix}

(18)

The term

I_{i j} (2)

is

I_{i j} (2) = {\begin{matrix} 0, & j < i, \\ \frac{1}{π} (\frac{1}{2} - 3 j) \sqrt{{(2 j + 1)}^{2} - {(2 i)}^{2}} + 6 j^{2} I_{i j} (0), & j = i, \\ \frac{1}{π} [(\frac{1}{2} - 3 j) \sqrt{{(2 j + 1)}^{2} - {(2 i)}^{2}} + (\frac{1}{2} + 3 j) \sqrt{{(2 j - 1)}^{2} - {(2 i)}^{2}}] + (4 j^{2} + 2 i^{2}) I_{i j} (0), & j > i . \end{matrix}

(19)

Detailed derivation is given in Appendix B. The order of numerical convergence for

D_{i j}^{(3)}

⁠, independent of the time-to-space mapping, is

O ({(Δ r)}^{2})

⁠.

III. SYNTHETIC EXPERIMENT

To test the validity of Eq. (11), synthetic experiments were conducted numerically.²⁰ The objective is to compare the use of the two- and three-point deconvolution operators, Eqs. (12), (13), and (16), which reconstructs the refraction index field directly from the velocity signal, i.e., the projected field, against the three-point deconvolution operator by Dasch,¹⁷ which reconstructs the index time derivative field from the spatial derivative of the velocity field, as in Eq. (8).

A. Friedlander pulse

A synthetic signal was constructed according to the form of a Friedlander pulse.¹⁶ The peak amplitude is assumed to decay nonlinearly, while the waveform propagates according to the speed of sound,

p (r, t) = \frac{p_{0}}{2} {(\frac{r_{0}}{r})}^{α} (1 - \frac{t - r / c_{0}}{τ}) (1 - tanh [\frac{2 (r / c_{0} - t)}{τ_{sh}}]) e^{- (t - r / c_{0}) / τ},

(20)

where p₀ is the peak pressure at a range r₀, α is the decay rate, τ is the positive phase duration, and

τ_{sh}

is the shock rise time. The shock front corresponds to that for a dissipative weak shock, such that the shock rise time is inversely proportional to the peak pressure.¹⁵ In assuming a constant propagation speed, the differences between each deconvolution approach can be evaluated independently of errors introduced by the linear time-to-space conversion, Eq. , which are associated with an amplitude dependence for the sound speed. However, other errors can arise in the time-to-space mapping due to neglecting decay in the waveform amplitude and portions of the waveform that may be mapped to negative range values.

By the isentropic gas law and the Dale–Gladstone relationship between density and index of refraction, the fluctuating index is related to acoustic pressure as

n (r, t) = ρ_{0} K [{(1 + \frac{p (r, t)}{P_{0}})}^{1 / γ} - 1],

(21)

where ρ₀ is the ambient air density, K is the Gladstone–Dale coefficient for 632.8 nm light,⁵ P₀ is the ambient air pressure, and γ is the ratio of specific heats. The time derivative of Eq. is straightforward to derive,

n_{t} (r, t) = p_{t} (r, t) \frac{ρ_{0} K}{γ P_{0}} {(1 + \frac{p (r, t)}{P_{0}})}^{1 / γ - 1},

(22)

where

p_{t} (r, t) = p (r, t) [- \frac{1}{τ - (t - r / c_{0})} - \frac{1}{τ} + \frac{2}{τ_{sh}} (1 + tanh [\frac{2 (r / c_{0} - t)}{τ_{sh}}])] .

(23)

B. Constants and signal processing parameters

For the synthetic experiments, several constants are determined according to the International Standard Atmosphere, and pulse parameters are approximately based on past experiments.²¹ Table I lists every constant specified. The choices for $t_{min}, t_{max}$ ⁠, and N ensure that the temporal resolution gives approximately 10 points over the shock rise time for probe beam offset ranges greater than or equal to 70 mm.

TABLE I.

Synthetic experiment constants for ambient air conditions, pulse properties, and signal processing parameters.

Description	Symbol	Value
Ambient temperature	T_C	15 °C
Ambient speed of sound	c₀	340 m/s
Ambient density	ρ₀	1.225 kg/m³
Ambient pressure	P₀	101 325 Pa
Gladstone–Dale constant	K	2.256 × 10^–4 m³/kg
Ratio of specific heats	γ	7/5
Pulse positive duration	τ	See Eq. (25)
Peak pressure	p₀	3000 Pa
Radius for peak pressure	r₀	0.1 m
Peak amplitude decay rate	α	1.23
Number of signal samples	N	65 536
Minimum time	$t_{min}$	$r_{p} / c_{0} - 6 τ_{sh}$
Maximum time	$t_{max}$	$2 r_{p} / c_{0} + τ$

Description	Symbol	Value
Ambient temperature	T_C	15 °C
Ambient speed of sound	c₀	340 m/s
Ambient density	ρ₀	1.225 kg/m³
Ambient pressure	P₀	101 325 Pa
Gladstone–Dale constant	K	2.256 × 10^–4 m³/kg
Ratio of specific heats	γ	7/5
Pulse positive duration	τ	See Eq. (25)
Peak pressure	p₀	3000 Pa
Radius for peak pressure	r₀	0.1 m
Peak amplitude decay rate	α	1.23
Number of signal samples	N	65 536
Minimum time	$t_{min}$	$r_{p} / c_{0} - 6 τ_{sh}$
Maximum time	$t_{max}$	$2 r_{p} / c_{0} + τ$

C. Computational steps

The pressure signal is constructed as a function of time for a radial distance corresponding to the probe beam offset distance, i.e., r = r_p. The minimum time $t_{min}$ and maximum time $t_{max}$ are selected to ensure that a zero pressure amplitude region ahead of the shock front maps to the maximum virtual range via Eq. (6) and the tail of the waveform maps to the range origin when $t' = r_{p} / c_{0} + τ$ ⁠, respectively. The choice of $t'$ ensures the waveform zero crossing is positioned at the probe beam offset distance in virtual range. The fluctuating refractive index and index time derivative are computed according to Eqs. (21) and (22). For each discrete time, the forward Abel transform, Eq. (5), is applied to the index time derivative Eq. (22). This is accomplished in two parts: (1) numerical integration by adaptive quadrature applied over the limits of $r_{p} + Δ r$ to the end of the waveform $c_{0} (t + 6 τ_{sh})$ and (2) analytical solution of the integral over the limits of r_p to $r_{p} + Δ r$ ⁠. An analytical solution is required to handle the singularity at r = r_p, which is found by linearly interpolating the index time derivative over the spatial interval, similar to how $D_{i j}^{(2)}$ was derived. The result is the projected velocity field in time at the given probe beam offset distance. The time waveform is converted to a virtual spatial waveform by Eq. (6). No smoothing or truncation is applied to the waveform at the virtual range origin. The deconvolution operators, Eqs. (13) and (16), are applied to the spatial waveform to obtain the index spatial waveform, which is then converted back to a time waveform. Similarly, the deconvolution operator from Dasch,¹⁷ denoted as $D_{i j}^{(3 A)}$ in the figures, is applied to the spatial velocity waveform to obtain a spatial waveform for the index time derivative. Upon obtaining the time series for the index time derivative, time integration yields the index. Each time waveform is plotted against the original index time waveform.

Figure 2 shows the computational steps to reconstruct the fluctuating index of refraction at a specific probe beam position. In Fig. 2(a), the projected data (velocity) time series is shown. A time-to-space transformation brings the waveform into virtual coordinates, as shown in Fig. 2(b). The spatial waveform is deconvolved by Eqs. (12) and (16) to obtain the fluctuating index of refraction, as shown in Fig. 2(c). Finally, the time-to-space transformation is reversed to obtain a time waveform for the fluctuating index of refraction, as in Fig. 2(d).

FIG. 2.

$(Color online) Steps to reconstruct the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 70 mm from the origin. Steps are arranged in a clockwise fashion. (a) A velocity field is synthesized by numerical evaluation of Eq. (5). (b) A time-to-space conversion constructs a spatial waveform from the velocity time series. (c) Deconvolution by a three-point inverse integral transform reconstructs the spatial index field. (d) The time-to-space conversion is reversed to obtain a time series for the reconstructed index field. The original signal, shown in orange, is nearly overlapped by the reconstructed waveform shown in dashed blue.$

View large Download slide

(Color online) Steps to reconstruct the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 70 mm from the origin. Steps are arranged in a clockwise fashion. (a) A velocity field is synthesized by numerical evaluation of Eq. (5). (b) A time-to-space conversion constructs a spatial waveform from the velocity time series. (c) Deconvolution by a three-point inverse integral transform reconstructs the spatial index field. (d) The time-to-space conversion is reversed to obtain a time series for the reconstructed index field. The original signal, shown in orange, is nearly overlapped by the reconstructed waveform shown in dashed blue.

D. Results

Figure 3 shows the results of deconvolving the synthetic signal at a probe beam position 70 mm from the origin. Near the rise and peak of the Friedlander pulse, the deconvolved results are very similar. Very little difference exists between the two- and three-point deconvolution results, as evidenced by their overlapping curves in both the waveforms and errors. Differences exist in the decaying portions of the reconstructed waveforms between the three-point Abel deconvolution and either the two- or three-point deconvolutions. The relative error (⁠ $| [n (t_{i}) - \tilde{n} (t_{i})] / n (t_{i}) |$ ⁠) and absolute error (⁠ $| n (t_{i}) - \tilde{n} (t_{i}) |$ ⁠) are computed on a sample-by-sample basis, with errors associated with $D_{i j}^{(2)}$ and $D_{i j}^{(3)}$ visually overlapping. The relative error goes to infinity at the original waveform zero crossing. In a similar manner, the relative error goes to zero at the crossings between the reconstructed and original waveforms. Later in time, the relative error grows since the original waveform tends to zero. In addition to these extrema, the relative error is roughly 10% over the majority of the waveform. In general, the absolute error decreases over time. The exception is that the absolute error increases toward the tail end of the reconstructed waveforms, when using either $D_{i j}^{(2)}$ or $D_{i j}^{(3)}$ ⁠, as shown in Fig. 4. This is a consequence of assuming the projected field propagates as a plane wave; however, absolute errors do not exceed values near the shock front. Notable is the smaller magnitude in relative and absolute errors in the reconstructed waveform tails, except very near the end, which corresponds to points that map near the virtual range origin.

FIG. 3.

$(Color online) Reconstruction of the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 70 mm from the origin. Shown in (a) are the original and reconstructed waveforms, (b) the relative error, and (c) the absolute error for the reconstructed waveforms.$

View large Download slide

(Color online) Reconstruction of the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 70 mm from the origin. Shown in (a) are the original and reconstructed waveforms, (b) the relative error, and (c) the absolute error for the reconstructed waveforms.

FIG. 4.

$(Color online) Reconstruction of the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 70 mm from the origin, highlighting the tail of the waveform. Shown are (a) the original and reconstructed waveforms, (b) the relative error, and (c) the absolute error for the reconstructed waveforms.$

View large Download slide

(Color online) Reconstruction of the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 70 mm from the origin, highlighting the tail of the waveform. Shown are (a) the original and reconstructed waveforms, (b) the relative error, and (c) the absolute error for the reconstructed waveforms.

A summary of reconstruction errors, at ranges of 70–260 mm, are shown in Table II. The NRMSE is defined as

NRMSE = \frac{\sqrt{(1 / N) \sum_{j = 1}^{N} | n (t_{j}) - \tilde{n} (t_{j}) |^{2}}}{| n_{+} - n_{-} |},

(24)

where

n_{+}

is the positive peak fluctuating index, and

n_{-}

is the negative peak fluctuating index. At the reported numerical precision, errors for the two-point deconvolution method are identical to the three-point deconvolution method. Despite differences in numerical convergence for these two deconvolution methods, it is likely that errors in the time-to-space conversion dominate. In comparison to the relative error for the positive peak fluctuating index, the relative error for the negative peak index is greater for all reconstruction methods. This is to be expected since the deconvolution operators place greater weight on portions of the waveform nearer the spatial origin, which is the tail of the waveform in the time domain. However, the choice for the time shift

t' = r_{p} / c_{0} + τ

balances the errors in the front and back shocks.⁴ In comparing the errors across all three reconstruction methods, it is evident that both the two- and three-point deconvolution methods provide the most accurate estimates for negative peak index fluctuation and positive phase duration. Relative errors for positive peak index fluctuation are comparable for all three methods. This suggests that at the ranges considered, the portion of error due to the time-to-space conversion is much greater in relation to the error introduced by the plane wave approximation. The proportion of error due to the plane wave approximation is expected to increase closer to the origin. Waveform shapes are more accurately reconstructed by a two- or three-point deconvolution in comparison to the three-point Abel deconvolution as quantified by the NRMSE. As another point of comparison, the errors reported by Yuldashev et al.⁴ for a blast waveform reconstructed at a range of 100 mm are 1%, 6%, and 1.2% for positive peak pressure, negative peak pressure, and positive phase duration, respectively. The errors in reconstructing a Friedlander waveform using either a two- or three-point deconvolution at the same range are comparable.

TABLE II.

Synthetic experiment reconstruction errors for fluctuating index positive peak $n_{+}$ ⁠, negative peak $n_{-}$ ⁠, positive phase duration τ, and normalized root mean square error (NRMSE). Two deconvolution methods are compared: (1) three-point deconvolution $D_{i j}^{(3)}$ and (2) three-point Abel deconvolution $D_{i j}^{(3 A)}$ ⁠.

Method	Range (mm)	Relative error (%)			NRMSE
Method	Range (mm)	$n_{+}$	$n_{-}$	τ	NRMSE
$D_{i j}^{(3)}$	70	3.5	8.9	2.5	0.007
	80	3.2	8.1	2.3	0.006
	100	2.7	6.8	1.9	0.005
	130	2.2	5.6	1.6	0.003
	180	1.7	4.3	1.3	0.002
	260	1.3	3.2	1.0	0.001
$D_{i j}^{(3 A)}$	70	3.5	14.7	5.6	0.010
	80	3.2	13.5	5.2	0.009
	100	2.7	11.5	4.4	0.007
	130	2.2	9.5	3.7	0.005
	180	1.7	7.5	2.9	0.003
	260	1.3	5.6	2.2	0.002

Method	Range (mm)	Relative error (%)			NRMSE
Method	Range (mm)	$n_{+}$	$n_{-}$	τ	NRMSE
$D_{i j}^{(3)}$	70	3.5	8.9	2.5	0.007
	80	3.2	8.1	2.3	0.006
	100	2.7	6.8	1.9	0.005
	130	2.2	5.6	1.6	0.003
	180	1.7	4.3	1.3	0.002
	260	1.3	3.2	1.0	0.001
$D_{i j}^{(3 A)}$	70	3.5	14.7	5.6	0.010
	80	3.2	13.5	5.2	0.009
	100	2.7	11.5	4.4	0.007
	130	2.2	9.5	3.7	0.005
	180	1.7	7.5	2.9	0.003
	260	1.3	5.6	2.2	0.002

Figure 5 shows the results of deconvolving a synthetic signal at a probe beam position 10 mm from the origin and placing the shock front at the probe beam position, i.e., $t' = r_{p} / c_{0}$ ⁠. This example illustrates the presence of additional errors that arise when portions of a waveform are mapped to negative range values in the time-to-space conversion or when the waveform amplitude does not decay to zero at the virtual range origin. Apparent is the truncation of the waveform at approximately 60 μs, which corresponds to the virtual origin in spatial coordinates. Any finite-duration waveform that is partially mapped to negative ranges is to be truncated at the origin prior to being mapped back to time (the deconvolution operations only apply to positive range values), and as a consequence, portions of the reconstruction will be in error if the waveform amplitude does not decay to zero at the virtual origin. In particular, the plane wave assumption for the velocity data amplifies the relative error at the end of the reconstructed waveform, relative to deconvolution by the Abel inverse integral transform. However, the initial portion of the Friedlander pulse is reproduced fairly well, such that the positive peak fluctuating index and positive phase duration are determined with a relative error of 0.002%–0.01% and 9.1%–18.9%, respectively. Use of the three-point Abel deconvolution gives the lowest relative error for amplitude, while use of either the two- or three-point deconvolution gives the lowest relative error for positive phase duration in this case. To reconstruct an entire finite-duration waveform with reasonable accuracy, the waveform must map to positive range values upon application of the time-to-space transformation. Conservatively, assume a finite-duration pulse propagates at the speed of sound with a pulse duration T. For the end of the temporal waveform to map to positive ranges, the probe beam range must be greater than or equal to the product of the pulse duration and the speed of sound, $r_{p} \geq c_{0} T$ ⁠.

FIG. 5.

$(Color online) Reconstruction of the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 10 mm from the origin.$

View large Download slide

(Color online) Reconstruction of the fluctuating index of refraction for a Friedlander pulse arriving at a LDV probe beam position 10 mm from the origin.

IV. LASER-INDUCED SHOCKWAVE EXPERIMENT

Measurements of laser-induced shockwaves by Hart and Lyons²¹ are reanalyzed using the three-point inverse integral transform, Eqs. (12) and (16). Reconstruction of the acoustic pressure time series is shown, along with statistics for the peak positive pressure, peak negative pressure, positive phase duration, and shock rise time.

Laser-induced shockwaves were generated by a combination of a pulsed Nd:YAG laser [Quantel Q-Smart 850 (Lumibird, Lannion, France); 1064 nm, 5.5 ns, 850 mJ], a Galilean beam expander, and a 100 mm focal length plano-convex lens. A discrete number of ranges (10–200 mm, in 10 mm steps) from the shockwave center, along the direction of the laser beam, were sampled by a heterodyne Mach–Zehnder interferometer. The interferometer is part of a LDV system [NLV-2500-2, Polytec (Irvine, CA), 3 MHz, 100 (mm/s)/V] that measures time series of vibrational velocity. A digitizer card [M2p.5943-x4, Spectrum (Grosshansdorf, Germany)] sampled analog signals from the LDV at a sample rate of 80 MHz. Further details of the experiment are provided in Hart and Lyons.²¹

A. Signal processing

The signal processing follows essentially the same steps as described for the synthetic experiments in Sec. III. The ambient speed of sound is determined for moist air,¹⁵ and moist air density is computed from the measured ambient pressure, temperature, and relative humidity. A linear trend in the raw velocity data is removed, which effectively subtracts low-frequency vibrations in system optics (reflecting mirror and LDV head) and thermal perturbations. Furthermore, a nine-point moving average filter is applied to the velocity signal so as to reduce random noise and preserve a sharp step response in the signal.²²

B. Shock waveform reconstructions and time-domain statistics

Figure 6 shows averaged reconstructed pressure waveforms measured at ranges of 30–200 mm in 10 mm increments. At each range, 200 waveforms were measured, shifted in time to synchronize their peak pressures (onto the average of the unshifted arrival times), and then averaged. For the 30 mm range, reconstruction error in the tail of the waveform is evident near in time to the 60 mm range waveform zero crossing. As shown in Fig. 4, such error is to be expected. Data collected at 10 and 20 mm ranges sometimes exceeded the measurement range corresponding to the sensitivity setting used with the LDV. Therefore, those waveforms are omitted. The shape of each waveform more closely matches a blast waveform¹⁶ rather than a Friedlander pulse. Furthermore, the shape of the shock front is indicative of spatial averaging due to the finite beam width of the LDV beam (2.5–3 mm). The influence of the digitizer bandwidth and nine-point moving average filter would be minimal since a step response with a duration as low as 0.113 μs [ $9 / (80 MHz$ ⁠)] would be preserved, while the bandwidth of the LDV permits a time resolution of 0.333 μs.

FIG. 6.

View large Download slide

(Color online) Averaged shock pressure waveforms measured at ranges of 30–200 mm. Range labels shown for every other waveform.

Figure 7 shows the median peak positive $p_{+}$ and negative $p_{-}$ pressure as a function of range. Regression for the positive peak pressure, over the ranges of 50–200 mm, gives a range decay rate of 1.30, i.e., $p_{+} (r) = p_{+} (r_{0}) {(r / r_{0})}^{- α}$ ⁠, where $r_{0} = 100$ mm, $p_{+} (r_{0}) = 3098$ Pa, and $α = 1.30$ ⁠. As a point of comparison, Yuldashev et al.⁴ obtain $α = 1.2$ for spark-generated shockwaves measured over ranges of 100–1000 mm. Differences in both range and peak pressure amplitude likely contribute to the different estimated decay rates.⁴ Peak pressure at ranges of 10 and 20 mm fall below the extrapolated region of the regression line. As noted earlier, this can be attributed partly to clipping in some velocity measurements. Microphone measurements at ranges of 160–200 mm show comparable peak pressures and waveforms.²¹

FIG. 7.

View large Download slide

(Color online) Median peak positive $p_{+}$ and negative $p_{-}$ pressure as a function of range. Regression for peak positive pressure is over 50–200 mm ranges.

Positive phase duration and shock rise time medians are shown in Fig. 8. Regression for the positive phase duration takes the form²³^,

τ (r) = a_{0} r^{a_{1} - a_{2} ln r},

(25)

where r is the range in meters, τ is the positive phase duration in seconds, and by regression

a_{0} = 2.48 \times 10^{- 5}, a_{1} = 1.27 \times 10^{- 1}

⁠, and

a_{2} = 2.00 \times 10^{- 2}

⁠. Although it can be expected that the squared positive phase duration would scale linearly according to the logarithm of range beyond 200 mm,²⁴ it is not observed in the data due to shorter measurement ranges. The shock rise time median over all ranges is 2.85 μs. Given the peak pressure at ranges of 100 and 200 mm, theoretically, the shock rise time would be 0.05 and 0.11 μs, respectively.¹⁵ The measured rise time corresponds to an effective experimental bandwidth of 351 kHz and is a clear indicator that spatial averaging by the LDV probe beam limits the time resolution of these measurements. As an additional consequence, the measured peak pressure is less than what would be measured by a system with a greater effective bandwidth. As Smeets¹⁰ noted, the width of the probe beam dictates the time resolution. The probe beam is approximately 5–6 mm in diameter (1/e² width) at the output of the LDV head. Since the origin of the shock front was approximately equidistant between the LDV head and reflecting mirror, it can be expected that the beam width traversed by the initial portions of the shock front was anywhere from 2.5 to 3 mm. This is close to the diameter of a 1/8 in. microphone (3.175 mm), which has a time constant of 2.5–2.9 μs,⁴ similar to the measured median shock rise time.

FIG. 8.

View large Download slide

(Color online) Median positive phase duration τ and shock rise time $τ_{sh}$ as a function of range.

V. CONCLUSIONS

The reconstruction of shock waveforms is an essential element in quantifying shockwave characteristics from interferometer measurements. For a heterodyne Mach–Zehnder interferometer, the measurand is a projection of the fluctuating index of refraction field time derivative, which is proportional to the time derivative for the change in optical phase difference. In practice, a LDV can measure an apparent velocity signal induced by a passing shockwave or an acoustic pulse with a sufficiently sensitive instrument, which can be interpreted as a projection of the fluctuating index time derivative.

One approach to reconstruct the fluctuating index field is to apply an Abel inverse transform integral on a spherically symmetric spatial waveform. When a velocity signal, corresponding to a shockwave, is processed in this manner, the fluctuating index time derivative is reconstructed. Once converted to a time series, this must be integrated in time to obtain the fluctuating index, which can then be used to estimate the fluctuating density and pressure. However, by assuming the projected data, local to the probe beam, obeys a plane wave propagation equation for short distances, it is possible to reconstruct the fluctuating index directly as in Eq. (11), an integral inverse transform. This form of reconstruction is advantageous in several respects: (1) it avoids amplifying noise in real projected data, and (2) it eliminates errors introduced by time integration.

Numerical approximations to Eq. (11) by two- and three-point interpolation were derived and compared to a three-point Abel integral inverse transform. Relative error in the two- and three-point integral inverse transforms, for sufficiently far ranges, is lower than in the three-point Abel integral inverse transform. This is a consequence of eliminating reconstruction error by time integration.

Shock pressure waveforms were reconstructed from measurements of laser-induced shockwaves. Peak pressures, positive phase duration, and shock rise times were statistically quantified. The decay rate of the positive peak pressure is indicative of an intermediate strength shockwave, and the rise times indicated limited time resolution due to the finite width of the probe beam. A potential avenue for future research is accounting for probe beam spatial averaging by an appropriate form of deconvolution.

SUPPLEMENTARY MATERIAL

See the supplementary material for data tables enumerating first, second (median), and third quartiles for measured quantities plotted in Figs. 7 and 8.

ACKNOWLEDGMENTS

We express our thanks to Cody M. Best [U.S. Army Engineer Research Development Center, Cold Regions Research and Engineering Laboratory (ERDC-CRREL)] for assisting in the experiment. This study was conducted for the Office of the Assistant Secretary of the Army for Acquisition, Logistics, and Technology [ASA(ALT)] under Research, Development, Test, and Evaluation Program Element 601102A/AB2/02, Fundamental Adaptive Protection and Projection Research. Permission to publish was granted by the Director, U.S. Army Engineer Research Development Center, Cold Regions Research and Engineering Laboratory.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: TWO-POINT INTEGRAL INVERSE TRANSFORM

A two-point inverse integral transform of Eq. follows a similar approach to that taken by Dasch¹⁷ in deriving a two-point Abel inverse integral transform. First, the transform is subdivided into discrete portions,¹⁷^,

\tilde{n} (r_{i}) = \frac{1}{π c_{0}} \sum_{j = i}^{N - 1} \int_{0}^{Δ r} \frac{\tilde{v} (r_{j} + δ)}{\sqrt{{(r_{j} + δ)}^{2} - r_{i}^{2}}} d δ .

(A1)

Let the velocity signal be approximated by a two-point, linear interpolation,

\tilde{v} (r_{j} + δ) = \frac{δ}{Δ r} [\tilde{v} (r_{j + 1}) - \tilde{v} (r_{j})] + \tilde{v} (r_{j}) + O ({(Δ r)}^{2}),

(A2)

which can be readily derived from the Lagrange interpolation formula. Taking the derivative of Eq. with respect to δ gives a finite-difference approximation used by Dasch.¹⁷ Substituting Eq. into Eq. and making a change of variable

u = δ / Δ r

gives for the integral portion

J_{i j} (m) = \frac{1}{π} \int_{0}^{1} \frac{u^{m}}{\sqrt{{(j + u)}^{2} - i^{2}}} d u,

(A3)

for m = 0, 1. Solution for m = 0 is provided by Eq. . By Eq. , the integral is singular for m = 0, and

j = i = 0

⁠; however, as Dasch¹⁷ notes, the projected data are expanded quadratically across two data points, giving a result of

2 / π

⁠. For m = 1, Eq. is solved analytically [Ref. 25, Eqs. (2.264.2) and (2.261)] and given by Eq. . Substituting Eqs. and into Eq. results in

\tilde{n} (r_{i}) = \frac{1}{c_{0}} \sum_{j = i}^{N - 1} {J_{i j} (1) [\tilde{v} (r_{j + 1}) - \tilde{v} (r_{j})] + J_{i j} (0) \tilde{v} (r_{j})} .

(A4)

By rearranging the terms in the sum, this can be cast in the form of Eq. ,

\tilde{n} (r_{i}) = \frac{1}{c_{0}} \sum_{j = 0}^{N - 1} D_{i j}^{(2)} \tilde{v} (r_{j}),

(A5)

where the deconvolution operator is defined in Eq. .

APPENDIX B: THREE-POINT INTEGRAL INVERSE TRANSFORM

A three-point inverse integral transform of Eq. follows a similar approach to that taken by Dasch¹⁷ in deriving a three-point Abel inverse integral transform. First, the transform is subdivided into discrete portions,¹⁷^,

\tilde{n} (r_{i}) = \frac{1}{π c_{0}} \sum_{j = i}^{N - 1} \int_{\begin{matrix} 0, j = i \\ - Δ r / 2, j > i \end{matrix}}^{Δ r / 2} \frac{\tilde{v} (r_{j} + δ)}{\sqrt{{(r_{j} + δ)}^{2} - r_{i}^{2}}} d δ,

(B1)

where the lower limit of the integral is interpreted as δ = 0 if j = i or

δ = - Δ r / 2

if j > i. Let the velocity signal be approximated by a three-point, quadratic interpolation,

\tilde{v} (r_{j} + δ) = \frac{δ^{2}}{2 {(Δ r)}^{2}} [\tilde{v} (r_{j - 1}) - 2 \tilde{v} (r_{j}) + \tilde{v} (r_{j + 1})] - \frac{δ}{2 Δ r} [\tilde{v} (r_{j - 1}) - \tilde{v} (r_{j + 1})] + \tilde{v} (r_{j}) + O ({(Δ r)}^{3}),

(B2)

which, like Eq. , is readily derived from the Lagrange interpolation formula. Taking the derivative of Eq. with respect to δ gives a finite-difference approximation used by Dasch.¹⁷ Substituting Eq. into Eq. and making a change of variable

u = δ / (Δ r / 2)

gives, for the integral portion,

I_{i j} (m) = \frac{1}{π} \int_{\begin{matrix} 0, j = i \\ - 1, j > i \end{matrix}}^{1} \frac{u^{m}}{\sqrt{{(2 j + u)}^{2} - {(2 i)}^{2}}} d u,

(B3)

for m = 0, 1, 2, which is simply two times the definition given by Dasch.¹⁷ Solution for m = 0, 1 is given by Eqs. and , which corrects the sign error in

I_{i j} (1)

[Ref. 26, Appendix 2]. For m = 2, Eq. is solved analytically [Ref. 25, Eqs. (2.264.3) and (2.261)] and given by Eq. . Substituting Eqs. into Eq. results in

\tilde{n} (r_{i}) = \frac{1}{c_{0}} \sum_{j = i}^{N - 1} {I_{i j} (2) [\frac{\tilde{v} (r_{j - 1}) - 2 \tilde{v} (r_{j}) + \tilde{v} (r_{j + 1})}{8}] - I_{i j} (1) [\frac{\tilde{v} (r_{j - 1}) - \tilde{v} (r_{j + 1})}{4}] + I_{i j} (0) \tilde{v} (r_{j})} .

(B4)

By rearranging the terms in the sum, this can be cast in the form of Eq. ,

\tilde{n} (r_{i}) = \frac{1}{c_{0}} \sum_{j = 0}^{N - 1} D_{i j}^{(3)} \tilde{v} (r_{j}),

(B5)

where the deconvolution operator is defined in Eq. .

REFERENCES

1.

H. M.

Glaz

,

P.

Colella

,

I. I.

Glass

, and

R. L.

Deschambault

, “

A numerical study of oblique shock-wave reflections with experimental comparisons

,”

Proc. R. Soc. A: Math. Phys. Eng. Sci.

398

(

1814

),

117

–

140

(

1985

).

https://doi.org/10.1098/rspa.1985.0028

Google Scholar

Crossref

2.

H.

Kleine

,

H.

Olivier

,

K.

Tsuji

,

K.

Etoh

,

K.

Takehara

, and

T. G.

Etoh

, “

Time-resolved Mach–Zehnder interferometry of shock waves

,” in

28th International Symposium on Shock Waves

, edited by

K.

Kontis

(

Springer

,

Berlin

,

2012

), pp.

577

–

583

.

Google Scholar

Crossref

3.

R. M.

Mersereau

and

A. V.

Oppenheim

, “

Digital reconstruction of multidimensional signals from their projections

,”

Proc. IEEE

62

(

10

),

1319

–

1338

(

1974

).

https://doi.org/10.1109/PROC.1974.9625

Google Scholar

Crossref

4.

P.

Yuldashev

,

M.

Karzova

,

V.

Khokhlova

,

S.

Ollivier

, and

P.

Blanc-Benon

, “

Mach–Zehnder interferometry method for acoustic shock wave measurements in air and broadband calibration of microphones

,”

J. Acoust. Soc. Am.

137

(

6

),

3314

–

3324

(

2015

).

https://doi.org/10.1121/1.4921549

Google Scholar

Crossref

PubMed

5.

W.

Merzkirch

,

Flow Visualization

(

Academic

,

New York

,

1974

), pp.

64

–

70

, 102–126.

Google Scholar

6.

G. S.

Settles

,

Schlieren and Shadowgraph Techniques: Visualizing Phenomena in Transparent Media

(

Springer

,

Berlin

,

2001

), pp.

25

–

38

.

Google Scholar

Crossref

7.

P.

Yuldashev

,

S.

Ollivier

,

M.

Averiyanov

,

O.

Sapozhnikov

,

V.

Khokhlova

, and

P.

Blanc-Benon

, “

Nonlinear propagation of spark-generated N-waves in air: Modeling and measurements using acoustical and optical methods

,”

J. Acoust. Soc. Am.

128

(

6

),

3321

–

3333

(

2010

).

https://doi.org/10.1121/1.3505106

Google Scholar

Crossref

PubMed

8.

M. M.

Karzova

,

P. V.

Yuldashev

,

V. A.

Khokhlova

,

S.

Ollivier

,

E.

Salze

, and

P.

Blanc-Benon

, “

Characterization of spark-generated N-waves in air using an optical schlieren method

,”

J. Acoust. Soc. Am.

137

(

6

),

3244

–

3252

(

2015

).

https://doi.org/10.1121/1.4921026

Google Scholar

Crossref

PubMed

9.

K.

Ishikawa

,

Y.

Shiraki

,

T.

Moriya

,

A.

Ishizawa

,

K.

Hitachi

, and

K.

Oguri

, “

Low-noise optical measurement of sound using midfringe locked interferometer with differential detection

,”

J. Acoust. Soc. Am.

150

(

2

),

1514

–

1523

(

2021

).

https://doi.org/10.1121/10.0005939

Google Scholar

Crossref

PubMed

10.

G.

Smeets

, “

Laser interference microphone for ultrasonics and nonlinear acoustics

,”

J. Acoust. Soc. Am.

61

(

3

),

872

–

875

(

1977

).

https://doi.org/10.1121/1.381339

Google Scholar

Crossref

11.

H.

Zhang

,

J.

Lu

,

Z.

Shen

, and

X.

Ni

, “

Investigation of 1.06 μm laser induced plasma in air using optical interferometry

,”

Opt. Commun.

282

(

9

),

1720

–

1723

(

2009

).

https://doi.org/10.1016/j.optcom.2009.01.036

Google Scholar

Crossref

12.

S. Q.

Cao

,

M. G.

Su

,

Z. H.

Jiao

,

Q.

Min

,

D. X.

Sun

,

P. P.

Ma

,

K. P.

Wang

, and

C. Z.

Dong

, “

Dynamics and density distribution of laser-produced plasma using optical interferometry

,”

Phys. Plasmas

25

(

6

),

063302

(

2018

).

https://doi.org/10.1063/1.5028414

Google Scholar

Crossref

13.

M. M.

Karzova

,

T.

Lechat

,

S.

Ollivier

,

D.

Dragna

,

P. V.

Yuldashev

,

V. A.

Khokhlova

, and

P.

Blanc-Benon

, “

Irregular reflection of spark-generated shock pulses from a rigid surface: Mach–Zehnder interferometry measurements in air

,”

J. Acoust. Soc. Am.

145

(

1

),

26

–

35

(

2019

).

https://doi.org/10.1121/1.5084266

Google Scholar

Crossref

PubMed

14.

M. M.

Karzova

,

T.

Lechat

,

S.

Ollivier

,

D.

Dragna

,

P. V.

Yuldashev

,

V. A.

Khokhlova

, and

P.

Blanc-Benon

, “

Effect of surface roughness on nonlinear reflection of weak shock waves

,”

J. Acoust. Soc. Am.

146

(

5

),

EL438

–

EL443

(

2019

).

https://doi.org/10.1121/1.5133737

Google Scholar

Crossref

PubMed

15.

A. D.

Pierce

,

Acoustics: An Introduction to Its Physical Principles and Applications

, 3rd ed. (

Springer

,

Cham, Switzerland

,

2019

), pp.

21

, 32, 653, 676–677.

Google Scholar

Crossref

16.

J. W.

Reed

, “

Atmospheric attenuation of explosion waves

,”

J. Acoust. Soc. Am.

61

(

1

),

39

–

47

(

1977

).

https://doi.org/10.1121/1.381266

Google Scholar

Crossref

17.

C. J.

Dasch

, “

One-dimensional tomography: A comparison of Abel, onion-peeling, and filtered backprojection methods

,”

Appl. Opt.

31

(

8

),

1146

–

1152

(

1992

).

https://doi.org/10.1364/AO.31.001146

Google Scholar

Crossref

PubMed

18.

A.

Torras-Rosell

,

S.

Barrera-Figueroa

, and

F.

Jacobsen

, “

Sound field reconstruction using acousto-optic tomography

,”

J. Acoust. Soc. Am.

131

(

5

),

3786

–

3793

(

2012

).

https://doi.org/10.1121/1.3695394

Google Scholar

Crossref

PubMed

19.

R.

Bracewell

,

The Fourier Transform and Its Applications

(

McGraw-Hill

,

New York

,

1965

), pp.

262

–

266

.

Google Scholar

20.

C. R.

Hart

, “

Synthetic reconstruction of shock waveforms

,” https://hdl.handle.net/11681/47127 (Last viewed June 5, 2023) (

2023

).

Google Scholar

21.

C. R.

Hart

and

G. W.

Lyons

, “

Interferometric measurements of laser-induced shockwaves in air

,”

Proc. Mtgs. Acoust.

48

(

1

),

045003

(

2022

).

https://doi.org/10.1121/2.0001647

Google Scholar

Crossref

22.

S. W.

Smith

,

The Scientist and Engineer's Guide to Digital Signal Processing

, 2nd ed. (

California Technical

,

San Diego, CA

,

1999

), p.

277

.

Google Scholar

23.

R. D.

Ford

,

D. J.

Saunders

, and

G.

Kerry

, “

The acoustic pressure waveform from small unconfined charges of plastic explosive

,”

J. Acoust. Soc. Am.

94

(

1

),

408

–

417

(

1993

).

https://doi.org/10.1121/1.407052

Google Scholar

Crossref

24.

W. M.

Wright

, “

Propagation in air of N waves produced by sparks

,”

J. Acoust. Soc. Am.

73

(

6

),

1948

–

1955

(

1983

).

https://doi.org/10.1121/1.389585

Google Scholar

Crossref

25.

I. S.

Gradshteyn

and

I. M.

Ryzhik

,

Table of Integrals, Series, and Products

, 8th ed. (

Elsevier

,

Waltham, MA

,

2015

).

Google Scholar

26.

K. M.

Martin

, “

Acoustic modification of sooting combustion

,” Ph.D. thesis,

University of Texas at Austin

,

Austin, TX

,

1992

.

Google Scholar

2024

U.S. Government

Spherical shock waveform reconstruction by heterodyne interferometry

I. INTRODUCTION

II. ACOUSTO-OPTIC RECONSTRUCTION THEORY

A. Inversion for fluctuating index of refraction field

B. Inverse transform by two-point and three-point interpolation

III. SYNTHETIC EXPERIMENT

A. Friedlander pulse

B. Constants and signal processing parameters

C. Computational steps

D. Results

IV. LASER-INDUCED SHOCKWAVE EXPERIMENT

A. Signal processing

B. Shock waveform reconstructions and time-domain statistics

V. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

APPENDIX A: TWO-POINT INTEGRAL INVERSE TRANSFORM

APPENDIX B: THREE-POINT INTEGRAL INVERSE TRANSFORM

REFERENCES

Supplementary Material

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Spherical shock waveform reconstruction by heterodyne interferometry

I. INTRODUCTION

II. ACOUSTO-OPTIC RECONSTRUCTION THEORY

A. Inversion for fluctuating index of refraction field

B. Inverse transform by two-point and three-point interpolation

III. SYNTHETIC EXPERIMENT

A. Friedlander pulse

B. Constants and signal processing parameters

C. Computational steps

D. Results

IV. LASER-INDUCED SHOCKWAVE EXPERIMENT

A. Signal processing

B. Shock waveform reconstructions and time-domain statistics

V. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

APPENDIX A: TWO-POINT INTEGRAL INVERSE TRANSFORM

APPENDIX B: THREE-POINT INTEGRAL INVERSE TRANSFORM

REFERENCES

Supplementary Material

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only