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.^{3} Shock waveform reconstruction, by the inverse Abel transform, was demonstrated earlier for point homodyne interferometry.^{4} 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^{7} and schlieren^{8} 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^{5} 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.^{9}

Little prior work exists on reconstructing shock waveforms from interferometer measurements. Notably, Smeets^{10} was the first to reconstruct pressure waveforms from spark-generated shockwaves. Two assumptions were made by Smeets:^{10} (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.^{11} Air density profiles, with time delays of 50–2650 ns after initiation of an in-air plasma, were reconstructed from interferograms as well.^{12} 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.^{4} 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.^{4} 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,^{15} 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^{16} 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.^{4} 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,^{17} 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).^{18} 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

*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 = 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 = 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

^{18}

^{,}

*n*is the fluctuating index of refraction field; and $\varphi $ 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. (1) gives

^{18}

^{,}

*r*, 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. (4) results in

_{p}^{19}

*et al.*

^{4}remark, there is no analytical solution for $ n t ( r , t )$, since it varies in space and time.

*c*

_{0}results in

^{4}

^{19}

^{15}

*y*. In other words, under the assumption of plane wave propagation for the fluctuating index field, Eq. (10) shows that the projected field propagates as a plane wave, which upon application of the time-to-space mapping yields Eq. (9). 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. (9) into Eq. (8) and integrating in time yields

*r*, as in Eq. (4); (2) a space-time transformation, Eq. (6), is applied to the projected field to place the waveform in virtual range; (3) an inversion integral transform, Eq. (11), 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.

_{p}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

^{17}

*D*contains weighting coefficients appropriate to the approximation made for the inverse transform integral, $ r i = i \Delta r , \u2009 r j = j \Delta r , \u2009 i = 0 , 1 , \u2026 , N \u2212 1$,

_{ij}*N*is the sample count of the LDV velocity-time signal, and $ \Delta r = c 0 \Delta t$ is the spatial resolution. Notably, Eq. (12) does not explicitly depend on the spatial resolution $ \Delta r$, but rather is a weighted sum of the velocity data divided by the ambient speed of sound. While Eq. (11) 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. (12) 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.

^{8}It can be expected that higher orders of interpolation will lead to higher orders of numerical convergence.

^{17}

^{17}The terms $ I i j ( 0 )$ and $ I i j ( 1 )$ are

^{17}

## III. SYNTHETIC EXPERIMENT

To test the validity of Eq. (11), synthetic experiments were conducted numerically.^{20} 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,^{17} which reconstructs the index time derivative field from the spatial derivative of the velocity field, as in Eq. (8).

### A. Friedlander pulse

^{16}The peak amplitude is assumed to decay nonlinearly, while the waveform propagates according to the speed of sound,

*p*

_{0}is the peak pressure at a range

*r*

_{0},

*α*is the decay rate,

*τ*is the positive phase duration, and $ \tau 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.

^{15}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. (6), 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.

*ρ*

_{0}is the ambient air density,

*K*is the Gladstone–Dale coefficient for 632.8 nm light,

^{5}

*P*

_{0}is the ambient air pressure, and

*γ*is the ratio of specific heats. The time derivative of Eq. (21) is straightforward to derive,

### 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.^{21} Table I lists every constant specified. The choices for $ t min , \u2009 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.

Description . | Symbol . | Value . |
---|---|---|

Ambient temperature | T _{C} | 15 °C |

Ambient speed of sound | c_{0} | 340 m/s |

Ambient density | ρ_{0} | 1.225 kg/m^{3} |

Ambient pressure | P_{0} | 101 325 Pa |

Gladstone–Dale constant | K | 2.256 × 10^{–4} m^{3}/kg |

Ratio of specific heats | γ | 7/5 |

Pulse positive duration | τ | See Eq. (25) |

Peak pressure | p_{0} | 3000 Pa |

Radius for peak pressure | r_{0} | 0.1 m |

Peak amplitude decay rate | α | 1.23 |

Number of signal samples | N | 65 536 |

Minimum time | $ t min$ | $ r p / c 0 \u2212 6 \tau sh$ |

Maximum time | $ t max$ | $ 2 r p / c 0 + \tau $ |

Description . | Symbol . | Value . |
---|---|---|

Ambient temperature | T _{C} | 15 °C |

Ambient speed of sound | c_{0} | 340 m/s |

Ambient density | ρ_{0} | 1.225 kg/m^{3} |

Ambient pressure | P_{0} | 101 325 Pa |

Gladstone–Dale constant | K | 2.256 × 10^{–4} m^{3}/kg |

Ratio of specific heats | γ | 7/5 |

Pulse positive duration | τ | See Eq. (25) |

Peak pressure | p_{0} | 3000 Pa |

Radius for peak pressure | r_{0} | 0.1 m |

Peak amplitude decay rate | α | 1.23 |

Number of signal samples | N | 65 536 |

Minimum time | $ t min$ | $ r p / c 0 \u2212 6 \tau sh$ |

Maximum time | $ t max$ | $ 2 r p / c 0 + \tau $ |

### 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 \u2032 = r p / c 0 + \tau $, respectively. The choice of $ t \u2032$ 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 + \Delta r$ to the end of the waveform $ c 0 ( t + 6 \tau sh )$ and (2) analytical solution of the integral over the limits of

*r*to $ r p + \Delta r$. An analytical solution is required to handle the singularity at

_{p}*r*=

*r*, 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,

_{p}^{17}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).

### 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 ) \u2212 n \u0303 ( t i ) ] / n ( t i ) |$) and absolute error ( $ | n ( t i ) \u2212 n \u0303 ( 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.

^{4}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.*

^{4}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.

Method . | Range (mm) . | Relative error (%) . | NRMSE . | ||
---|---|---|---|---|---|

$ n +$ . | $ n \u2212$ . | τ
. | |||

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

$ n +$ . | $ n \u2212$ . | τ
. | |||

$ 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 \u2032 = 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 \u2265 c 0 T$.

## IV. LASER-INDUCED SHOCKWAVE EXPERIMENT

Measurements of laser-induced shockwaves by Hart and Lyons^{21} 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.^{21}

### 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,^{15} 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.^{22}

### 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^{16} 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 \u2009 MHz$)] would be preserved, while the bandwidth of the LDV permits a time resolution of 0.333 *μ*s.

Figure 7 shows the median peak positive $ p +$ and negative $ p \u2212$ 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 ) \u2212 \alpha $, where $ r 0 = 100$ mm, $ p + ( r 0 ) = 3098$ Pa, and $ \alpha = 1.30$. As a point of comparison, Yuldashev *et al.*^{4} obtain $ \alpha = 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.^{4} 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.^{21}

^{23}

^{,}

*r*is the range in meters,

*τ*is the positive phase duration in seconds, and by regression $ a 0 = 2.48 \xd7 10 \u2212 5 , \u2009 a 1 = 1.27 \xd7 10 \u2212 1$, and $ a 2 = 2.00 \xd7 10 \u2212 2$. Although it can be expected that the squared positive phase duration would scale linearly according to the logarithm of range beyond 200 mm,

^{24}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.

^{15}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

^{10}noted, the width of the probe beam dictates the time resolution. The probe beam is approximately 5–6 mm in diameter (1/e

^{2}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,

^{4}similar to the measured median shock rise time.

## 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

^{17}in deriving a two-point Abel inverse integral transform. First, the transform is subdivided into discrete portions,

^{17}

^{,}

*δ*gives a finite-difference approximation used by Dasch.

^{17}Substituting Eq. (A2) into Eq. (A1) and making a change of variable $ u = \delta / \Delta r$ gives for the integral portion

*m*= 0, 1. Solution for

*m*= 0 is provided by Eq. (14). By Eq. (A3), the integral is singular for

*m*= 0, and $ j = i = 0$; however, as Dasch

^{17}notes, the projected data are expanded quadratically across two data points, giving a result of $ 2 / \pi $. For

*m*= 1, Eq. (A3) is solved analytically [Ref. 25, Eqs. (2.264.2) and (2.261)] and given by Eq. (15). Substituting Eqs. (14) and (15) into Eq. (A1) results in

### APPENDIX B: THREE-POINT INTEGRAL INVERSE TRANSFORM

^{17}in deriving a three-point Abel inverse integral transform. First, the transform is subdivided into discrete portions,

^{17}

^{,}

*δ*= 0 if

*j*=

*i*or $ \delta = \u2212 \Delta r / 2$ if

*j*>

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

*δ*gives a finite-difference approximation used by Dasch.

^{17}Substituting Eq. (B2) into Eq. (B1) and making a change of variable $ u = \delta / ( \Delta r / 2 )$ gives, for the integral portion,

*m*= 0, 1, 2, which is simply two times the definition given by Dasch.

^{17}Solution for

*m*= 0, 1 is given by Eqs. (17) and (18), which corrects the sign error in $ I i j ( 1 )$ [Ref. 26, Appendix 2]. For

*m*= 2, Eq. (B3) is solved analytically [Ref. 25, Eqs. (2.264.3) and (2.261)] and given by Eq. (19). Substituting Eqs. (17)–(19) into Eq. (B1) results in

## REFERENCES

*28th International Symposium on Shock Waves*

*Schlieren and Shadowgraph Techniques: Visualizing Phenomena in Transparent Media*

*N*-waves in air: Modeling and measurements using acoustical and optical methods

*N*-waves in air using an optical schlieren method

*Acoustics: An Introduction to Its Physical Principles and Applications*

*The Fourier Transform and Its Applications*

*The Scientist and Engineer's Guide to Digital Signal Processing*

*N*waves produced by sparks

*Table of Integrals, Series, and Products*