1. Introduction

Phase Shifting has been widely adopted in optical interferometry for retrieving phase information encoded in the interference fringes [1–5]. The technique functions primarily by acquiring a number of intensity images with phase increments between successive frames. These phase increments are generally applied using a piezoelectric device (PZT). Although three frames are normally sufficient to compute the phase distribution, the measurement process is sensitive to various systematic [6–18] and random errors [19–24]. The large number of frames has been shown to decrease the sensitivity to these errors.

This method, commonly known as temporal technique for phase measurement, can be broadly classified to fall into two categories: conventional and generalized. While the step widths are typically multiples of π/2, in the former case, these are arbitrary in the latter case [25–26]. As far as conventional techniques are concerned, several algorithms have been reported to compensate for errors that are introduced during phase evaluation. The systematic errors occurring in the phase shift device have been investigated. The corresponding algorithms insensitive to these errors are described in the literature. Most of these algorithms aim at minimizing the errors arising due to the nonlinear characteristic inherent in the piezoelectric device PZT. These algorithms termed as self-compensating [1,8–10,12,22,27] have proved to be effective while compensating deterministic phase shift error, even until third order nonlinearities. Some of the algorithms minimize the effect of other systematic errors such as optical system aberrations [5,18], parasitic fringes [6,28], photodetection errors [29,30] and quantization errors [17,31]. Random errors due to thermal or shot noise have also been observed to influence the phase measurement [21 ]. Random errors in measurement of phase can also crop up because of turbulent and laminar air flows in the optical path [30], vibrations [19–20,22], electronic noise [32] originating in the photodetection during amplification, and optical noise commonly known as speckle [33]. The statistical properties of phase-shift algorithms have been investigated for the case of additive Gaussian intensity noise. It has been shown that some advanced systematic-error-compensating algorithms are sensitive to random noise [21].

The error in phase measurement can also arise if the fringe profile is not purely sinusoidal because of detector nonlinearity or multiple-beam interference. The algorithm proposed by Hibino et al [14] offers the possibility to reduce the calibration error in the presence of higher order harmonics. Hibino et al [14] and Surrel [11] showed that 2κ + 3 and 2κ + 2 samples are respectively necessary to minimize the effect of the κ ^th order harmonic, with the phase step 2π/(κ+2) between samples. The restriction on the free use of phase steps represents a drawback and a possible limitation for the methods wide spread application. The other important limitations of the algorithm reside in its inability to function in configurations requiring multiple PZTs.

While conventional phase shifting has demonstrated its ability to overcome various sources of systematic errors, generalized phase shifting technique though flexible as far as the selection of phase steps is concerned hasn’t been successful in overcoming most of the systematic errors discussed earlier. Although the algorithm proposed by Carré [1] is insensitive to first order calibration error and enables the use of arbitrary phase steps, the study imposes a restriction on the choice of the phase step value in order to minimize certain other errors referred in [34]. The algorithm proposed by Lai and Yatagai [35] requires an additional optical setup for the generation of Fizeau fringes and can only handle fringes that are not straight and equally spaced [36]. Farrell and Player [37–38] describe an ellipse fitting technique that uses Bookstein algorithm [39] - a forced fit method - to fit the data set into an ellipse. The data set fits into a hyperbole for noisy data [40]. The statistical approach proposed by Cai et al [41] to extract reference phase steps imposes requirements on the spatial resolution of the CCD camera (which in turn dictates the sampling number of the speckle intensity). The max-min algorithm proposed by Chen et al is computationally exhaustive and requires a large number of data frames (15 or more) for reliable operation [42]. The algorithm based on iterative least squares estimation proposed by Han and Kim is suitable only for compensating linear first order errors in PZT and requires the initial guess of the phase step to be as close as possible to its true value [43]. The concept of identifying the reference phase steps proposed by Jüptner et al is suitable for only linear miscalibration error and can produce large number of outliers for noisy interferograms [44]. Although the five frame algorithm proposed by Stoilov and Dragostinov gives the flexibility of using arbitrary phase steps, its use is limited only to linear errors encountered in the PZT [45]. The stochastic approach [46] to detect the PZT nonlinearity in PZT proposed by Patil et al [47] is sensitive to random sources of error.

To the best of authors knowledge no generalized phase shifting technique has been reported so far which at a time caters to the presence of multiple PZTs in an optical configuration, is insensitive to errors due to nonsinusoidal waveforms, to PZT miscalibration and to additive white Gaussian intensity noise, and has the ability to work with diverging as well as converging beam. The use of multiple PZTs has been introduced in an optical configuration by Rastogi [48–49] for automating the extraction of information conveyed by both carrier and moiré in a holographic moiré setup. With the field of phase shifting interferometry advancing rapidly, the possibility of incorporating multiple PZTs in an interferometric configuration is a significant conceptual progress which should be useful for developing multiple channel measurement capabilities in moiré based interferometers.

The aim of this paper is to propose an integral generalized approach which first caters to the presence of multiple PZTs in an optical configuration, second is capable of handling harmonics, third is insensitive to errors arising from PZT miscalibration, fourth is applicable to spherical beams, and fifth is a robust performer even in the presence of white Gaussian intensity noise. We intend to apply the super-resolution frequency estimation approach which basically takes the advantage of the fact that a polynomial can be associated to a periodic intensity fringe pattern recorded by a CCD. Another polynomial termed as annihilating filter [50] is designed which has zeroes at the frequencies that are associated with the polynomial for the intensity fringes. The discrete convolution of the filter and the intensity fringe yields zero. The spectral information embedded in the signal corresponds to the phase steps imparted by the multiple PZTs in the interferometric setup. Hence, the parametric estimation of this annihilating filter yields the desired spectral information embedded in the signal, which in our case are the phase steps. Although the discrete Fourier transform is an efficient tool for the estimation of well separated frequencies, the separation of closely spaced frequencies in the presence of noise and less number of samples can be handled efficiently with high resolution techniques [51]. This method for frequency estimation is also sometimes referred to as super-resolution technique because of its ability to resolve spectral lines (frequencies) separated in frequency f = ω/2π by less than 1/N cycles per sampling interval (here, ω and N refers to angular frequency and number of samples, respectively), which is the resolution limit for the classical periodogram-based methods, such as Fourier transform. As a case study, we show the effectiveness of our algorithm to extract the dual phase steps in holographic moiré. Since, the proposed technique functions by extracting the phase step at each pixel location of the acquired frames, this method allows the use of diverging as well as converging beams. The advantage of using any arbitrary phase step between 0 and π radians overcomes the limitation exhibited by previously reported methods. We also show that the robustness of our algorithm to additive white Gaussian noise can be enhanced by incorporating a denoising step using the concept of truncated singular value decomposition [52–53].

In what follows, we first present a brief introduction to holographic moiré. Section 3 presents the design of annihilating filter for phase step determination. Section 4 deals with the identification of the harmonic content of the waveform. Section 5 shows the denoising procedure. Section 6 shows the simulation results for phase step extraction and retrieval of phases in holographic moiré in the presence of noise.

2. Holographic moiré

Although holographic interferometry has proved its effectiveness for measuring displacements along the line of sight, it does not provide means for obtaining whole field distribution of displacements orthogonal to the line of sight. This limitation has led to the development of more complex interferometers to obtain direct visualization of displacement components and their derivatives. The method referred to as holographic moiré has proved its usefulness in determining in-plane displacements [54], slopes and curvature [55], and difference displacements [56] of objects.

The application of phase-stepping to holographic moiré is a relatively complex procedure in which a dual-phase stepping is applied in the two arms of a holographic moiré setup [48–49]. Besides providing automation to the method, one of the significant advances of holographic moiré is that it offers retrieving simultaneously the phase information carried on moiré and carrier at any point on the object surface. The recorded fringe intensity at a point (x,y) for m^th phase step is given by

where, a_k and b_k is the complex Fourier coefficient of the k^th order harmonic, $i=\sqrt{-1}$ and I_dc is the local average value for intensity; and pairs φ ₁ and φ ₂, and, α and β, represent phase differences and phase shifts, respectively, in the two arms of the holographic moiré setup. The coefficients a_k and b_k are in fact real and are related to the appropriate choice of phase origin at a point where the intensity reaches a maximum.

Phase shifted holographic moiré offers the advantage of extracting simultaneously the phase values φ ₁ and φ ₂. Phase information corresponding to moiré, Φ, and carrier, Ψ, are consequently obtained by and which, for example, in a holographic moiré configuration are directly related to in-plane and out-of-plane displacement components, respectively.

3. Theory of estimation of phase steps using a super-resolution method

In order to explain the basis of the spectral estimation procedure, let us write Eq. (1) as

where, ℓ_k = a_k exp(ikφ ₁), u_k = exp(ikα), ℘_k = b_k exp(ikφ ₂), v_k = exp(ikβ); superscript * denotes the complex conjugate; φ ₁ and φ ₂, and, α and β, are the unknown parameters; I_n(x,y;m) corresponds to the n^th frame, where n = 0 refers to the data frame corresponding to the first phase shifted intensity pattern. Equation (4) represents the complex valued sinusoidal signals where α,….,κα; -α,….,-κα; β,….,κβ; -β,….,-κβ represent the frequencies embedded in the signal, and by estimating them, the phase steps α and β, imparted by the PZTs can be determined. Frequency estimation of sinusoids is a classical problem in spectral estimation and as mentioned in Section 1 can be better handled using a super-resolution method. In this method we first transform the discrete time domain signal I_n in Eq. (4) into a complex frequency domain by taking its Z-transform. Let the Z-transform of I_n be denoted by I(z). The objective here is to design another polynomial P(z) termed as annihilating filter which has zeros at frequencies associated with I(z), which in turn would result in I(z)P(z) = 0. Since the frequencies translate into zeros, spectral estimates feature high resolution. In what follows we will derive the expression for P(z) and explore the condition for which the multiplication of P(z) and I(z) is zero.

The Z-transform of Eq. (4) is written as

which after expansion and simplification can be written as

where,

Taking into account the terms D’ s and P’ s and substituting them in Eq. (6), we obtain

Equation (7) can be represented as

where, D(z) represents the numerator and P(z) the denominator in the right hand side of Eq. (7), with the latter being given by

Given the complex form of D(z), let us consider for the sake of simplicity writing only its first term,

From Eq. (10) we observe that some coefficients of C₁(z), specifically, C1₀, C1₁, C1₂……C1_4κ, and C1_N, C1_N+1,C1_N+2……C1_N+4κ depend upon the unknown amplitude and frequency. We observe that the coefficients C1_4κ+1, C1_4κ+2, C1_4κ+3……, C1_N-1 are identically zero for N ≥ 4κ + 1. This condition is at the essence of the design of annihilating filter P(z). In what follows we show that the (4κ+1)^th degree polynomial P(z), when discretely convolved with I(z) yields zero.

Since the multiplication of two signals in frequency domain corresponds to its discrete convolution in time domain, inverse Z-transform of Eq. (8) gives

where, ⊗ represents the convolution operator. In Eq. (11), p_n is the inverse Z-transform of P(z) defined by

where, δ(n) is the unit impulse signal. Equation (11) can be further written as

Using the fact,

we can express Eq. (13) in a matrix form as

We partition the first matrix on the left hand side of Eq. (15), select the middle N-4κ-2 rows equations corresponding to zero row values in the matrix on the right hand side, and form a new matrix as follows:

We can infer from Eq. (16) that the polynomial P(z) is an annihilating filter which when convolved with the moiré intensity signal yields zero. It can further be deduced that at least N ≥ 8κ+2 samples are required to extract the roots of the polynomial P(z). This enables us to find the unknown values u_k and v_k. The phase steps α and β, can hence be computed using α = ℜ(lnu ₁/i) and β = ℜ(lnv ₁/i).

4. Detection of nonsinusoidal waveform and the corresponding harmonic content

The occurrence of nonsinusoidal waveform - a consequence of detector nonlinearity or multiple reflections in laser cavity - results in an error in the computation of phase. This section explains how to deduce the number of harmonics κ present in the signal, which can subsequently be applied in the design of the Vandermonde system of equations in Section 6 for the determination of phase values φ ₁ and φ ₂. Let us rewrite Eq. (1), for the case κ = 1 (i.e. pure sinusoidal wave) and noiseless signal, as

Initially we assume T1≠0 and T2 = 0, a case in which the intensity I(x,y;m) in Eq. (17) corresponds to two wave interferometry. Fourier transform of I(x,y;m) for m = 1,2,……,N should result in three peaks in frequency domain corresponding to I_dc, α and -α. Similarly for T1≠0 and T2≠0, a case in which the intensity I(x,y;m) in Eq. (17) corresponds to holographic moiré interferometry, Fourier transform of I(x,y;m) for m = 1,2,……,N should result in five peaks in frequency domain corresponding to I_dc, α, -α, β, and -β. However, with limited data samples, the resolution of closely spaced frequencies is troublesome in the presence of noise. Also because of the “leakage” phenomenon, energy in the main lobe can leak into the side lobes obscuring and distorting other spectral responses.

Interestingly, the singular value decomposition of R in Eq. (16), represented as R = USV ^T, yields more reliable information than the Fourier transform method regarding the number of harmonics present in the signal, where S is a diagonal matrix with M nonzero and N - M zero singular values; U and V are unitary matrices [57]. When T1 ≠ 0 and T2 = 0, in Eq. (17), and for a noiseless case, the number of nonzero diagonal entries in S is M = 3 (corresponding to I_dc, α, and -α). On the other hand, for T1 ≠ 0 and T2 ≠ 0, the number of nonzero diagonal entries is M = 5 (corresponding to I_dc, α, -α, β, and -β). Hence, if H is the number of PZTs used in the optical setup, the number of harmonics κ in the signal can be determined using M = 2Hκ+1. In the presence of noise, M principal values of S would still tend to be larger than the N-M values which were originally zero. In addition, the M eigenvectors corresponding to the M eigen values of R ^T R are less susceptible to noise perturbations in comparison to the remaining N-M eigenvectors.

Figure 1 illustrates typical singular values for S when the number of samples (data frames) used are N = 25. In this example we assume the presence of second order harmonics in the signal, κ = 2, and the presence of two PZTs in the optical configuration, H = 2. In this case and for a noiseless signal we should expect M = 2Hκ + 1 = 9 diagonal entries of S in Fig. 1 to be significantly larger in magnitude; the N - 9 diagonal entries being zero. Even in presence of noisy data for the case with (SNR = 10 dB), M = 9 principal values are still larger in magnitude than the N - 9 diagonal values. Hence by selecting only those significant values of S (which in this example, corresponds to M = 9), the number of harmonics κ is determined to be κ = 2. The plot thus allows for a reliable estimation of the number of harmonics present in the signal.

 

Fig. 1. Plot showing the magnitude of diagonal values in matrix S versus N data points. From the plot the number of harmonics κ in the signal can be computed using M = 2 Hκ+1. In figure (o) and (+) represent the diagonal values for noiseless and noisy (SNR=10 dB) signals, respectively.

Download Full Size | PDF

5. Denoising the signal

There are various sources of noise in optical interferometry and the relative contribution of each noise source depends on the particular system and its application [58]. Denoising is thus an important step towards reducing the effect of noise. We consider white Gaussian noise to contaminate the signal; the reason being that many distributions can be approximated by a Gaussian function. It should also be noted that many noise sources are known to follow a Gaussian noise source in accordance to the central limit theorem [21]. To enhance the signal in the presence of white Gaussian noise, we apply the concept of truncated singular value decomposition. In the present case, number of data frames ≥ 8κ+2 are needed for applying the denoising procedure.

The concept of TSVD involves the following steps [52–53]. The matrix R is first written in Hankel matrix form, say R̂, and its singular value decomposition R̂ = ÛS̄V̂ ^T shows the nonsingular principal values of S̄ which are significantly different from zero. After setting the non significant N-M singular values of S̄ to zero, a matrix Ŝ is formed. A denoised matrix Z _M =ÛŜV̂ ^T which approximates R̂ in the least squares sense, is then obtained by using the first M columns of Û, Ŝ and V̂; Û and V̂ being unitary matrices. Finally, a denoised signal Ī_n is retrieved by arithmetic averaging along the anti-diagonals (or diagonals) of Z _M using

where, r = max[1, n - number of rows (R̂)+1) and q = min[number of rows (R̂), n]. The denoised signal Ī_n is subsequently applied in Eq. (16).

6. Evaluation of phase steps and phase distributions in presence of noise

The next step after denoising the signal is to extract the phase steps. We present the application of super-resolution technique for the extraction of phase steps by simulating moiré fringes given by

where, (x ₀,y ₀) is the origin of the X×Y pixels for fringe pattern with pitch λ ₁; (p ₀,y ₀) is the origin of X×Y pixels for fringe pattern with pitch λ ₂; ξ is some arbitrary constant; κ = 2 and the phase steps are selected as α = π/4 and β = π/3. A white Gaussian noise with SNR from 0 to 100 dB is added to test the robustness of the proposed concept. Fringes shown in Fig. 2(a)–(c) correspond to Eq. (19) for κ = 1 and for three different noise levels. On the other hand, fringes shown in Figs. 3(a)–(c) corresponds to Eq. (19) for second order harmonic, κ = 2, but for same noise levels as in Figs. 2(a)–(c), respectively. These fringes have been generated under the assumption a _{0, ± 1, ± 2} = b _{0, ± 1, ± 2} = 0.5.

 

Fig. 2. Holographic moiré for κ = 1 and for noise levels (a) 0 dB (b) 10 dB , and (c) 60 dB additive white Gaussian noise.

Download Full Size | PDF

 

Fig. 3. Holographic moiré for κ=2 and for noise levels (a) 0 dB (b) 10 dB , and (c) 60 dB additive white Gaussian noise.

Download Full Size | PDF

Figure 4(a) shows the plot for phase step extraction when the noise level is considered from 0 to 100 dB; eighteen data frames are used to obtain this plot. As can be seen from the plot, the phase steps cannot be accurately determined for SNR level below 40 dB. In order to extract the phase steps reliably at a lower SNR, additional data frames are mandatory so that denoising procedure explained in Section 5 can be successfully applied. Figures 4(b) and 4(c) show a plot for the phase step extraction when twenty seven and thirty six data frames are used in the computation of phase steps and denoising procedure is applied. The plot shows an increase in reliability in the computation of phase steps with a larger number of frames. Large data frames can induce unwanted noise in the measurement, if controlled environment is not used. Our denoising technique works well when random noise follows a Gaussian profile, which is the case most of time. For noise which does not follow a Gaussian profile (air turbulence, or vibrations), reference 30 provides an insight into the measures which can be taken to minimize their influence on the measurement of phase distribution.

 

Fig. 4. Plot of phase steps α and β versus noise when computation is done over (a) eighteen (b) twenty seven, and (c) thirty six frames. Plots in. (b) and (c) are obtained after applying the denoising procedure. Large number of frames are due to the presence of two PZTs (H = 2) and two harmonics κ = 2, which in turn impose lower limit on data samples (4Hκ+ 2) as eighteen for phase step estimation.

Download Full Size | PDF

Once the phase steps are estimated the parameters ℓ_k and ℘_k can be solved using the linear Vandermonde system of equations obtained from Eq. (4). The Vandermonde system of equations always has full rank as long as α_i’s and β_i’s are distinct. The matrix thus obtained can be written as

where (α ₁,β ₁), (α ₂,β ₂), (α ₃,β ₃), .. and (α_N,β_N) are the phase steps for frames I ₀, I ₁, I ₂, .., and I _N-1, respectively. The phases φ ₁ and φ ₂ are subsequently computed from the argument of ℓ₁ and ℘ ₁. Figures 5(a) and 5(b) show typical errors in the computation of phase φ ₁ without and with the application of the denoising procedure, respectively.

 

Fig. 5. Plot of error in the computation of phase φ ₁ when SNR = 30 dB obtained a) without and b) with the application of the denoising procedure.

Download Full Size | PDF

The noise level is assumed to be 30 dB. In this simulation phase steps obtained in Fig. 4(c) are used for the computation of phase. Similarly Figs. 6(a) and 6(b) show errors in the computation of phase φ ₂ without and with the application of the denoising step, respectively, all the other parameters remaining the same. Figure 7 shows the wrapped phases φ ₁ and φ ₂ obtained with the denoising step.

 

Fig. 6. Plot of error in the computation of phase φ ₂ when SNR = 30 dB obtained a) without and b) with the application of the denoising procedure.

Download Full Size | PDF

 

Fig. 7. Wrapped phase distributions φ ₁ (solid line) and φ ₂ (broken line) as functions of pixel position when SNR = 30 dB. Denoising procedure has been applied for obtaining these results.

Download Full Size | PDF

7. Conclusion

To conclude, we have presented a novel integral approach - based on the use of a super-resolution frequency estimation method - to phase shifting interferometry, starting from phase step estimation to phase evaluation at each point on the object surface. A significant advantage of this approach is that it also offers the possibility to identify the harmonic content of the signal. The described method holds promise for determining simultaneously multiple phase distributions in the presence of higher order harmonics and in optical configurations with multiple PZTs. The proposed technique works well with the diverging as well as converging beams since the phase steps are retrieved point by point before being applied to the Vandermonde system of equation. The robustness of the algorithm to white Gaussian noise is enhanced by the incorporation of denoising step before the extraction of phase step. Further research will focus on estimating the lower bounds of the distance between phase steps that can be resolved using the high resolution method.