## Abstract

A subspace-based method is applied to phase shifting interferometry for obtaining in real time values of phase shifts between data frames at each pixel point. A generalized phase extraction algorithm then allows for computing the phase distribution. The method is applicable to spherical beams and is capable of handling nonsinusoidal waveforms in an effective manner. Numerical simulations demonstrate phase measurement with high accuracy even in the presence of noise.

©2005 Optical Society of America

## 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 *σ*
^{2}; 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, *σ*
^{2} 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**(*ω*_{0}
) . . **a**(*ω2κ*)] where for instance element **a**(*ω*_{0}
)consists of *m*×1 matrix with unity entries corresponding to *I*_{dc}
; **a**(*ω*_{1}
)=[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 *λ*
_{1}≥*λ*
_{2}≥….≥*λ*
_{n}
≥….*λ*_{m}
. Let S
_{m×n}
=[s_{1},s_{2},…..*s*_{n}
] be the orthonormal eigenvectors associated with *λ*
_{1}≥*λ*
_{2}≥…..*λ*_{n}
. The space spanned by {s_{1},s_{2},…..*s*_{n}
} is known as *signal subspace*. The set of orthonormal eigenvectors **G**
_{m×(m-n)}
=[g_{1},g_{2},…..g
_{m-n}
] associated with eigen values *λ*
_{n+1}≥*λ*
_{n+2}≥…..*λ*_{m}
spans a subspace known as *noise subspace*. Since **APA**c∈**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}
+*σ*
^{2}, where *t* spans from 1, 2,3,…,*m*. We observe that *λ*
_{1}≥*λ*
_{2}≥…..*λ*
_{n}≥*σ*
^{2} and *λ*
_{n+1}=…..=*λ*_{m}
=*σ*
^{2}. 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 ($\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**(*ω*)‖^{2}=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.

## 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*).

## 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, *α*
_{0},…,*α*_{N}
-1are the phase steps for frames *I*
_{0},..,*I*
_{N-1}, respectively. The phase *φ* is computed from the argument of ℓ_{1}. 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.

## 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.

## Appendix A

Let us consider a signal

for *t*=0,1,2,…*m*,..,*N*-1

Here, *η (t)* represents white Gaussian noise. The covariance of a function *I (t)* is defined as [10]

For simplicity, let us consider *κ*=1 and rewrite Eq. (A1) as

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

Substituting Eqs. (A3) and A(4) in Eq. (A2), we obtain the following

Equation (A5) can be written in the following compact form

where, *c*
_{1}=*I*_{dc}*a1e-*^{iφ}*e-*_{iαt}
+*I*_{dc}*a1eiφe*^{iαt}
, *c*
_{2}=*I*_{dc}*a1e-iφe-*^{iαt}
${\mathit{+}a}_{1}^{2}$
*e*^{−iφt} *e*^{−iαt}, and *c*_{3}=*I*_{dc}*a1e-iφe*^{iαt}
+${a}_{1}^{2}$
*e*^{2iφ} *e*^{2iαt}.

Let, **E**{${I}_{\mathit{\text{dc}}}^{2}$
+*c*
_{1}}=${A}_{0}^{2}$, *E*{${a}_{1}^{2}$+*c*2}=${A}_{1}^{2}$, and *E*{${a}_{1}^{2}$+*c*
_{3}}=${A}_{2}^{2}$. Therefore,

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

In practice, expectation *E* in Eq. (2) 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*
_{1}, *c*
_{2}, and *c*
_{3}, will oscillate uniformly between 0 and 2*π*. In this limit, the expectation of *c*
_{1}, *c*
_{2}, and *c*
_{3} will approach zero because

However, if finite number of frames are taken for averaging, expectation *c*
_{1}, *c*
_{2}, and *c*
_{3} 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

## Acknowledgements

This research is funded by the Swiss National Science Foundation.

## References and Links

**1. **Y. Surrel, “Phase stepping: a new self-calibrating algorithm,” Appl. Opt. **32**, 3598–3600 (1993). [CrossRef] [PubMed]

**2. **Y. Surrel, “Design of algorithms for phase measurements by the use of phase stepping,” Appl. Opt. **35**, 51–60 (1996). [CrossRef] [PubMed]

**3. **K. Hibino, B. F. Oreb, D. I. Farrant, and K. G. Larkin, “Phase shifting for nonsinusoidal waveforms with phase-shift errors,” J. Opt. Soc. Am. A **12**, 761–768 (1995). [CrossRef]

**4. **K. G. Larkin and B. F. Oreb, “Design and assessment of symmetrical phase-shifting algorithms,” J. Opt. Soc. Am. A **9**, 1740–1748 (1992). [CrossRef]

**5. **A. Patil, R. Langoju, and P Rastogi, “An integral approach to phase shifting interferometry 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] [PubMed]

**6. **R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” in *Proceedings RADC, Spectral Estimation Workshop*, Rome, NY, (243–258) 1979.

**7. **G. Bienvenu, “Influence of the spatial coherence of the background noise on high resolution passive methods,” in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing, Washington, DC, 306–309 (1979)

**8. **C. Rathjen, “Statistical properties of phase-shift algorithms,” J. Opt. Soc. Am. A **12**, 1997–2008 (1995). [CrossRef]

**9. **B. D. Rao and K. V. S. Hari, “Weighted subspace methods and spatial smoothing: analysis and comparison,” IEEE Transactions on Signal Processing **41**, 788–803 (1993). [CrossRef]

**10. **T. Söderström and P. Stoica, “Accuracy of higher-order Yule-Walker methods for frequency estimation of complex sine waves,” IEE Proceedings-F **140**, 71–80 (1993).

**11. **A. J. Barabell, “Improving the resolution performance of eigenstructure-based direction-finding algorithms,” in *Proceedings of the International Conference on Acoustics, Speech, and Signal Processing*, Boston, MA, 336–339 (1983).

**12. **J. J. Fuchs, “Estimating the number of sinusoids in additive white noise,” IEEE Transactions on Acoustics, Speech, and Signal Processing **36**, 1846–1853 (1988) [CrossRef]