## Abstract

The paper presents approaches based on traditional phase shifting, flexible least-squares, and signal processing methods in dual phase shifting interferometry primarily applied to holographic moiré for retrieving multiple phases. The study reveals that these methods cannot be applied straightforward to retrieve phase information and discusses the constraints associated with these methods. Since the signal processing method is the most efficient among these approaches, the paper discusses significant issues involved in the successful implementation of the concept. In this approach the knowledge of the pair of phase steps is of paramount interest. Thus the paper discusses the choice of the pair of phase steps that can be applied to the phase shifting devices (PZTs) in the presence of noise. In this context, we present a theoretical study using Cramér-Rao bound with respect to the selection of the pair of phase step values in the presence of noise.

© 2006 Optical Society of America

## 1. Introduction

Holographic moiré is an efficient technique for the nondestructive evaluation of rough objects. The method functions by integrating two holographic interferometers within one system. The result is produced in the form of a moiré. However, the complex nature of the interferometer and the beat inherent in the process precludes the use of standard phase shifting procedures to holographic moiré [1–3]. It has been shown previously that the design of special phase shifting procedures to holographic moiré allows for simultaneously determining the information regarding the out-of-plane and in-plane displacement components [1–3]. In this configuration, the phase terms corresponding to the sum of phases (carrier) and difference of phases (moiré) carry information regarding the out-of-plane and in-plane displacement components, respectively. This capability adds new dimension to the measurement of displacement components and in the nondestructive testing of rough objects.

Information carried by the carrier and moiré can be decoded by either the traditional phase shifting approach or by least-squares fit techniques or, by the introduction of signal processing concepts in phase shifting interferometry. In the traditional approach, we refer to three-frame, five-frame, and seven-frame algorithms [1–2] which allow for accommodating dual PZTs. By least-squares fit technique, we extend the least squares fit data reduction method proposed by Morgan [4] and Grievenkamp [5] to accommodate dual phase shifting devices PZTs. On the other hand, the signal processing approach allows for applying high resolution frequency estimation techniques such as polynomial rooting (annihilation method) [6], MUltiple-SIgnal Classification (MUSIC) [7], Minimum-Norm (min-norm) [8], Estimation of Signal Parameters via Rotational Invariance Technique (ESPRIT) [9], and the Maximum-Likelihood Estimation (MLE) [10] to optical configurations including two PZTs.

However, these approaches have inherent limitations and cannot be applied straightforward to configurations such as holographic moiré for the simultaneous extraction of two orthogonal displacement components. For instance, in the traditional approach, the retrieval of information is not straightforward and the phase map corresponding to the carrier yields a phase pattern corrupted by moiré and vice-versa. Hence, additional step needs to be associated to the measurement to yield uncorrupted phase maps. On the other hand, in both the least squares fit technique and the signal processing approach, though the phase steps imparted by the PZTs can be arbitrary, the selection of phase steps in the presence of noise is crucial to the successful implementation of the concept.

In this context, the objective of this paper is to discuss in detail the constraints encountered while following these approaches. We believe that understanding the conditions that are necessary to employing these algorithms is of paramount importance for the successful implementation of the proposed concept. Since, the signal processing approach is the most efficient among these approaches, we present a statistical analysis using Cramér-Rao lower bound [11] for the selection of the optimal pair of phase steps in the presence of noise. Once the phase steps have been estimated within an allowable accuracy, the Vandermonde system of equations [12] can be designed for the determination of phases.

As a case study, we will consider the five-frame algorithm (traditional approach) [2] and the annihilation method (spectral estimation approach) [6]. The performance of other estimators such as MUSIC, MIN-NORM and ESPRIT can also be studied in a similar way. The present study is expected to provide the guidelines for applying dual phase shifting procedures to holographic moiré.

The paper is organized as follows. Section 2 presents the three approaches and discusses the constraints associated with each of these approaches. Section 3 presents the Cramér-Rao bound (CRB) analysis for the estimation of phase steps and studies the performance of the annihilation method in estimating the phase steps with respect to noise.

## 2. Dual phase shifting interferometry: methods and their limitations

#### 2.1. Tradition approach: a five-frame algorithm

Figure 1 shows the configuration for the holographic moiré setup.

The intensity equation for moiré is given by [2]

where, *I _{dc}* is the mean intensity,

*V*is the visibility, and φ

_{1}(

*P*) and φ

_{2}(

*P*) are the interference phases at a point

*P*on the object surface along the two arms of the interferometer. The sum of phases can be extracted by introducing appropriate phase shifts in the two arms of the interferometer and recording the corresponding intensities. Equation (1) thus becomes

where, α and β are the phase steps introduced in the two arms of the interferometer. In order to recover the sum of phase term Φ_{+} = φ_{1}(*P*)+ φ_{2}(*P*), we apply pairs of phase steps (-2α, -2β), (-α,-β), (0,0), (α,β), and (2α,2β) to the PZTs. Frames *I*
_{1},*I*
_{2},*I*
_{3}, *I*
_{4}, and *I*
_{5} corresponding to these phase steps are recorded in the computer [1, 2]. Assuming, α = β = π/2, the term corresponding to the sum of phases is given by

where, χ is some integer constant. Similarly, in order to recover the difference of phase term, Φ_{-} = φ_{1}(*P*) - φ_{2}(*P*), we apply pairs of phase steps (-2α,2β), (-α,β), (0,0), (α, -β), and (2α, -2β) to the PZTs. Frames *I*
_{1}, *I*
_{2}, *I*
_{3}, *I*
_{4}, and *I*
_{5} corresponding to these phase steps are recorded in the computer. Assuming, α = β = π/2, the term corresponding to the difference of phases is given by

### 2.1.1. Constraints in tradition approach

Although α = β = π/2 in Eqs. (3) and (4) minimizes the first order calibration errors in PZTs, the conditions in Eqs. (3) and (4) are the main constraints when it comes to extracting the wrapped phases. For instance, in Eq. (3) whenever, φ_{1}(*P*) - φ_{2}(*P*) = (2_{χ} + 1)π, the equation becomes indeterminate. Because of the arctangent operator this indeterminacy manifests itself as a discontinuity and Φ_{+} jumps by either +2π or - 2π depending upon the sign. Hence, while extracting the wrapped sum of phases, the fringes corresponding to moiré are also seen to modulate the wrapped phase pattern, Φ_{+}.

The above mentioned phenomenon is better understood by simulating the moiré fringes. Let the phase terms recorded at a pixel (*nʹ*,*jʹ*) on the CCD be written as

where, (*n*
_{0}, *j*
_{0}) is the origin for the intensity image of *n* × *j* pixels corresponding to phase in Eq. (5a). In Eq. (5), *n*
_{0} is identical in both equations and the centers of the concentric fringes are offset by *p*
_{0} in y-direction only. Here, Φ_{ran1} and Φ_{ran2} represent the random phase terms (0 to 2π) because of the rough nature of the object surface. Assuming the visibility *V* to be unity, Figs. 2(a) and 2(b) show the fringe pattern (512 × 512 pixels) corresponding to Eq. (2), under the assumption Φ_{ran1} = Φ_{ran2} = 0 and Φ_{ran1} = Φ_{ran2} ≠ 0, respectively. The wrapped phase maps corresponding to Φ_{+} using Eq. (3) for Figs. 2(a) and 2(b) are shown in Figs. 3(a) and 3(b), respectively. The figures show that the information carried by the sum of phases is corrupted by the moiré fringes. The figures are plotted without taking into consideration the constraints in Eq. (3). Figure 4 shows a typical plot along a row in Fig. 3(a). From the plot it can be observed that whenever, φ_{ran1} - φ_{ran2} = (2_{χ} + 1)π, the discontinuity (shown by R1) is ±2π in Φ_{+}. However, this discontinuity can be removed by processing the wrapped phase term Φ_{+} using a computer. An efficient way to perform this task is to resample the wrapped phase using the cosine operator. Since, at φ_{1}(*P*) - φ_{2}(*P*) = (2_{χ} + 1)π, there is a jump of ±2π, from basic trigonometry we get cos [φ_{1}(*P*) - φ_{2}(*P*) ± [2_{χ} +1)π] = cos [φ_{1}(*P*) - φ_{2}(*P*)]. Hence, the discontinuity is removed. In Fig. 4
*f*(*x*) represents the continuous function obtained from *f*(*x*) = cos(Φ_{+}). However, to make the procedure compatible with most of the commercially available unwrapping softwares, the computer generated phase steps can be imparted to *f*(*x*). For instance, phase shifts of 0, π/2, π, and 3π/2 can be applied to extract the wrapped phase. Figures 5(a)–5(b) shows that the discontinuities due to constraints in Eq. (3) are removed.

Similarly, when it comes to extracting the information carried by moiré using Eq. (4), we observe that the phase map is corrupted by carrier fringes. The phase maps corresponding to Φ_{-} in Figs. 2(a)– 2(b) are shown in Figs. 6(a)– 6(b), respectively. The plot in Fig. 7 drawn across a central row shows that whenever φ_{1}(*P*)+ φ_{2}(*P*) = (2_{χ} + 1)π, the discontinuity (shown by R2) is ±2π in Φ_{-}. This discontinuity is removed in a similar way as explained in the previous paragraph. Figure 7 shows the continuous function *g*(*x*) obtained from *g*(*x*) = cos(Φ_{-}). Applying phase shifts of 0, π/2, π, and 3π/2 to *g*(*x*) and solving the four phase shifted images results in the estimation of wrapped difference of phases. Figures 8(a)–8(b) show the wrapped difference of phases obtained for fringes in Figs. 2(a)–2(b), respectively.

The same phenomenon is also observed for the three-frame [1] and seven-frame [2] algorithms. The point which needs to be emphasized here is that the straightforward adaptation of these algorithms to dual phase stepping is not possible. Moreover, these algorithms are sensitive to nonsinusoidal wavefronts (a consequence of detector nonlinearity or multiple reflections inside the laser cavity, or the phase shifter itself). These algorithms also do not offer any flexibility either in the selection of phase step values α and β or in the use of non-collimated wavefronts for phase shifting.

#### 2.2. Flexible least-squares method

The use of least squares fitting techniques has assumed greater significance in simple phase shifting interferometry [4, 5] because of its ability to allow arbitrary phase steps. The concept of least-squares fit has primarily been shown to be effective in configurations involving a single PZT. It can, however, be extended to holographic moiré shown in Fig. 1, albeit, with some constraints. The equation for holographic moiré defined by Eq. (2) can be written as

where, γ = *I _{dc}*

*V*cosφ

_{1}(

*P*), δ = -

*I*

_{dc}*V*sinφ

_{1}(

*P*), μ =

*I*

_{dc}*V*cosφ

_{2}(

*P*), and

*v*= -

*I*

_{dc}*V*sinφ

_{2}(

*P*). Here,

*N*is the number of data frames. Since Eq. (6) is linear with respect to unknown coefficients

*I*, γ, δ, μ, and

_{dc}*v*, we can use least-squares technique to minimize

*E*(

*P*), defined as

For the best fit, the error function in Eq. (7) should be minimized (min*E*(*P*)). This is done by setting the first derivative of *E*(*P*) with respect to the unknown coefficients *I _{dc}*, γ, δ, μ, and

*v*, written as

*∂E*/

*∂I*,

_{dc}*∂E*/

*∂*γ,

*∂E*/

*∂δ*,

*∂E*/

*∂μ*, and

*∂E*/

*∂v*, respectively, equal to zero. The resulting equations can be written in the matrix form as

**B**=

**A**

^{-1}

**X**, where,

**B** = [*I _{dc}* γ δ

*μ*

*v*]

^{T}, and

**X**= [∑

*I*∑

_{n}*I*cos α

_{n}_{n}∑

*I*sin α

_{n}_{n}∑

*I*cos β

_{n}_{n}∑

*I*sinβ

_{n}_{n}]

^{T}. Here the simulation is carried out from

*n*= 0 to

*n*=

*N*- 1. The solution to the above matrix equation results in the determination of the unknown coefficients

*I*, γ, δ,

_{dc}*μ*, and

*v*, and subsequently, in the determination of φ

_{1}and φ

_{2}. The sum and difference of phases can then be obtained using

### 2.2.1. Constraints in flexible least-squares method

The main constraint in flexible least-squares method is that the phase steps should be carefully selected such that matrix **A** is non-singular. It has also been observed that the use of sequential phase steps such as *n*α and *n*β in matrix **A** results in the determinant to be zero. In practical situations although the successive phase steps cannot be exact multiples of *n*α and *n*β, the resulting matrix in this case will be nearly singular or poorly conditioned. Therefore, if the least squares technique is applied, then only few arbitrary phase steps should be selected such that the matrix **A** is non-singular. For instance, if thirteen data frames are acquired and, α = 40° and β = 20°, then five data frames corresponding to frames *n* = 0,2,5,8 and 10 can be selected in a non-regular order. This would result in matrix **A** to be non-singular. It is also observed that even if the matrix **A** is not well conditioned, mathematical tools such as MATLAB [13] may still compute the inverse of the matrix **A**. Therefore, it is always advisable to check whether the matrix is well conditioned or not.

Hence, the well know generalized data reduction technique proposed by Morgan [4] and Grievenkamp [5] cannot be extended straightforward to holographic moiré, since utmost care is required in the selection of the pair of phase steps. The other method by which the generalized holographic moiré can be realized is by designing the Vandermonde system of equations. A Vandermonde matrix usually arises in the polynomial least squares fitting, Lagrange interpolating polynomials, or in the statistical distribution of the distribution moments [14–16]. The solution of *N* × *N* Vandermonde matrix requires *N*
^{2} operations. The advantage of Vandermonde matrix is that its determinant is always nonzero (hence invertible) for different values of *n*α and *n*β. Hence, the matrix for determining phases φ_{1} and φ_{2} can be written in the form

where, (α_{0},β_{0}), (α_{1},β_{1}),..,and (α_{N-1},β_{N-1}) are phase steps for frames *I*
_{0}, *I*
_{1}, *I*
_{2},…,and *I*
_{N-1}, respectively. The phase distribution φ_{1} and φ_{2} are subsequently computed from the argument of ℓ and ℘, where ℓ = 0.5*I _{dc}*

*V*exp(

*j*(φ

_{1}) and ℘ = 0.5

*I*

_{dc}*V*exp(φ

_{2}). Here, ∗ denotes the complex conjugate.

A point to note is that in phase shifting interferometry using the least-squares fit method, the first step involves the exact determination of the phase steps while the second step involves the determination of the phase distribution. The method proposed by Lai and Yatagai [16] for the estimation of phase steps in real time cannot be applied in a straight-forward manner here as the procedure is applicable to a single PZT based configuration only. In such a case the algorithm suggested by Patil *et al*. [17] using stochastic approach seems to be the most appropriate. The algorithm however is sensitive to noise and can only work well if the data is assumed to have a high signal to noise ratio (SNR). To overcome this problem, high resolution techniques have been proposed. Limitations associated with the use of high-resolution techniques is the topic of the next subsection.

#### 2.3. Signal processing approach

In the signal processing approach, a parallelism is drawn between the frequencies that are present in the spectrum and the phase steps imparted to the PZTs. The problem therefore reduces to estimating the frequencies or in other words the phase steps. Once the phase steps have been estimated, the Vandermonde system of equation shown in Eq. (10) can be applied for the extraction of phase distributions. Patil *et al*. [6–10] have introduced signal processing algorithms to configurations involving single and multiple PZTs. The advantage of the these approaches lies in their ability to compensate for the errors arising due to non-sinusoidal wavefronts and PZT miscalibrations. These algorithms also allow the use of spherical wavefronts and arbitrary phase steps between 0 and π radians. The number of data frames that these methods require is equal to at least twice the number of frequencies that are present in the spectrum. Hence, for dual PZTs and for harmonic κ = 1 (assuming sinusoidal wavefronts), we need at least *N* = 10 data frames. Since, discussion of these signal processing algorithms is beyond the scope of this paper, let us consider, as a case study, the case of one such algorithm based on the design of an annihilation filter.

### 2.3.1. Annihilation filter method

In the annihilation filter method, we transform the discrete time domain signal *I _{n}* in Eq. (2) into a complex frequency domain by taking its Z-transform. Let the Z-transform of

*I*be denoted by

_{n}**I**(

*z*). The objective is to design another polynomial

**P**(

*z*) termed as annihilation filter which has zeros at frequencies associated with

**I**(

*z*). This in turn would result in

**I**(

*z*)

**P**(

*z*) = 0. The phase steps α and β are estimated by extracting the roots of the polynomial

**P**(

*z*).

In brief, let us consider the presence of multiple-order harmonics in the wavefront and hence, Eq. (2) for κ^{th} order harmonics can be written as

$$\phantom{\rule{.5em}{0ex}}\sum _{k=-1}^{-\kappa}{b}_{k}\mathrm{exp}\left[\mathrm{jk}\left({\phi}_{2}+\mathrm{n\beta}\right)\right]\phantom{\rule{.2em}{0ex}}+\sum _{k=1}^{\kappa}{b}_{k}\mathrm{exp}\left[\mathrm{jk}\left({\phi}_{2}+\mathrm{n\beta}\right)\right]\phantom{\rule{.2em}{0ex}};$$

where, *a _{k}* and

*b*are the complex Fourier coefficients of the κ

_{k}^{th}order harmonic;

*j*= √-1. Let us rewrite Eq. (11) in the following form

where, ℓʹ_{k} = *a _{k}*exp(

*jk*φ

_{1}),

*u*= exp(

_{k}*jk*α), ℘

_{k}=

*b*exp(

_{k}*jk*φ

_{2}), and ϑ

_{k}= exp(

*jk*α). The task now reduces to designing an annihilation filter of the form

$$=\sum _{k=0}^{4\kappa +2}{P}_{k}{z}^{k}$$

The multiplication in the frequency domain corresponds to discrete convolution in the time domain. Thus, the discrete convolution of the polynomial **P**(*z*), represented as **P**
_{n}, with the intensity signal, represented as **I**
_{n}, vanishes identically, and can be written as

Hence, finding the roots of the polynomial **P**(*z*) yields the phase steps α and β. Since the measurement is sensitive to noise this method necessitates acquiring additional data frames as phase steps cannot be estimated reliably at lower SNR’s with minimum required (*N* = 10) data frames. To reduce the effect of noise, a denoising procedure is suggested in Reference [6] so that phase steps can be estimated even at lower SNR’s. The reader is referred to Reference [6] for details on the annihilation technique and the denoising method.

### 2.3.2. Constraints in signal processing approach

Of the three approaches, the signal processing approach seems to be the most promising as it is effective in minimizing some of the commonly occurring systematic and random errors during the measurement. However, the constraint in the signal processing approach is the selection of the pair of phase steps that are applied to dual PZT’s. In the presence of low signal-to-noise ratio (SNR), a careful selection of the pair of phase steps with limited number of data frames is necessary. Since, noise plays an important role in the successful implementation of the concept, a sound knowledge of the allowable phase steps is of prime importance. Hence, there is a need to study the statistical characteristics of phase shifted holographic moiré (which is independent of the estimator itself) using Cramér-Rao bound. In the next section we present the CRB for holographic moiré in the context of the selection of the phase steps. Once the CRB is obtained the study will focus on the performance evaluation of the annihilation filter method.

## 3. Optimizing phase shifts by signal processing approach

This section provides useful guidelines to optimizing the selection of phase step values obtained by signal processing approach. This study can be performed by deriving the CRB for holographic moiré. The CRB provides valuable information on the potential performance of the estimators. The CRBs are independent of the estimation procedure and the precision of the estimators cannot surpass the CRBs. In this context, we will first derive the CRB for the phase steps α and β as a function of SNR and *N*. Although the phase steps α and β can be arbitrary for pure intensity signal, the central question is: what is the smallest difference between the phase steps α and β that can be retrieved reliably by any estimator as a function of SNR and *N*? Hence, we need to determine all the allowable values of phase steps α and β at a particular SNR and *N*. The reason why we are interested in computing the smallest phase step difference is to provide an experimentalist with a fair idea of phase steps which must be selected at a particular SNR and *N*, so that the phase values φ_{1} and φ_{2} can be estimated reliably. Suppose that inadvertently, the two phase steps are chosen very close to one another and the measurements are performed at a low value of SNR. In such as case, the phase values φ_{1} and φ_{2} cannot be estimated reliably. However, for the same value of SNR, if the two phase steps α and β are far apart, the phase values φ_{1} and φ_{2} can be reasonably estimated. Hence, there is a need to set up guidelines for the selection of the phase steps.

In order to respond to this query, we will derive the mean square error (MSE) for difference in phase steps α – β as a function of SNR which will indicate the theoretical Cramér-Rao lower bounds to the MSE. The Cramér-Rao bound can then be compared with the MSE obtained by any estimator (in the present case we will study the annihilation filter method). For this, we perform 500 Monte-Carlo simulations at each SNR and the separation between the phase steps α and β is varied from 0 to 100%.

Two scenarios are possible. First, the MSE obtained from the estimator is below the CRB. We attribute this observation to the non-reliability in the estimated phase steps and discard those values of separation between α and β in which the MSE is below the CRB. Second, the MSE of the estimator is above the CRB. In such a case, we infer that closer the MSE of the estimator is to the CRB, the more efficient is the estimator. Of course, for an unbiased estimator, the MSE obtained by the estimator can never reach the theoretical CRB for limited number of sample points. In the present study, we derive the CRB for the phase steps and compare the performance of the annihilation filter method with the CRB.

#### 3.1. Cramér-Rao bound for holographic moiré

Let **I**ʹ = **I** + η) denote the vector consisting of measured intensities in the presence of noise η, where **I** = [*I*
_{0}
*I*
_{1}
*I*
_{2}∙*I*
_{N-1}]^{T}. The vector **I**ʹ is characterized by the probability density function *p*(**I**ʹ;Ψ) = *p*(**I**ʹ), where Φ is the set of unknown parameters in the moiré fringes. In the present case, Φ = (*I*
_{dc1}
*I*
_{dc2},*V*
_{1},*V*
_{2},φ_{1},φ_{2},φ_{1},φ_{2},α,β)^{T}, where subscripts 1 and 2 refer individually to the two arms of the interferometer. Note that in Eq. (2), it was assumed that *I*
_{dc1} = *I*
_{dc2} = *I _{dc}* and

*V*

_{1}=

*V*

_{2}=

*V*. Therefore,

*I*

_{dc1}and

*I*

_{dc2}and are the average intensities while

*V*

_{1}and

*V*

_{2}are the fringe visibilities in the two arms. Now, if $\widehat{\Phi}$ is an unbiased estimator of the deterministic Ψ, then the covariance matrix of the unbiased estimator,

*E*{$\widehat{\Psi}$ $\widehat{\Psi}$

^{T}}, where

*E*is an expectation operator, is bounded by its lower value given by [11]

where, **J** is the Fisher Information matrix. This matrix is defined as [9]

In other words, if ψ_{r} (any rth element of *r ^{th}*) is an unbiased estimator of deterministic ψ

_{r}, based on

**I**ʹ, whose CRB is given by

where, **J**
^{-1}
^{r,r} is the (*r*, *r*) element in matrix **J**
^{-1}. Assuming η is the additive white Gaussian noise with zero mean and variance σ^{2}, the probability density function *p*(**I**ʹ;Ψ) is defined by the mean and variance of the noise. Thus, the joint probability of the vector **I**ʹ is given by

It is well understood 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}. Simplifying Eq. (18) by taking the logarithmic function (note that log *p* is asymptotic to *p*(**I**ι; Ψ)), we obtain

Differentiating Eq. (19) with the *r ^{th}* element in Ψ, say ψ

_{r}, we obtain,

Similarly, Eq. (19) can be differentiated with respect to ψ_{s}. Therefore Eq. (15) for the *r ^{th}* and

*s*element in Ψ can be written as

^{th}Equation (21) can be simplified into the following compact form

Finally, the lower bounds are given by variance, var($\widehat{\psi}$
_{r}) ≥ **J**
^{-1},_{r,r}. In the present example the typical Fisher Information matrix will be

where, the subscripts in Eq. (23) represent derivatives with respect to the unknown parameters. Hence, to determine the allowable values of phase steps α and β as function of SNR and *N*, let us first look at the following simple calculation. Suppose that, we have a variable *q* = Ψ^{T}
**U**, which represents the difference between the phase steps α and β, then the variance of q is given by

$$=E\left[{\mathbf{U}}^{T}\left(\Psi -\hat{\Psi}\right){\left(\Psi -\hat{\Psi}\right)}^{T}\mathbf{U}\right]$$

$$={\mathbf{U}}^{T}E\left[\left(\Psi -\hat{\Psi}\right){\left(\Psi -\hat{\Psi}\right)}^{T}\right]\mathbf{U}$$

$$\ge {\mathbf{U}}^{T}{\mathbf{J}}^{-1}\mathbf{U}$$

The problem, therefore, narrows down to selecting a matrix **U** such that the minimum difference between α and β is computed. Since all the parameters are unknown, the matrix **U** is given by **U** =[0 0 0 0 0 0 - 1 1]^{T} .

#### 3.2. Results of CRB analysis

We now perform 500 Monte-Carlo simulations at each SNR and study the performance of the annihilation filter method in retrieving the difference of the phase steps at a pixel point. During the analysis the phases φ_{1} and φ_{2} are kept the same as described earlier in Eqs. (5) and (5b), respectively. Additive white Gaussian noise having SNR between 0 and 70 dB is considered during the simulations. Our analysis consists of two parts.

During the first part, we determine the MSE of β – α as a function of the SNR and the difference between the phase steps α and β as percentage of α. During the analysis we choose α = π/6 and vary the difference of α and β between 0 to 100% of α. Given the fixed value of α, this analysis indicates the reliable values of β which can be selected for a particular SNR and *N*. It is important to note that the value of α can be selected anywhere between 0 and π radians. However, to keep the analysis simple we have selected just one value of α in the present study. The allowable values of β can be determined by comparing the plot of MSE of β – α, obtained using the annihilation filter method for 500 iterations at each SNR, with respect to the mean square error obtained from the CRB given by Eq. (23). The values in the plot where the MSE of β – α, obtained from the annihilation filter method, goes below the MSE obtained using the CRB, are considered as prohibited for the selection of separation between the phase step values α and β. The test is carried out for data frames *N* = 11, 15, and 20. The study shows that separation between the phase steps can be reduced as the number of data frames increases. Moreover, the incorporation of denoising procedure substantially improves the allowable lower range of phase step separation. For instance, Fig. 9(f) shows that for the separation between α and β as 70% of a, and for SNR 30 dB, the MSE is 0.01 rad^{2} using the annihilation method.

In the second analysis, we perform 1000 simulations at a pixel to compare the bounds in the retrieval of phase steps α and β obtained using the annihilation filter method, with theoretical bounds obtained by the CRB for various phase step values in the presence of noise. Since the previous analysis showed that the denoising procedure is effective, the present study is performed with the incorporation of a denoising step and for *N* = 20. Figure 10 shows that as the separation between the phase steps α and β increases, their values can be reliably estimated at lower SNRs. Figures 10(e)–10(f) show that as the separation between the phase steps increases the bounds obtained by annihilation filter method reaches the theoretical bounds given by the CRB at much lower SNR as compared to that obtained in Figs. 10(a)–10(d). This analysis once again reemphasizes the fact that larger the separation between the phase steps, more reliably is the phase obtained at lower SNRs. Similar analysis can also be performed for other algorithms.

## 4. Conclusion

To conclude, this paper has presented three approaches namely, the traditional approach, the least-squares method, and the signal processing (annihilation filter method) approach in dual phase shifting interferometry. The traditional approach suggests that an additional processing step needs to be incorporated to extract the wrapped phases. The study of the least squares fit and the signal processing approaches reveals that a proper choice of the pair of phase steps is of paramount importance. The analysis using Cramér-Rao bound can act as a guideline to select the optimal pair of phase steps in the presence of noise. We believe that a thorough understanding of the issues associated with each of these approaches will pave the way for a real-time and automated simultaneous measurement of two components of displacement in holographic moiré.

## Acknowledgments

The financial support of the Swiss National Science Foundation is gratefully acknowledged.

## References and links

**1. **P. K. Rastogi, “Phase shifting applied to four-wave holographic interferometers,” Appl. Opt. **31**, 1680–1681 (1992). [CrossRef]

**2. **P. K. Rastogi, “Phase-shifting holographic moiré: phase-shifter error-insensitive algorithms for the extraction of the difference and sum of phases in holographic moiré,” Appl. Opt **32**, 3669–3675 (1993). [CrossRef]

**3. **P. K. Rastogi, “Phase shifting four wave holographic interferometry,” J. Mod. Opt. **39**, 677–680 (1992). [CrossRef]

**4. **C. J. Morgan, “Least squares estimation in phase-measurement interferometry,” Opt. Lett. **7**, 368–370 (1982). [CrossRef]

**5. **J. E. Grievenkamp, “Generalized data reduction for heterodyne interferometry,” Opt. Eng. **23**, 350–352 (1984).

**6. **A. Patil, R. Langoju, and P. Rastogi, “An integral approach to phase shifting interferome-try using a super-resolution frequency estimation method,” Opt. Express **12**, 4681–4697 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-20-4681. [CrossRef]

**7. **A. Patil and P. Rastogi, “Subspace-based method for phase retrieval in interferometry,” Opt. Express **13**, 1240–1248 (2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-4-1240. [CrossRef]

**8. **A. Patil and P. Rastogi, “Estimation of multiple phases in holographic moiré in presence of harmonics and noise using minimum-norm algorithm,” Opt. Express **13**, 4070–4084 (2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-11-4070. [CrossRef]

**9. **A. Patil and P. Rastogi, “Rotational invariance approach for the evaluation of multiple phases in interferometry in presence of nonsinusoidal waveforms and noise,” J. Opt. Soc. Am. A **9**, 1918–1928 (2005). [CrossRef]

**10. **A. Patil and P. Rastogi, “Maximum-likelihood estimator for dual phase extraction in holographic moiré,” Opt. Lett. **17**, 2227–2229 (2005). [CrossRef]

**11. **D. C. Rife and R. R. Boorstyn, “Single-tone parameter estimation from discrete-time observations,”IEEE Transactions on Information Theory **IT–20**, 591–598 (1974). [CrossRef]

**12. **M. Marcus and H. Minc, “Vandermonde Matrix,”in *A Survey of Matrix Theory and Matrix Inequalities* (Dover, New York), pp. 15–16 (1992).

**13. **http://www.mathworks.com/.

**14. **http://mathworld.wolfram.com/VandermondeMatrix.html.

**15. **K. M. Hoffman and R. Kunze, “Linear Algebra, 2nd ed.,” (Englewood Cliffs, New Jersey: Prentice Hall), 1971.

**16. **G. Lai and T. Yatagai, “Generalized phase-shifting interferometry,” J. Opt. Soc. Am. A **8**, 822–827 (1991). [CrossRef]

**17. **Abhijit Patil, Benny Raphael, and Pramod Rastogi, “Generalized phase-shifting interferometry by use of a direct stochastic algorithm for global search,” Opt. Lett. **12**, 1381–1383 (2004). [CrossRef]