1. Introduction

In interferometry, the incorporation of a phase shifting device PZT has become a common practice for determining the phase of the interference beams. Typically, while using a single PZT, one has access to only one information from the phase, for instance, the out-of-plane displacement. However, if multitude of information is desired, incorporation of multiple PZTs seems to be one of the possible solutions. This concept was first put forward by Rastogi [1–2] who showed that two orthogonal displacement components can be measured by placing dual phase shifting devices in an appropriate holographic moiré configuration.

However, the estimation of multiple phase information is still prone to errors that may arise during measurements. For instance, the phase computation is sensitive to the nonlinear response of the PZT to the applied voltage. The error in the measurement of phase can compound further due to the nonsinusoidal intensity profile and noise. Furthermore, if multiple phases are to be estimated then incorporation of multiple PZTs is necessary, and in such a case there is a need to propose a method which addresses these concerns. Moreover, one desired feature while using a phase shifting device is to compute the phase step values in real time, rather than relying on the previously calibrated phase shift values. This step ensures that any error that may arise due to ageing of PZT or hysteresis is minimized.

In this context, the paper shows an experimental investigation using a generalized log-likelihood function for retrieving multiple phase information in holographic moiré in the presence of nonlinear response of the PZTs to the applied voltages. The proposed method is also effective in compensating the nonsinusoidal waveforms that can occur due to detector nonlinearity, or multiple reflections inside the laser cavity, or the phase shifter itself. The other salient feature of the proposed algorithm lies in its ability to allow the use of arbitrary phase steps and non-collimated waveforms.

2. Dual phase shifting in holographic moiré

Holographic moire is essentially an incoherent superposition of two interferometric systems which result in the formation of a moiré pattern. Incorporation of two phase shifting devices in a holographic moiré configuration enables determining two orthogonal displacement components simultaneously. For this, two distinct phase steps are applied simultaneously to the two arms of the moiré setup. The reader is referred to Reference 1 for details of the optical setup. The method consists of acquiring data while voltages are simultaneously applied to the two PZTs. The fringe intensity recorded at any point (x,y) of the n^th data frame is given by

where, I _dc1 and I _dc2 represent local background intensities and, γ ₁ and γ ₂ represent fringe visibilities for each interferometric system. In Eq. (1), pairs φ ₁, φ ₂ represent the phase differences, and ${\alpha }^{\prime }$_n, ${\beta }^{\prime }$_n represent phase shifts in the two systems, respectively. Since, the acquired intensity on the CCD is also sensitive to the random errors that can occur during the measurement, the term η_n represents additive white Gaussian noise with mean zero and variance σ ². For N data frames, the true phase shifts {${\alpha }^{\prime }$_n, ${\beta }^{\prime }$_n} applied to the PZTs due to the nonlinear characteristics of the PZTs can be written as

where, ${\alpha }^{\prime }$_n and ${\beta }^{\prime }$_n are the polynomial representation of the phase steps α and β, respectively, and ε ₁₂ and ε ₂₂ are the nonlinear error coefficients of the two PZTs response to the applied voltages. In what follows, we present a nonlinear maximum likelihood method for the estimation of phase steps and nonlinear coefficients [3]. For this let us rewrite Eq. (1) in a compact form

where, I = [I ₀ I ₁ ...I _N-1]^T, η = [η ₀ η ₁ ...η _N-1]^T, T is the transpose operator, and matrix S(ξ) is expressed as S(ξ) = [1 s _α s*_α s _β s*_β] where ξ = {α, ε ₁₂,β, ε ₂₂} is the parameter set, and signal vectors s _α and s _β are given by s _α(n) = exp(j${\alpha }^{\prime }$_n), s _β(n) = exp(jβ_n ) where $j=\sqrt{-1}$. The matrix C in Eq. (3) is the coefficient matrix given by C = [c ₀ c ₁ c*₁ c ₂ c*₂]^T with c ₀ = I _dc1 + I _dc2 and elements c ₁ and c ₂ are given by c ₁ = I _dc1 γ ₁ exp(jφ ₁) and c ₂ = I _dc2 γ ₂ exp(jφ ₂). Since the noise is assumed to be additive white Gaussian with variance σ ², the probability density function of the data vector I¯, parameterized by ξ, is given by p(I¯; ξ)

where, H is the Hermitian transpose. Since the likelihood function of the data L(I¯, ξ) is proportional to p(I¯; ξ), the maximum likelihood estimate of ξ, is obtained by maximizing L(I¯, ξ). From Eq. (4) we observe that, this maximization is equivalent to minimizing the expression D(ξ)=[I¯-S(ξ)C]^H [I¯-S(ξ)C]. We also note that, this expression is in the form of a standard least squares if S(ξ) is assumed to be a known matrix. Therefore, D(ξ) is minimized if C is

Substituting Eq. (5) into the expression for D(ξ), we obtain

The vector ξ, which appears in the matrix S(ξ) is nonlinearly related to I¯. The estimation of nonlinear parameters in ξ, is obtained by minimizing the expression for D(ξ) in Eq. (6), which in turn leads to maximizing the expression

It can be observed from Eq. (7) that obtaining the optimum ξ, will require a multi-dimensional search over the parameter space. One can use any global search method for the multi grid search. Once the parameter set ξ, is estimated, we can determine the values of ${\alpha }^{\prime }$_n and ${\beta }^{\prime }$_n, followed by estimation of phases φ ₁ and φ ₂ by any least square estimation method.

3. Evaluation of the method

The method is tested experimentally in a holographic moiré setup [1]. We first calibrate the two PZTs responses to the applied voltages and identify the linear regions in both the PZTs. In this linear region we sample nonuniform points to simulate nonlinear response. The phase steps in the two arms of the holographic moiré are selected as α = 25°(0.4363rad), β = 50°(0.8727rad); and nonlinear coefficients ε ₁₂, and ε ₂₂ are selected as 0.15 and 0.20, respectively. Thus, the robustness of our method can be tested for known magnitude of phase steps and nonlinearity. We acquire twenty five data frames on a 576×768 pixel CCD camera, while the phase steps are applied simultaneously to both the PZTs. Our objective is to first find the parameter set ξ = {α, ε ₁₂,β,ε ₂₂}by maximizing the expression in Eq. (7), so that the phase distributions φ ₁ and φ ₂ can be obtained. For this any global optimization technique can be applied. In the present evaluation, we apply the Probabilistic Global Search Lausanne (PGSL)[4] to maximize the expression in Eq. (7). The advantage of PGSL lies in its ability to estimate the global optimum without the requirement to make initial guesses for the parameter set. However, one needs to fix the lower and upper bounds of the elements in the parameter set ξ. In this analysis the lower and upper bounds were set as ξ_lb = {0,0,0,0} and ξ_ub = {π/2,0.3, π/2,0.3}, respectively.

Figure 1(a) shows a typical moire fringe pattern obtained experimentally. Figure 1(b) shows the corresponding holographic fringe pattern in one arm of the interferometer. It is well known that the retrieved phase steps at each pixel will not be constant because of noise and other errors, and hence studying the distribution of the retrieved phase steps can be one way to compute the phase step and nonlinear coefficient. We thus first randomly select a pixel on a 576 × 768 CCD camera, and for N = 25 data frames, we sample the intensity values and compute the phase steps at this pixel location. This process is repeated for 1600 randomly selected pixels on the data frame and a histogram is plotted to study their distribution. The histogram for phase steps α and β is drawn by dividing the phase steps between 0° and 180° into 721 equal intervals/bins. Similarly, the histogram for nonlinear coefficients ε ₁₂ and ε ₂₂ is drawn between 0 and 0.3 with 100 equal bins. The number of occurrence of each parameter is put in the corresponding bin. Figures 2(a)–2(b) show the histograms for α and β obtained using our method. The phase steps corresponding to bins consisting of maximum appearances in Figs. 2(a) and 2(b) are α = 24.7°(0.4311 rad) and β = 48.9deg(0.8535 rad), respectively. The corresponding values of nonlinear coefficient ε ₁₂ and ε ₂₂ are 0.145 and 0.187, respectively, as shown in Figs. 2(c) and 2(d). The plots show that the parameter set ξ, is obtained efficiently. Once these parameters are obtained, we apply the least squares method to obtain C as given in Eq. (5) from which the phases φ ₁ and φ ₂ are obtained from the argument of c ₁ and c ₂, respectively, since recall that c ₁ = I _dc1 γ ₁ exp(jφ ₁) and c ₂ = I _dc2 γ ₂ exp(jφ ₂) in Eq. (3). The wrapped phases φ ₁ and φ ₂ are shown in Figs. 3(a) and 3(b), respectively.

The efficiency our method is further tested by comparing the variance on each parameter of the proposed estimator with the Cramér-Rao lower Bound (CRB). The following section discusses in detail the statistical efficiency of the method in the estimation of phases φ ₁ and φ ₂.

 

Fig. 1. Fringe patterns obtained experimentally from optical set up in Reference 1: a) Moiré, and b) fringes corresponding to the phase difference φ ₁.

Download Full Size | PDF

4. CRB for interferometry model with multiple PZTs

Let us derive the CRB for Eq. (1). The probability density function is given by p(I¯;ψ) = p(I¯) where, I¯ = [I¯₀ I¯₁ ...I¯_N-1]^T, ψ is the set of unknown parameters given by ψ = {I _dc1,γ ₁,α,ε ₁₂,φ ₁,I _dc2,β,ε ₂₂,φ ₂}. It is well known that the CRB for each unknown parameter can be determined by observing the diagonal elements of the inverse of the Fisher Information matrix, J ^-1, which is given by [5]

 

Fig. 2. Histogram representations for phase steps a) α = 25°, and b) β = 50°; and for nonlinear coefficients c) ε ₁₂ = 0.15, and ε ₂₂ = 0.20.

Download Full Size | PDF

 

Fig. 3. Wrapped phases a) φ ₁, and b) φ ₂.

Download Full Size | PDF

where E is the expectation operator. The (r, s)^th element in the Fisher information matrix is given by

The CRB of each parameter ψ_r in the parameter set ψ is given by the r^th diagonal element of the inverse of the Fisher information matrix, i.e., ${\mathbf{J}}_{\mathit{\text{rr}}}^{-1}$. From the theory of the CRB, if $\widehat{\psi }$ is an unbiased estimator of the deterministic ψ, then the covariance matrix of the unbiased estimator, E{$\widehat{\psi }$ $\widehat{\psi }$ ^T} is bounded by CRB and is given by E{$\widehat{\psi }$ $\widehat{\psi }$ ^T} ≥ J ^-1. Hence, the mean square error (MSE) of each estimated parameter is lower bounded by

Closer the value of the E{$\widehat{\psi }$_r ²} to ${\mathbf{J}}_{\mathit{\text{rr}}}^{-1}$, the more efficient is the estimator. Finally, to observe the efficiency of the proposed method we compare the MSE of the phase φ ₁ with the CRB given for φ ₁ versus signal-to-noise ratio (SNR). We have selected same phase steps as considered in the experiment, and phase differences φ ₁ and φ ₂ are selected as 20°(0.3490 rad) and 45°(0.7854rad), respectively. Figure 4 shows the comparison of MSE in the estimation of phase φ ₁ obtained by our method with the CRB of φ ₁ for various levels of nonlinearity in Eq. (2). We have chosen the value of nonlinear calibration errors ε ₁₂ = ε ₂₂ = 0.01,0.1,0.2 and 0.4 respectively. From Figs. 4(a)-4(d) we observe that the MSE of phase φ ₁ obtained using our method is close to the CRB. Similar results are obtained for phase φ ₂ and hence the corresponding plots are not shown. Although the maximum-likelihood method is known to be optimally efficient and reaches CRB, with a limited number of data samples, the MSE cannot reach CRB but can only be close to the CRB. Figures 4(c) and 4(d) respectively, show a certain bias in the MSE for phase φ ₁ for SNR values from 46dB and 40dB, onwards. This error can be attributed to the inefficiency of the global search method in estimating the correct phase step response for highly nonlinear cases. Similarly, we can also demonstrate the efficiency of the proposed method in the presence of higher order harmonics.

 

Fig. 4. a) Comparison of the MSE in the estimation of phase φ ₁ of the our proposed method with CRB for various levels of nonlinearity. Same nonlinearity is assumed in both the PZTs a) ε ₁₂ = ε ₂₂ = 0.01, b) ε ₁₂ = ε ₂₂ = 0.1, c) ε ₁₂ = ε ₂₂ = 0.2, and d) ε ₁₂ = ε ₂₂ = 0.4.

Download Full Size | PDF

5. Conclusion

To conclude, we have proposed a novel approach for estimating multiple phases in an interferogram in the presence of PZT nonlinearities, harmonics, and noise. The method allows the flexibility of applying arbitrary phase steps. Since the non-linear phase steps are estimated pixel wise in real time, the errors occurring due to the ageing of the PZTs and the hysteresis associated with them can be minimized. In addition we have demonstrated the efficiency of the proposed method by comparing its MSE with the CRB. Further research will focus on designing a more efficient search technique so that the highly nonlinear case of PZT response could also be dealt with.