1. Introduction

In phase shifting interferometry (PSI), the phase is computed by acquiring few to several phase shifted intensity images on the CCD camera and solving the intensity equations to yield the phase information. The phase of interfering beams is varied by a piezoactuator device known as PZT. The adaptability of contemporary phase shifting algorithms to holographic interferometry has been rather straightforward. These algorithms have demonstrated their ability to minimize various systematic [1–18] and random [19–24] sources of errors occurring during measurement.

Holographic moiré [25–29] offers the possibility to measure simultaneously the out-of-plane and in-plane displacement components. The information regarding the out-of-plane and in-plane displacement components is carried by carrier and moiré, respectively. However, this necessitates design of special optical setups which can accommodate dual PZTs. In such cases, the direct adaptability of phase shifting algorithms to compute the phases in real time is a difficult task. Recently, multiple phase extraction using super-resolution technique has been suggested; however, this technique requires incorporation of additional denoising procedures and thus increases the computational cost [30].

The objective of this paper is to propose a novel approach to identify phase steps imparted to each PZT in an optical configuration in the presence of nonsinusoidal waveforms. Nonsinusoidal waveforms can occur because of multiple reflections inside the laser cavity, detector nonlinearity, or error in phase shifter [11]. The phase steps imparted to the PZTs are also not repeatable because of the hysteresis effect or ageing and hence careful calibration needs to be carried out. In such case, it is advisable that phase steps be computed once the intensity images are acquired so as to reduce reliance on previous calibration. Sources of random error which is known to follow Gaussian distribution [21] also make the estimation of phase steps difficult.

To show the feasibility of the proposed concept, in this paper, an algorithm known as minimum-norm [31–32] (commonly known as min-norm) is applied to holographic moiré which handles nonsinusoidal wavefronts and estimates phase steps pixel-wise in the presence of white Gaussian noise [21]. The algorithm provides the flexibility of using arbitrary phase steps, the feature most commonly attributed to generalized phase shifting algorithms [33–34]. The concept basically functions by designing an autocovariance matrix from the acquired moiré intensity recorded on the CCD for N data frames [31–32, 35]. The eigendecomposition of the autocovariance matrix yields the signal- and noise-subspaces. The phase steps are estimated from the noise subspace. In the following section, we show an important property in which noise subspace is orthogonal to sinusoidals present in the spectrum. Since we draw a parallelism between the frequency present in the spectrum and the phase steps employed by the PZT’s, this property is exploited to extract the latter. Once the phase step values of the two PZT’s are estimated pixelwise, the Vandermonde system of equations are applied to estimate the phase distribution.

The next section describes the min-norm algorithm for estimation of phase steps followed by its evaluation in Section 3. In Section 4 the method for phase extraction is presented.

2. Min-Norm algorithm for holographic moiré

The min-norm algorithm has been used extensively in frequency estimation of sinusoidals buried in white noise. Although, both the min-norm algorithm and the multiple signal classification technique (MUSIC) derive frequency estimates from noise subspace, the former offers substantial computational advantage [36–37] with almost the same statistical accuracy. Hence, in the present study min-norm algorithm has been preferred for phase extraction. The algorithm functions by designing an autocovariance matrix from N phase shifted moiré images measured at a point (x, y) on the CCD. Usually, an autocovariance matrix is formulated from small overlapping fragments of the data sets using sliding window technique known as spatial smoothing. Equation governing the spatial smoothing technique will be presented in Section 3. Subsequently, an autocovariance matrix is eigen-decomposed to yield the signal- and noise subspaces. The derivation [36–37] in the following section shows that the signal and noise spaces are orthogonal to each other. Therefore, the phase step values of the two PZT’s are estimated pixel-wise from the noise subspace. The following paragraph explains in detail the procedure.

Two distinct phase steps are applied simultaneously to the two illumination arms of the moiré setup shown in Fig. 1. The method consists of acquiring images while voltages are simultaneously applied to the two PZTs. The recorded fringe intensity at a point (x, y) of the t^th frame is given by [25–29]

where, a_k and b_k are complex Fourier coefficients, κ is the order of harmonics, i=√-1, I_dc is the local average value of 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.

 

Fig. 1. Schematic of the optical setup in holographic moiré.

Download Full Size | PDF

The intensity for holographic moiré in Eq. (1) is assumed to be buried in white Gaussian noise since studies have shown that noise influencing the phase measurement follows a Gaussian profile [21]. Equation (1) in the presence of additive white Gaussian noise can thus be written 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 η the additive white Gaussian noise with mean zero and variance σ ². In Eq. (2), 0 (corresponding to I_dc ); ±α,±2α…..,±κα; ±β,±2β…..,±κβ ; are the frequencies ${\left\{{\omega }_n\right\}}_{n=0}^{4\kappa }$, and by evaluating them the phase steps α and β can be identified. Here, n represents the number of frequency components present in the signal. The first step, common to all the methods [31–32, 35–37] seeking to derive frequency estimates using the subspace based methods, is to form an autocovariance matrix. In the present case, an autocovariance matrix is designed from N recorded phase shifted sequences. The autocovariance of signal I(t) in Eq. (2) is defined as [38]

where, ere, E [·] represents the expectation operator which averages over the ensemble of realizations. For simplicity, let us consider κ=1 and rewrite Eq. (2) as

Similarly, let us write I *(t-p) for κ=1 as

Substituting Eqs. (4) and (5) in Eq. (3), we obtain the following

where, $c_1=I_{\mathit{dc}}{a}_1{e}^{-i{\phi }_1}{e}^{-i\alpha t}+I_{\mathit{dc}}{a}_1{e}^{i{\phi }_1}{e}^{i\alpha t}+I_{\mathit{dc}}{b}_1{e}^{-i{\phi }_2}{e}^{-i\beta t}+I_{\mathit{dc}}{b}_1{e}^{i{\phi }_2}{e}^{i\beta t}$,

$c_2=I_{\mathit{dc}}{a}_1{e}^{-i{\phi }_1}{e}^{-i\alpha t}+a_1^2{e}^{-2i{\phi }_1}{e}^{-2i\alpha t}+a_1{b}_1{e}^{-i({\phi }_1-{\phi }_2)}{e}^{-it(\alpha -\beta )}+a_1{b}_1{e}^{-i({\phi }_1+{\phi }_2)}{e}^{-it(\alpha +\beta )}$,

$c_3=I_{\mathit{dc}}{a}_1{e}^{i{\phi }_1}{e}^{i\alpha t}+a_1^2{e}^{2i{\phi }_1}{e}^{2i\alpha t}+a_1{b}_1{e}^{i({\phi }_1+{\phi }_2)}{e}^{it(\alpha +\beta )}+a_1{b}_1{e}^{i({\phi }_1-{\phi }_2)}{e}^{it(\alpha -\beta )}$,

$c_4=I_{\mathit{dc}}{a}_1{e}^{-i{\phi }_2}{e}^{-i\beta t}+b_1^2{e}^{-2i{\phi }_1}{e}^{-2i\alpha t}+a_1{b}_1{e}^{i({\phi }_1-{\phi }_2)}{e}^{it(\alpha -\beta )}+a_1{b}_1{e}^{-i({\phi }_1+{\phi }_2)}{e}^{-it(\alpha +\beta )}$, and

$c_5=I_{\mathit{dc}}{b}_1{e}^{i{\phi }_2}{e}^{i\beta t}+a_1^2{e}^{2i{\phi }_2}{e}^{2i\beta t}+a_1{b}_1{e}^{i({\phi }_1+{\phi }_2)}{e}^{it(\alpha +\beta )}+a_1{b}_1{e}^{-i({\phi }_1-{\phi }_2)}{e}^{-it(\alpha -\beta )}$.

Let, E[$I_{\mathit{\text{dc}}}^2$ +c ₁]=$A_0^2$, E[$a_1^2$+c ₂]=$A_1^2$, E[$a_1^2$+c ₃]=$A_2^2$, E[$b_1^2$+c ₄]=$A_3^2$, and E[$b_1^2$+c ₅]=$A_4^2$. Therefore,

In Eq. (7), σ ² δ_{p ,0} is the expectation for Gaussian noise η (t) and is given by

In Eq. (8), δ_g,h =1 if g=h ; and δ_g,h =0 otherwise. In practice, expectation E in Eq. (3) is computed by averaging over finite number of frames. If a large number of frames is taken for averaging, the exponential terms containing t in c ₁, c ₂, c ₃, c ₄, and c ₅ will oscillate uniformly between 0 and 2π. In this limit, the expectation of c ₁, c ₂, c ₃, c ₄, and c ₅ will approach zero because $\underset{0}{\overset{2\pi }{\int }}$ e ^iψ d ψ=0. However, if finite number of frames are taken for averaging, expectation c ₁, c ₂, c ₃, c ₄, and c ₅ will have a small finite value different from zero. Hence, for κ harmonics in the intensity, the final derivation of covariance of I(t) is given by

An autocovariance of function I(t) is assumed to depend only on the lag between the two averaged samples. The autocovariance matrix can thus be written as [35–38]

where, I(t)=[I(t-1),…..,I(t-m)], m is the autocovariance length, and r*(-p)=r(p). An autocovariance matrix R _I can be shown to have the form

where, R _s and R _ε are the signal and noise contributions, A _m ×(4 _κ +1)=[a(ω ₀) . . a(ω _4κ)] where for instance element a(ω ₀)consists of m×1 matrix with unity entries corresponding to I_dc ; a(ω ₁)=[1 e^iα . . eiα(m-1)]^T; (·) ^c and (·)^T is the conjugate transpose and transpose operator, respectively, I is the m×m identity matrix; and P _m×m matrix is

In Eq. (11), since R _I is a Hermitian matrix, we have, R _I =R* _I ∈C ^m×m, where C ^m×m is a complex matrix. The matrix R _I is said to be positive semidefinite since its eigenvalues are nonnegative. The singular value decomposition of R _I yields a P matrix whose diagonal elements (also known as eigenvalues) can be ordered as $A_0^2$≥$A_1^2$≥….≥A _4κ² ≥….$A_m^2$ . This decomposition yields significant results as far as identifying the eigenvectors associated with signal and noise subspaces are concerned. The following paragraph explains the signal and noise subspaces.

Let S _m×n=[s₁,s₂,…..s _n ] be the orthonormal eigenvectors associated with $A_0^2$≥$A_1^2$≥…..$A_{4\kappa }^2$. Let G _m×(m-n)=[g ₁,g ₂,…..g _m-n] be the set of orthonormal eigenvectors associated with eigen values $A_{4\kappa +1}^2$ ≥$A_{4\kappa +2}^2$≥…..$A_m^2$. Since APA ^c ∈C ^m×m has rank n (n<m), it has n eigen values and the remaining m-n eigen values are zero. If we further suppose that (A ²,r) is an eigenpair of L∈C ^m×m and W=L+ρ I with ρ∈C, then (A ²+ρ,r) is an eigenpair of W. In consequence we obtain $A_t^2$ =${\mathit{Ã}}_t^2$ +σ ², where t spans from 1,2,3,…,m. We observe that $A_0^2$≥$A_1^2$≥….≥$A_{4\kappa }^2$≥σ ² and $A_{4\kappa +1}^2$ = $A_{4\kappa +2}^2$ =…..=$A_m^2$ =σ ². Following this corollary and from Eq. (11) we get

Last equality in Eq. (13) means that APA ^c G=0 and since AP has full column rank we have

This means that sinusoidals {a(ω_k )${\}}_k^n$ ₌₀ are orthogonal to noise subspace. This also indicates that {g k${\}}_{k=1}^{m-n}$ of G belong to the null space of A ^c (can be written as g k∈N(A ^c ). The dimension of N(A ^c ) is equal to m-n which is also the dimension of range space of G, written as R(G). This fact combined to the observation that A ^c G=0 leads us to

Since, by definition

and R(G)=N(S ^c ); we have, N(S ^c )=N(A ^c ). We further deduce that since R(S) and R(A) are orthogonal complements to N(S ^c ) and N(A ^c ), respectively, it follows that

Observing Eqs. (16) and (17); the subspaces R(S) and R(G) are called signal subspace and noise subspace, respectively.

In practice, since only the estimate R̂_I of R_I is available, only the estimate Ĝ of G can be determined. MUSIC uses m-n linearly independent vectors in R(Ĝ) to obtain the frequency estimates. Using Eq. (14), in root MUSIC, the frequencies are estimated as angular positions of n roots of equation [39]

In Eq. (18), a(z) is obtained from a(ω), and by replacing e^iω by z we get a(z)=[1 z ^-1 . z ^-(m-1)]^T. Since the minimum possible value of m is n+1, it can be observed from Eq. (11) that we need N data frames that have at least twice the number of sinusoidal components in the signal. For example, for κ=2, we have n=9 (since n=2H_κ+1, where H is the number of PZT’s); we need at least eighteen data frames (the dc component is also counted as dc frequency). Hence, the minimum number of data frames required for phase extraction is 2n. If only one PZT is employed for phase shifting in an optical setup, then number of data frames required is 4κ+2.

From Eq. (18), one can also see that since any vector in R(Ĝ) is orthogonal to {a(ω_k )}n _k=0, therefore, without sacrificing the accuracy too much, only one such vector can be used. This would result in substantial computational savings. The statistical accuracy of the min-norm algorithm is similar to that obtained using MUSIC. Hence, the performance of MUSIC algorithm is achieved at reduced computational cost [32]. Min-norm algorithm computes the frequencies as the angular positions of the n roots of the polynomial

where, $\left[\begin{array}{c}1\\ \hat{\mathbf{g}}\end{array}\right]$ is the vector in R(Ĝ) with first element equal to unity, that has minimum Euclidean norm. The following explains the rational behind this specific selection [30–31].

Let the Euclidean norm of a vector be denoted by ‖·‖. Since, in reality only the estimate R̂ _I of R _I is available, only the estimate Ŝ of S can be determined. Let the matrix Ŝ be partitioned as

As, $\left[\begin{array}{c}1\\ \hat{\mathbf{g}}\end{array}\right]\in R\left(\hat{\mathbf{G}}\right)$, Eq. (16) can be written as

Equation (21) can be rewritten using Eq. (20) as

The minimum-norm solution to Eq. (22) is given by

under the assumption that the inverse exists. Given that the identity matrix I can be written as

and that one eigenvalue of I-χχ ^c is equal to 1-‖χ‖2 and the remaining n-1 eigenvalues of I-χχ ^c are equal to unity, the inverse in Eq. (23) exists if and only if

If Eq. (25) is not satisfied, there will be no vector of the form of [1 ĝ]^T in R(Ĝ). Condition in Eq. (25) is equivalent to rank(S̄ ^c S̄)=n which, in turn, holds if and only if

From Eq. (17) it can be observed that any block of matrix S made from more than n consecutive rows should have rank equal to n. Therefore, Eq. (26) is valid as long as N is sufficiently large. This completes the derivation of min-norm frequency estimator defined in Eq. (19). In addition, it has been shown that by using the min-norm vector R(Ĝ) as defined in Eq. (19), the spurious frequency estimates can be reduced, a problem sometimes associated with the MUSIC method.

To apply min-norm for estimating the phase steps, the number of harmonics present in the signal must first be computed so that appropriate value of n is determined. The details on selection of m and N will be explained in next Section. Determining the number of harmonics is a typical problem in signal processing and more details can be obtained from any standard signal processing text book. One way [40], by which the harmonics can be determined is by observing the Singular Value Decomposition (SVD) of R _I matrix in Eq. (10).

3. Simulation of the min-norm algorithm to holographic moiré

The concept is tested by simulating holographic moiré fringe pattern in Eq. (1). Let us rewrite Eq. (1), for the case κ=2 and noise η as

The phase steps are selected as α=π/4, β=7π/18, and phases φ ₁,φ ₂ at (x,y) on the intensity map are defined by

and,

where (p′,p′) is the center for circular fringe pattern corresponding to phase φ ₁, (p″, p′) is the center for circular fringe pattern corresponding to phase φ ₂ and, φ_{R 1} and φ_{R 2} are the random phases. Figure 2 shows the fringe pattern corresponding to Eq. (1) for two different values of κ for different noise levels. Figures 2(a) and (b) show fringe maps for κ=1 ; and, no noise and 10dB SNR, respectively. Figures 2(c) and (d) show fringe maps for κ=2 ; and, no noise and 10dB SNR, respectively.

In practice, an autocovariance matrix R _I is formulated from small overlapping fragments of data sets using sliding window technique known as spatial smoothing. Equation governing the spatial smoothing technique for the design of an autocovariance matrix is given by [35]

where, R̂ _I is the estimate of R _I . The approach which obtains frequency estimates from R̂ _I in Eq. (30) is also called the forward approach. We study the retrieval of phase steps α and β at a pixel (x, y) using this approach in the presence of additive white Gaussian noise with SNR between 0 and 70 dB. Let us assume that the number of frequencies n present in the signal is determined to be nine using the method suggested in Section 2. Hence, the number of harmonics present in the moiré fringes is κ=2, which sets lower limit on the number of data frames as eighteen. Once n is identified, an appropriate value of m must be selected such that m>n . The minimum value of m is n+1 and it is observed that though m>n increases the accuracy of frequency estimates, this is at a higher computational cost. On the other hand, m too close to N does not yield an autocovariance matrix R̂ _I similar to R _I . This in turn results in spurious frequency estimates. Performing eigendecomposition of R̂ _I gives estimates for eigenvectors Ĝ. Phase steps α and β are then estimated using Eq. (19).

 

Fig. 2. Fringe map corresponding to a) κ=1 and pure signal, a) κ=1 and 10 SNR, a) κ=2 and pure signal, a) κ=2 and 10 SNR.

Download Full Size | PDF

Figure 3(a) shows the plot for a case in which the number of data frames N=18 and m=12 is selected. The phase step values α and β at any arbitrary pixel location on the data frame cannot be estimated from this plot. In the second case, we choose the number of data frames to be N=22. Figures 3(b)–(d) show typical plots for phase steps α and β for m=12, 15, and 19, respectively. From these plots it can be observed that even if the number of frames is increased the phase steps can be estimated only when appropriate value for m is selected. Figure 3(b) shows that phase steps α and β can be estimated reliably from 35 dB onwards. In the third case we choose N=26, and plots for m=18 and 23 are shown in Figs. 3(e) and 3(f), respectively. Figure 3(e) shows that phase steps can be estimated reliably for values of SNR 35 dB and above. This shows that with small increase in data frames, a substantial improvement in estimating the phase steps is not obtained. On the other hand, Fig. 3(f) shows that value of m too close to N does not yield phase step estimates.

A better estimate for the values of frequencies can be obtained if the sample autocovariance matrix is modified using the following expression [41]

The approach which obtains frequency estimates from R̂ _I in Eq. (31) is called the forward-backward approach. It has been shown that the accuracy of frequency estimates can be enhanced if this sample autocovariance matrix is applied.

Figure 4(a) shows the plot for data frames N=18 and m=12. Phase step values α and β can be estimated from this plot for values of SNR of 40dB and above. In the second case, we choose the number of data frames N to be 22. Figures 4(b)–(d) show typical plots for phase steps α and β with m=12, 15, and 19, respectively. Figures 4(b)–(c) shows that the values of α and β can be estimated at much lower SNR (25 dB and above). Figure 4(d) shows that value of m is critical and its value too close to N does not yield result. In the third case, we choose N=26 ; and corresponding plots for m=18 and 23 are shown in Figs. 4(e) and 4(f), respectively. From these plots, we observe that phase steps α and β can be reliably estimated from Fig. 4(e) as compared to Fig. 4(f). We also observe from Fig. 4(e) that phase steps can be estimated at much lower SNR (20dB and above) as compared to that obtained using other plots. Again, Fig. 4(f) shows that the values of m too close to N do not yield reliable phase step estimates. From these three cases we can conclude that phase step values α and β can be reliably estimated at lower SNR’s with increase in data frames N and m not too close to N and n.

 

Fig. 3. Plot for phase step values α and β (in degrees) obtained at an arbitrary pixel location on a data frame for different values of N and m using the forward approach. During the simulation the phase steps are assumed to be α=45° and β=70°.

Download Full Size | PDF

 

Fig. 4. Plot for phase step values α and β (in degrees) obtained at an arbitrary pixel location on a data frame for different values of N and m using the forward-backward approach. During the simulation the phase steps are assumed to be α=45° and β=70°.

Download Full Size | PDF

4. Phase extraction

Phase distributions φ ₁ and φ ₂ can be obtained by solving the linear Vandermonde system of equations. To design the equations, the number of harmonics κ present in the signal and values of phase steps α and β pixel-wise must be known. Once, these values are known, the parameters ℓ_k and ℘_k can be solved using the linear Vandermonde system of equations obtained from Eq. (2). Phase distributions φ ₁ and φ ₂ are subsequently computed from the arguments of ℓ ₁ and ℘ ₁. The matrix thus obtained can be written as

where (α₀,β ₀), (α₁,β ₁), .., and (α_N-1,β _N-1)are phase steps for frames I ₀, I ₁, I ₂,.., and I _N-1, respectively. The advantage of Vandermonde system of equations is that the matrix shown in Eq. (32) can always be inverted as long as different values of α and β are used in the design of equations. Figures 5(a) and 5(b) show typical errors in the computation of phase values φ ₁ and φ ₂ under the assumption that SNR=30, N=22, and m=15. For obtaining the errors in typical phase maps, the results obtained in Fig. 4(c) using the forward-backward approach is considered. Figure 6 shows wrapped phase maps for φ ₁ and φ ₂, all the other parameters remaining the same.

 

Fig. 5. Plots show typical error in computation of phase distribution a) φ ₁ (in radians), and b) φ ₂ (in radians), for phase step obtained from Fig. 4(c) for 30dB noise.

Download Full Size | PDF

 

Fig. 6. Plot shows wrapped phase for a) φ ₁ and b) φ ₂ for phase step obtained from Fig. 4(c) for 30dB noise.

Download Full Size | PDF

5. Conclusion

We have proposed a novel approach for recovering phase distribution in optical setups accommodating multiple PZT’s. The successful extraction of dual phase distributions from a simulated holographic moiré pattern shows the feasibility of our proposed concept. A significant advantage of the proposed algorithm lies in its ability to measure phase steps in the presence of noise. The algorithm enables to choose arbitrarily phase shifts between 0 and π. The proposed technique works well with both diverging and converging beams since it relies on retrieving the phase step values pixelwise before applying them to the Vandermonde system of equation. The accuracy in the measurement of the phase steps in the presence of additive white Gaussian noise has been shown to increase with increasing number of data frames. Further research will focus on studying the statistical behavior of the algorithm using Cramér-Rao bound analysis in retrieving two phase steps.