1. Introduction

In phase shifting interferometry, phase increments between the object and reference beams are usually introduced by piezo-actuator device (PZT). One of the sources of error associated with the use of PZT is the phase shift calibration error [1]. This causes the user to rely upon the shift value obtained by the calibration as being exact. Non-sinusoidal waveform of the signal arising due to laser cavity or CCD nonlinearity is another probable source of error in the computation of phase. Several phase shifting algorithms have been proposed which overcome these errors [2–4]. However, these algorithms impose conditions on phase step values in order to minimize the errors due to nonsinusoidal waveforms which is assumed to be known a priori [2–4]. Recently a phase shifting algorithm has been proposed which accommodates higher order harmonics and facilitates use of arbitrary phase steps [5]. However, the algorithm requires incorporation of additional denoising procedure in order to extract phase steps in the presence of noise. This procedure introduces an intermediate step which unnecessarily adds to the computational cost of the method. The step could also be a source of potential error.

In this paper an algorithm that can be used to determine the phase of the wavefront in presence of systematic error sources such as phase shift miscalibration and nonsinusoidal waveforms is developed. A salient feature of the algorithm is that it enables for computing in real-time values of phase steps at each pixel point. The algorithm based on MUltiple SIgnal Classification technique (MUSIC) [6–7] thus besides allowing to handle miscalibration error and nonsinusoidal waveforms offers the flexibility of choosing spherical beams and arbitrary phase step values between 0 and π radians. The paper also investigates the influence of white Gaussian noise on the performance of the algorithm [8]. The advantage of the proposed algorithm as compared to the one previously suggested [5] lies in its ability to determine the phase step values pixelwise without the incorporation of a denoising procedure. The proposed method does not only reduce the computational burden but also eliminates a source of potential error in the determination of phase step values.

MUSIC, a subspace based method, has been successfully used to estimate the frequencies of sinusoids corrupted with white noise. Although spectrum estimation can be handled efficiently by Fourier transform, resolution of closely spaced frequencies in the presence of noise is troublesome with limited number of data samples. Drawing a parallelism between the frequencies present in the spectrum and the phase shifts applied to the PZT, the proposed method offers an innovative means of estimating the phase step values at each pixel point in real time. The method is based on a canonical decomposition of a positive definite Toeplitz matrix formed by an estimated covariance sequence. The sequence corresponding to the phase shifted intensity images buried in white noise is acquired temporally. This decomposition yields a very important result as far as the estimation of frequency is concerned. Since the frequencies are estimated using the eigendecomposition of the covariance matrix, the method is referred as subspace based method.

Section 2 outlines the MUSIC algorithm. Section 3 shows simulation results corresponding to extraction of phase step values using forward-backward approach [9] applied in the design of covariance sequence. Section 4 presents a generalized approach by which the phase distribution can be obtained. This section also demonstrates typical errors that may arise in computing the phase pixel-wise in the presence of noise using the MUSIC algorithm.

2. Subspace-based method

The recorded fringe intensity at a point (x, y) on the t^th frame is given by

where, I_dc is the local average value of intensity, a_k is the complex Fourier coefficient, i=√-1, u_k =exp(ikα), superscript ∗ denotes the complex conjugate, η is the white Gaussian noise with mean zero and variance σ ²; and φ, α, and k represent phase distribution, phase step, and the order of harmonics, respectively.

First step consists of forming the covariance matrix from N recorded phase shifted sequences. The covariance of signal I(t) in Eq. (1) is defined as [10]

where, E[·] represents the expectation operator which averages over the ensemble of realizations; the terms $A_0^2$ ,$A_1^2$ ,…,$A_k^2$ ….$A_{2\kappa }^2$ are explained in Appendix A, σ ² is the variance, and δ _p,0 is the Kronecker delta (δ _g,h =1 if g=h ; and δ_g,h =0 otherwise). The reader is referred to Appendix A for the derivation of Eq. (2). The covariance of function I(t) is assumed to depend only on the lag between the two averaged samples. The covariance matrix can thus be written as [6–7, 10]

where, I(t)=[I(t-1),…..,I(t-m)], m is the covariance length, and (·) ^c is the conjugate transpose of a vector or matrix. The covariance matrix R_I can be shown to have the form

where, R _s and R _ε are the signal and noise contributions, A _m ×(2κ+1)=[a(ω₀ ) . . a(ω2κ)] where for instance element a(ω₀ )consists of m×1 matrix with unity entries corresponding to I_dc ; a(ω₁ )=[1 e^jα . . e^jα (m-1)]^T; I is the m×m identity matrix; and P _{(2κ+1)×(2κ+1)} matrix is

Since R _I is positive semidefinite, its eigen values are nonnegative. The eigen values of R _I can be ordered as λ ₁≥λ ₂≥….≥λ _n ≥….λ_m . Let S _m×n =[s₁,s₂,…..s_n ] be the orthonormal eigenvectors associated with λ ₁≥λ ₂≥…..λ_n . The space spanned by {s₁,s₂,…..s_n } is known as signal subspace. The set of orthonormal eigenvectors G _m×(m-n) =[g₁,g₂,…..g _m-n ] associated with eigen values λ _n+1≥λ _n+2≥…..λ_m spans a subspace known as noise subspace. Since APAc∈C ^m×m (C represents complex matrix) has rank n (n<m), it has n eigen values and the remaining m-n eigen values are zero. If we further suppose that (λ, r) is an eigenpair of L∈C ^m×m and W=L+ρI with ρ∈ C, then (λ+ρ, r) is an eigenpair of W, and in consequence we obtain λ_t =λ̃ _t +σ ², where t spans from 1, 2,3,…,m. We observe that λ ₁≥λ ₂≥…..λ _n≥σ ² and λ _n+1=…..=λ_m =σ ². Following this corollary and from Eq. (4) we get

Last equality in Eq. (6) means that APAcG=0 and since AP has full column rank we have A ^c G=0. This means that sinusoidals ${\{\mathbf{a}({\omega }_k)\}}_{k=0}^n$ are orthogonal to noise subspace. This can be stated as (${\displaystyle \sum _{k=m-n}^m{\mathbf{a}}^c(\omega ){\mathbf{G}}_k=0}$. Hence, true frequencies ${\left\{{\omega }_k\right\}}_{k=0}^n$ are the only solutions to the equation a ^T(ω)GG ^c a(ω)=‖G ^c a(ω)‖²=0 for m>n. Here, (·)^T represents transpose of a matrix. Since, in practice, only the estimate R̂ _I of R _I is available, only the estimate Ĝ of G can be determined. In the present study we employ root MUSIC [11] which computes the frequencies as angular positions of n roots of equation

that are nearest and inside the unit circle. Here, a(z) is obtained from a(ω), and by replacing e^jω by z we get a(z)=[1 z ^-1 . z-^(m-1)]^T. Since the minimum possible value of m is n +1, from Eq. (7) it can be observed that we need data frames (N) that are at least twice the number of sinusoidal components in the signal. Hence, if κ=2, which means n=5 (since n=2κ+1), we need at least ten data frames (the dc component is also counted as dc frequency). Hence, the minimum number of data frames required while using MUSIC method for phase extraction is 4κ+2.

To apply MUSIC for estimating the phase step, one first needs to find out the number of harmonics present in the signal so that appropriate value of n is determined. The details on selection of m and N will be explained in next Section.

In many cases κ is unknown and can be determined [12] by observing the Singular Value Decomposition (SVD) of R _I matrix in Eq. (3). For a noiseless signal the SVD of R _I =USV ^T results in a diagonal matrix S with 2κ+1 nonzero and N-2κ-1 zero singular values, where U and V are unitary matrices. If the data is noisy, the M=2κ+1 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 ${\mathbf{R}}_I^{\mathrm{T}}$ R _I are less susceptible to noise perturbations in comparison to the remaining N-M eigenvectors. Figure 1 illustrates typical singular values of S obtained from SVD of matrix R _I for noise, at SNR of 10 dB, and without noise, and when κ=2 in Eq. 1. Although eight frequencies were assumed to be present during the estimation, only five principal values of S (corresponding to I_dc , α, -α, 2α, -2α) for noisy and noiseless signals show a distinctly larger magnitude as compared to the remaining values. The plot thus allows a reliable estimation of the number of harmonics.

 

Fig. 1. Plot of the diagonal of S versus frequency for noiseless signal and signal with SNR=10 dB.

Download Full Size | PDF

3. Evaluation of the algorithm

The concept is tested by simulating the fringe pattern in Eq. (1) with phase step α=π/4, κ=2, and phase φ given by

where x′×x′ is the origin of the fringe pattern. In practice, only the estimate of R _I , represented as R̂ _I , of a covariance matrix is known and the sample covariance matrix is designed using [9]

which is as close as possible to R _I in Eq. (3) in least squares sense. The methods which obtain the frequency estimates from R̂ _I given in Eq. (9) are called forward-backward approaches [9]. We first discuss retrieving the phase step in the presence of additive white Gaussian noise with SNR between 0 and 70 dB. First, the number of frequencies present in the signal must be determined. In the present example, the value of n (number of frequencies) is determined to be n=5 using the method suggested in Section 2. Hence, the minimum number of data frames N for extracting the phase step values is 10 (4κ+2). Subsequently, an appropriate value of m must be selected such that N>m>n. It is observed that m>n increases the accuracy of frequency estimates. This is achieved at a higher computational cost and also m too close to N does not yield the covariance matrix R̂ _I similar to R _I , which in turn results in spurious frequency estimates. Performing eigendecomposition of R̂ _I gives estimates for eigenvectors G and S, represented as Ĝ and Ŝ respectively. Finally using Eq. (7) the frequencies or the phase step values α are estimated pixelwise.

Figures 2(a)–(b) shows the plot for the case when data frame N=10, and m=7 and m=9, respectively. As expected, the phase step α at any arbitrary pixel location on the data frame cannot be estimated at lower SNRs (below 25 dB) from this plot. Figure 2(b) shows that value of m too close to N does not yield result. In the second case the number of data frame N=14 is selected. Figures 2(c)–(e) show typical plot for phase step α with m=7, 9 and 12, respectively. From Figs. 2(c)–2(e), it can be observed that m=9 yield better results as compared to those obtained with m=7 or 12. In the third case N=18 is selected, and plot for m=12 is shown in Fig. 2(f). From the plot it can be observed that phase steps α can be more reliably estimated than that from Fig. 2(d) at lower SNR’s. From these three cases, it can be concluded that phase step α can be reliably estimated even at lower SNR’s with increase in data frames N, and for m not too close to n and N (as a thumb rule, midway between n and N).

 

Fig. 2. Plots of phase step values α (in degrees) obtained using forward-backward approach at an arbitrary pixel point with different values of N and m.

Download Full Size | PDF

4. Phase distribution measurement

Once the phase step values α has been estimated pixelwise, the parameter ℓ _k can be solved using a linear Vandermonde system of equations obtained from Eq. (1), and which can be written as

where, α ₀,…,α_N -1are the phase steps for frames I ₀,..,I _N-1, respectively. The phase φ is computed from the argument of ℓ₁. Figure 3 shows typical phase obtained using Eq. (10). For computing the phase φ, the phase step values α obtained using fourteen data frames is used. Figure 4 shows typical error in the computation of phase φ and for the SNR of 30 dB.

 

Fig. 3. Plot shows typical wrapped phase φ (in radians) for phase step values obtained from Fig. 2(d) for 30dB noise.

Download Full Size | PDF

 

Fig. 4. Plot shows typical absolute error obtained in computation of phase φ (in radians) for phase step values determined from Fig. 2(d) for 30dB noise.

Download Full Size | PDF

5. Conclusion

To conclude, we have proposed a new generalized approach for recovering phase distribution in the presence of higher order harmonics. The method offers the flexibility to select arbitrary phase shifts between 0 and π. The accuracy in the measurement of the phase step in the presence of additive white Gaussian noise has been shown to increase with large data frames. The proposed technique works well with both diverging and converging beams since it first retrieves the phase step values pixel wise before applying them to the Vandermonde system of equation. The advantage of the proposed method also lies in its ability to measure phase step values in real time. Simulated results demonstrate the effectiveness of the proposed method.