## Abstract

Electronic speckle pattern interferometry is useful for the qualitative depiction of the deformation profile of harmonically vibrating objects. However, extending the process to achieve quantitative results requires unwrapping the phase in the interferogram, which contains significant noise due to the speckle. Two methods to achieve accurate phase information from time-averaged speckle pattern interferograms are presented. The first is based on a direct inverse of the regions within corresponding phase intervals, and the second is based on optimization of four independent parameters. The optimization method requires less time than more commonly used algorithms and shows higher precision of the resulting surface displacement.

© 2016 Optical Society of America

## 1. INTRODUCTION

Time-averaging electronic speckle pattern interferometry (ESPI) is widely used in nondestructive vibration analysis, providing information on vibration amplitude with high spatial resolution. There are many variations of this interferometric technique, such as phase stepping [1–4], pulsed ESPI [5], amplitude fluctuation [6,7], and max–min scanning [8]. In most cases, they require complicated experimental arrangements and precise adjustments in order to enhance the quality of the acquired images and get the necessary parameters for postprocessing of the image. Reference-updating ESPI [9,10] has shown that it is possible to retrieve the out-of-plane vibration amplitude using a simpler setup without any additional requirements such as high stability and the four-frame phase-stepping algorithm.

The goal of this work is to extend the existing methodology of the reference-updating ESPI to real-world applications using novel approaches of phase retrieval, which result in high-resolution full-field vibration data. The main achievement of this study is the development of two image processing methods for time-averaging ESPI. One of them uses a common model of data processing via direct the inverse of the function describing spatial intensity distribution. The second method is based on global multi-objective parameter optimization and, to the best of our knowledge, represents a novel approach of the data processing in ESPI.

## 2. THEORY

In this section, theory from [11] is presented with subsequent derivations of the formulas for the intensity of the subtracted frames. The experimental setup is shown in Fig. 1.

A piezoelectric transducer (PZT) attached to a mirror is introduced in the reference arm in order to produce a varying optical path difference. Light reflected from the object passes through a semi-transparent glass and, combined with the reference beam, illuminates the detector. The object is sinusoidally excited by a shaker. Polarizing beam splitters control the purity of the polarization, since only the beams with collinear polarization vectors will interfere.

The ESPI system comprises two light beams, a beam reflected from the object and a reference beam, produced by the same laser source. The object beam phase ${\varphi}_{\mathrm{o}}(\mathrm{r\u20d7})$ is random due to the roughness of the object, while the reference beam phase ${\varphi}_{\mathrm{r}}(\mathrm{r\u20d7})$ may vary randomly or may be spatially smooth. Ground glass is introduced before the beam splitter; therefore, the reference beam is speckled. Since the imaging lens collects light from both beams with the same field of view, the mean sizes of the captured speckles for both beams are identical and the intensity statistics are the same. When the two beams are superimposed on the image plane and interfere, the resulting intensity before the deformation is

Time-averaging ESPI is generally used to determine the amplitude values of harmonically vibrating objects. If the excitation force is directly linked to the vibration of the object, the phase change due to a simple harmonic motion with frequency ${\omega}_{0}$ is described as

where $t$ is the time variable and $\xi $ is the contribution of the harmonic motion which contains the following information about the vibration amplitude of the object $A$:The instantaneous intensity of the image at a single point on the detector is given by

Integration of Eq. (4) results in the Bessel function only if the exposure time is an integer multiple of the vibration period $\frac{2\pi}{{\omega}_{0}}$, otherwise it has additional components. Integration errors can be neglected for high-frequency vibration measurements [12]. As a result, there is no longer a dependency on frequency, only on the amplitude of vibration.

In classical time-averaging ESPI, the vibration amplitude is obtained by subtraction of an image of the object under stress from an image of a reference frame [13]. If the time between capturing the reference frame and the deformation frame is relatively large, decorrelation effects can significantly reduce the contrast and the quality of the image. In this case, subtraction of the subsequent frames will result in no intensity modulation in the interferogram. However, introducing a linear phase shift into the reference beam can provide information about the vibration amplitude in the form of dark and bright fringes in the subtracted image [10]. These fringes are often referred to as “correlation fringes,” since they are not physically interferometric fringes.

With the addition of a linear phase shift, normally provided by an oscillating mirror, Eq. (2) is modified to

where the multiplication coefficient $\gamma $ can be regarded as the “phase speed,” since it is measured in ${s}^{-1}$. Taking a series of images, the time-averaged intensity of the $n$th frame measured by the detector isIn this paper, we propose to represent Eq. (9) as the sum of a zeroth-order Bessel function and a zeroth-order Struve function, each multiplied by a coefficient:

where $M$ and $N$ are coefficients depending on the parameters given by the linear phase shift and camera characteristics. Curve fitting of Eq. (9) using the least squares method results in a small fitting error of the order of ${10}^{-8}$. Since the Bessel function ${J}_{0}(\xi )$ equals unity when its argument is zero while the Struve function vanishes, $M$ is defined by the value of ${I}_{\text{sub}}$ when $\xi =0$ asThe second coefficient $N$ can be found only numerically. Elimination of the second term in Eq. (10) increases the error from the order of ${10}^{-8}$ to ${10}^{-4}$ and thus reduces the precision of the method, yet facilitates digital processing.

Consequently, $M$ defines the maximum value of the subtraction intensity, $\mathrm{max}({I}_{\text{sub}})$. The phase coefficient $\gamma $ significantly affects the intensity. Pixels for which the maximum intensity tends to zero can be regarded as low modulated pixels and do not carry any information. Therefore, they may be substituted by neighboring pixel values during image processing. In the assumption above, the initial phase ${\varphi}_{n}$ has been omitted, but when it is taken into account as a random component, it results in drastically changed intensity for a given linear phase shift. In that case, Eq. (11) transforms into

In Fig. 2, the gray vertical line passing through different curves for a given value of $\gamma {t}_{\text{cam}}$ shows that there is no ideal operating mode for various values of ${\varphi}_{n}$, which means that speckle noise is unavoidable.

To conclude, the intensity of the differential image can be approximated by

as suggested in [14]. This simplification gives a good approximation of the intensity profile and following normalization, which eliminates coefficients $\beta $ and $M$, leaves only the Bessel function. Due to the quasi-periodicity of the Bessel function, the phase map resulting from the inversion of the function is “wrapped” over a limited interval. This interval is defined by the distance between 0 and the first root of the Bessel function.Phase unwrapping is a necessary procedure to obtain the full-field data. There are various methods of phase unwrapping, including techniques based on a Hilbert transform [15], Zernike polynomials [16], and the Fourier method [17,18]. More advanced methods require iterative procedures [19], which are reliable but time consuming. Another method is to apply the so-called branch cut method [20], which is generally faster and can unwrap the phase of discontinuous objects. Most of these algorithms are applied to the wrapped phase values obtained using the phase stepping, in which the calculated phase is represented by an arctangent function and is wrapped in the interval from 0 to $2\pi $. Therefore, these algorithms cannot be applied to the inverse of the Bessel function, and another method for the phase retrieval must be found. In what follows, two methods of processing of the experimental data are presented and discussed. The procedure can be called a “phase unwrapping,” although the processing of the phase itself is not performed; the correct phase values are obtained directly from the intensity image.

## 3. EXPERIMENT

To demonstrate the effectiveness of the phase unwrapping methods, ESPI images of a Brazilian Cuíca drum, shown in Fig. 3, were produced. A wooden stick is attached to the center of the membrane, fixed perpendicularly to the drumhead. A shaker sinusoidally excited the free end of the stick along the cylindrical body of the instrument at 1971 Hz. The PZT mirror was driven by a sawtooth waveform, and the camera recorded several frames during each oscillation.

Image acquisition and processing for vibration analysis using subtraction time-averaging ESPI is accomplished in several steps:

- 1. Recording frames of the static object with a given linear phase shift in the reference arm provided by a moving PZT mirror. Subtraction of subsequent frames. Temporal averaging of the images for noise reduction.
- 2. Recording frames of the vibrating object with a given linear phase shift in the reference arm provided by a moving PZT mirror. Subtraction of subsequent frames in order to obtain correlation fringes. Temporal averaging of the images for noise reduction.
- 3. Spatial filtering and normalization by the reference image.
- 4. Frequency filtering and smoothing of the resulting single differential frame.
- 5. Phase unwrapping of the image.

## 4. FILTERING

Averaging $n$ fringe patterns decreases the speckle noise $\sqrt{n}$ times [21]. The multiplication coefficient $\beta $ can be found from the recording of the static object, resulting in ${I}_{\text{sub}}(0)$. Division of the filtered frame from the vibrating state by the filtered frame from the fixed state results in elimination of the coefficient $\beta M$ and normalization of the intensity.

Though spatial filtering decreases the image resolution, it is useful for digital processing of continuous objects. Frequency filters are also applied for speckle noise removal in order to enhance an image and are widely used, despite the fact that the processing time drastically increases. Most frequency filters are based on Fourier filtering. In general, the Fourier transform of a small part of an image is multiplied by a filter function pixel-by-pixel, and then transformed back to the spatial domain. This procedure is done for the whole image, and the filter function defines the method and the effect of filtering [22].

The filtering algorithm applied here involves temporal filtering, spatial filtering, background correction for the image and the reference, normalization by the reference, and Fourier filtering. Temporal filtering is performed by taking the mean value for every pixel through all captured and subsequently subtracted frames, resulting in a single image from the vibrating state and a single reference image with reduced speckle contrast. Further smoothing is accomplished by the spatial filtering. Both images are spatially processed by applying an averaging filter and a median filter on a $3\times 3$ pixel window. Background correction is based on the extraction of the average intensity from a dark fringe area or from the background. Division of the subtracted frame by the reference results in a normalized image free from diffraction patterns caused by dust on the optical components. To further enhance the quality of the image, the background correction is repeated and the image is normalized on the maximum value.

Subsequent Fourier filtering further improves contrast and removes speckles from the image. The proposed method of Fourier filtering is a modification of the algorithm presented in [23]: the image is inverted, so that dark fringes become bright and vice versa, but the Gaussian smoothing step proposed in [23] is omitted, making the algorithm faster.

The image is divided into small overlapping blocks. Every block $P$ is inverted with maintained normalization

where $\mathcal{I}$ is a matrix of ones of the same size as $P$. For every ${P}_{\text{inv}}$, the filter is a Fourier transform of the block to the power $a$: ${|\mathcal{F}\{{P}_{\text{inv}}\}|}^{a}$, where $\mathcal{F}$ represents the Fourier transform. The image is multiplied by the filter in the Fourier domain asThe processing steps are shown in Fig. 4, where it is clearly seen how speckle effects are removed from the differential frame with additional contrast improvement. The processing is complicated by the rough surface of the drum and its inhomogeneity, which produces shadows and large intensity variations. The intensity profile is normalized from 0 to 1 and can be processed as a modulus of the Bessel function without examination of the intensity modulation. The central node on the membrane, where the stick was attached, is also eliminated from the intensity pattern which makes the processing easier.

## 5. PHASE UNWRAPPING

A high-resolution unwrapping of the Bessel fringes was pursued in [2]. However, it takes advantage of the phase stepping and can process only conveniently oriented correlation fringes, which makes it rather difficult to use in many cases.

In this section, two methods of phase unwrapping are presented. The first is based on a direct inverse of the regions of the image corresponding to particular phase intervals as a simple substitution of already existing methods. The second utilizes an optimization function that depends upon the continuity of the surface and intensity.

#### A. Peak Direct Inverse

The direct inverse function of $|{J}_{0}(\xi )|$ applied to every pixel of the image does not result in the correct phase values. The intensity can be processed unambiguously only if the corresponding phase values belong to a single interval between a root of the Bessel function and a neighboring extreme point. Therefore, knowledge of the phase interval for each pixel removes phase ambiguities in image processing. To assign phase intervals, bright and dark fringes are processed separately to exclude errors due to the random noise and background effects.

First, the intensity image is thresholded in order to separate bright correlation fringes corresponding to different phase intervals. The threshold value is based on the background noise. For each area, its local maximum is compared to maxima of $|{J}_{0}(\xi ))|$ function. Each area is numerated by $(2m-1)$, where $m$ is the number of the closest Bessel peak value. If the local maximum is lower than the maximum of the given Bessel lobe, the whole region is normalized on this value. This procedure enhances low-intensity fringes which are affected by the background noise.

Black fringes are defined through the inversion of the binary mask obtained via thresholding. Each region is expanded, and the order numbers $2m$ are given with respect to the neighboring bright fringes. The interval assignment is shown in Fig. 5. Bright fringes get odd numbers, the corresponding phase values are located in the intervals between zeros of the Bessel function, and dark fringes are numbered with even numbers so that the corresponding phase values are located in the intervals between the extremes of the Bessel function.

The separation between dark and bright fringes is performed because the algorithm picks only the local extreme of the image, and this procedure is relatively simple. This results in a map with assigned orders, as shown in Fig. 6. In this image, the topology of the output phase is visible and clear.

Intensity values of the local maxima and minima can be directly converted into phase values. The conversion results in a cloud of points which can be interpolated assuming continuity of the object without any ambiguity of the side lobe order, as shown in Fig. 7(a). The first-order region is processed pixel-by-pixel, since it does not contain ambiguities. The phase is obtained using the grid of recovered values by interpolation and is shown in Fig. 7(b). This method is robust and works for various shapes of fringe patterns, providing information about the absolute value of the vibration amplitude.

The three-dimensional profile of the surface vibration amplitude is shown in Fig. 8(a). The absolute displacement is defined by Eq. (3). The corresponding intensity pattern is shown in Fig. 8(b).

This method can be improved by a final pixel-by-pixel inversion applied while searching near the value of the obtained phase. However, this procedure significantly slows the process.

#### B. Optimization Method

The optimization method takes advantage of the continuity of the object surface as well as the Bessel order assignment described above. The optimization is accomplished on a grid of pixel values with arguments chosen based on the minimum fringe spacing. Small grid size results in increased evaluation time and required computer memory. The chosen pixels form a matrix ${I}_{\mathrm{s}}$ of the size $[m,n]$ of the experimental values. $\mathrm{\Phi}$ is the matrix of the phase values which are calculated and adjusted through the minimization of the objective function ${F}_{\text{obj}}$. Multiple constraint optimization is accomplished using the Rosenbrock optimization algorithm [24], which performs global optimization over the set of parameters and general constraints. Since there are several criteria to be assessed, the optimization is multi-objective and the objectives must be taken with different weight, according to their importance and contribution. Thus, the objective function is presented as the sum of four errors by

where ${C}_{1}$, ${C}_{2}$, ${C}_{3}$, and ${C}_{4}$ are weighting coefficients.The first error ${E}_{1}$ defines the convergence of the parameters to the experimental data by

The second error ${E}_{2}$ ensures that each parameter stays within a given phase interval and appears only if the value borders are violated and is given by

The third and the fourth errors track the continuity of the second derivative of the phase profile in horizontal and vertical directions as

In the current algorithm, the initial values of the weighting coefficients are $C=1$, ${C}_{2}=1000$, ${C}_{3},{C}_{4}=1$. When the interval condition is met and ${E}_{2}$ turns to zero, the contribution of smoothness and convergence become more important, ${C}_{1}$ becomes 5, and the coefficients ${C}_{3}$ and ${C}_{4}$ equal 90.

Figure 9(a) shows the grid of unwrapped points, followed by interpolation, presented in Fig. 9(b). In comparison with the direct inverse method, the phase points are homogeneously distributed over the pixel range, which provides good spatial resolution.

The profile of the surface vibration amplitude is shown in Fig. 10(a). The reconstructed intensity is shown in Fig. 10(b).

This method can be improved by adjusting the grid size to the fringe spacing in order to place an equal number of points between the dark fringes and thus to enhance spatial resolution.

#### C. Results

The optimization results in a smoother profile which is provided by the objective function and in an interferogram of a better quality with reference to the direct inversion. The initial intensity profile is noisy, which results in asperities in the resulting phase map derived through the direct inversion. Figure 11 shows the profiles, taken through the upper part of the intensity image for the initial normalized and spatially filtered image and for the two calculated intensity images. Zero displacement is depicted by high peaks, and the side lobes are enhanced in both cases. Zero transitions can get lost during optimization if the smoothness of the resulting surface prevails over the convergence factor. However, the peak direct inverse method allows simple and accurate tracking of the vibration nodes of the surface.

## 6. CONCLUSIONS

The work described here introduces two methods for determining the absolute values of the vibration amplitude from an electronic speckle pattern interferogram. The optimization procedure provides results spanning the entire area of an image, whereas the peak direct inverse method concentrates on the extreme of the interference pattern and nodal parts of the vibrating surface. This indicates that the direct inverse method should perform better on higher amplitude variations where fringes are close to each other. The optimization algorithm works at least twice as fast as the direct inverse algorithm, which makes it very attractive in real-time applications.

One of the main challenges faced by the phase unwrapping algorithm is the correct detection of the Bessel function orders based on the maximum values within the region. The issue of the sign ambiguity of the amplitude cannot be solved by the current methods and experimental arrangement, however, high-quality interferometric images can be obtained for further processing.

## Funding

Seventh Framework Programme; (FP7-PEOPLE-2013-ITN; 605867).

## Acknowledgment

This research work has been funded by the European Commission within the ITN Marie Curie Action project BATWOMAN under the Seventh Framework Programme (EC grant agreement no. 605867). The authors thank A. Mayer for technical support.

## REFERENCES

**1. **K. Creath, “Phase-shifting speckle
interferometry,” Appl. Opt. **24**, 3053–3058
(1985).

**2. **D. N. Borza, “Full-field vibration amplitude recovery
from high-resolution time-averaged speckle interferograms and digital
holograms by regional inverting of the Bessel
function,” Opt. Lasers Eng. **44**, 747–770
(2006). [CrossRef]

**3. **C. Tang, J. Zhang, C. Sun, Y. Su, and K. L. Su, “Electronic speckle pattern
interferometry for fracture expansion in nuclear graphite based on PDE image
processing methods,” Proc. SPIE **9525**, 952521 (2015). [CrossRef]

**4. **C. Song, M. V. Matham, and K. H. K. Chan, “Speckle referencing: digital speckle
pattern interferometry (SR-DSPI) for imaging of non-diffusive
surfaces,” Proc. SPIE **9524**, 952421 (2015). [CrossRef]

**5. **J. N. Petzing, “Vibration analysis using fast speckle
metrology,” Proc. SPIE **3745**, 134–140
(1999).

**6. **W. Wang and C. Hwang, *The Development and Applications of Amplitude
Fluctuation Electronic Speckle Pattern Interferometry Method*,
Recent Advances in Mechanics
(Springer,
2011).

**7. **X. Dai, X. Shao, Z. Geng, F. Yang, Y. Jiang, and X. He, “Vibration measurement based on
electronic speckle pattern interferometry and radial basis
function,” Opt. Commun. **355**, 33–43
(2015). [CrossRef]

**8. **E. Vikhagen, “Nondestructive testing by use of TV
holography and deformation phase gradient
calculation,” Appl. Opt. **29**, 137–144
(1990). [CrossRef]

**9. **B. Pouet, T. Chatters, and S. Krishnaswamy, “Synchronized reference updating
technique for electronic speckle interferometry,”
J. Nondestr. Eval. **12**, 133–138
(1993). [CrossRef]

**10. **T. R. Moore and S. A. Zietlow, “Interferometric studies of a piano
soundboard,” Acoust. Soc. Am. **119**, 1783–1793
(2006). [CrossRef]

**11. **T. R. Moore and J. J. Skubal, “Time-averaged electronic speckle pattern
interferometry in the presence of ambient motion. Part I. Theory and
experiments,” Appl. Opt. **47**, 4640–4648
(2008). [CrossRef]

**12. **W. C. Wang, C. H. Hwang, and S. Y. Lin, “Vibration measurement by the
time-averaged electronic speckle pattern interferometry
methods,” Appl. Opt. **35**, 4502–4509
(1996). [CrossRef]

**13. **R. Jones and C. Wykes, *Holographic and Speckle Interferometry*,
2nd ed. (Cambridge
University, 2001).

**14. **T. Statsenko, V. Chatziioannou, and W. Kausel, “Interferometric studies of the Brazilian
Cuíca,” in *Proceedings of the Third
Vienna Talk on Music Acoustics* (2015),
pp. 153–157.

**15. **L. Krzemien and M. Lukomski, “Algorithm for automated analysis of
surface vibrations using time-averaged DSPI,”
Appl. Opt. **51**, 5154–5160
(2012). [CrossRef]

**16. **L. Ernesto, M. Espinosa, H. M. Carpio Valadez, and F. J. Cuevas, “Demodulation of interferograms of closed
fringes by Zernike polynomials using a technique of soft
computing,” Eng. Lett. **15**, 99–104
(2007).

**17. **C. Trillo, A. F. Doval, F. Mendoza-Santoyo, C. Pérez-López, M. de la Torre-Ibarra, and J. L. Deán, “Multimode vibration analysis with
high-speed TV holography and a spatiotemporal 3D Fourier transform
method,” Opt. Express **17**, 18014–18025
(2009). [CrossRef]

**18. **J. Muñoz Maciel, F. J. Casillas-Rodríguez, M. Mora-González, F. G. Peña Lecona, V. M. Duran-Ramírez, and G. Gómez-Rosas, “Phase recovery from a single
interferogram with closed fringes by phase
unwrapping,” Appl. Opt. **50**, 22–27
(2011). [CrossRef]

**19. **J. M. Bioucas-Dias and G. Valadão, “Phase unwrapping via graph
cuts,” IEEE Trans. Image Process. **16**, 698–709
(2007). [CrossRef]

**20. **M. J. Huang and J. K. Liou, “Retrieving ESPI map of discontinuous
objects via a novel phase unwrapping algorithm,”
Strain **44**, 239–247
(2008). [CrossRef]

**21. **J. W. Goodman, *Laser Speckle and Related Phenomena*
(Springer,
1975).

**22. **C. Tang, L. Wang, H. Yan, and C. Li, “Comparison on performance of some
representative and recent filtering methods in electronic speckle pattern
interferometry,” Opt. Lasers Eng. **50**, 1036–1051
(2012). [CrossRef]

**23. **C. Li, C. Tang, H. Yan, L. Wang, and H. Zhang, “Localized Fourier transform filter for
noise removal in electronic speckle pattern interferometry wrapped phase
patterns,” Appl. Opt. **50**, 4903–4911
(2011). [CrossRef]

**24. **H. H. Rosenbrock, “An automatic method for finding the
greatest or least value of a function,”
Comput. J. **3**, 175–184
(1960). [CrossRef]