## Abstract

The objective of this paper is to describe an integral approach -based on the use of a super-resolution frequency estimation method - to phase shifting interferometry, starting from phase step estimation to phase evaluation at each point on the object surface. Denoising is also taken into consideration for the case of a signal contaminated with white Gaussian noise. The other significant features of the proposal are that it caters to the presence of multiple PZTs in an optical configuration, is capable of determining the harmonic content in the signal and effectively eliminating their influence on measurement, is insensitive to errors arising from PZT miscalibration, is applicable to spherical beams, and is a robust performer even in the presence of white Gaussian intensity noise.

© 2004 Optical Society of America

## 1. Introduction

Phase Shifting has been widely adopted in optical interferometry for retrieving phase information encoded in the interference fringes [1–5]. The technique functions primarily by acquiring a number of intensity images with phase increments between successive frames. These phase increments are generally applied using a piezoelectric device (PZT). Although three frames are normally sufficient to compute the phase distribution, the measurement process is sensitive to various systematic [6–18] and random errors [19–24]. The large number of frames has been shown to decrease the sensitivity to these errors.

This method, commonly known as temporal technique for phase measurement, can be broadly classified to fall into two categories: conventional and generalized. While the step widths are typically multiples of π/2, in the former case, these are arbitrary in the latter case [25–26]. As far as conventional techniques are concerned, several algorithms have been reported to compensate for errors that are introduced during phase evaluation. The systematic errors occurring in the phase shift device have been investigated. The corresponding algorithms insensitive to these errors are described in the literature. Most of these algorithms aim at minimizing the errors arising due to the nonlinear characteristic inherent in the piezoelectric device PZT. These algorithms termed as self-compensating [1,8–10,12,22,27] have proved to be effective while compensating deterministic phase shift error, even until third order nonlinearities. Some of the algorithms minimize the effect of other systematic errors such as optical system aberrations [5,18], parasitic fringes [6,28], photodetection errors [29,30] and quantization errors [17,31]. Random errors due to thermal or shot noise have also been observed to influence the phase measurement [21 ]. Random errors in measurement of phase can also crop up because of turbulent and laminar air flows in the optical path [30], vibrations [19–20,22], electronic noise [32] originating in the photodetection during amplification, and optical noise commonly known as speckle [33]. The statistical properties of phase-shift algorithms have been investigated for the case of additive Gaussian intensity noise. It has been shown that some advanced systematic-error-compensating algorithms are sensitive to random noise [21].

The error in phase measurement can also arise if the fringe profile is not purely sinusoidal because of detector nonlinearity or multiple-beam interference. The algorithm proposed by Hibino et al [14] offers the possibility to reduce the calibration error in the presence of higher order harmonics. Hibino et al [14] and Surrel [11] showed that 2*κ* + 3 and 2*κ* + 2 samples are respectively necessary to minimize the effect of the *κ*
* ^{th}* order harmonic, with the phase step 2

*π*/(

*κ*+2) between samples. The restriction on the free use of phase steps represents a drawback and a possible limitation for the methods wide spread application. The other important limitations of the algorithm reside in its inability to function in configurations requiring multiple PZTs.

While conventional phase shifting has demonstrated its ability to overcome various sources of systematic errors, generalized phase shifting technique though flexible as far as the selection of phase steps is concerned hasn’t been successful in overcoming most of the systematic errors discussed earlier. Although the algorithm proposed by Carré [1] is insensitive to first order calibration error and enables the use of arbitrary phase steps, the study imposes a restriction on the choice of the phase step value in order to minimize certain other errors referred in [34]. The algorithm proposed by Lai and Yatagai [35] requires an additional optical setup for the generation of Fizeau fringes and can only handle fringes that are not straight and equally spaced [36]. Farrell and Player [37–38] describe an ellipse fitting technique that uses Bookstein algorithm [39] - a forced fit method - to fit the data set into an ellipse. The data set fits into a hyperbole for noisy data [40]. The statistical approach proposed by Cai et al [41] to extract reference phase steps imposes requirements on the spatial resolution of the CCD camera (which in turn dictates the sampling number of the speckle intensity). The max-min algorithm proposed by Chen et al is computationally exhaustive and requires a large number of data frames (15 or more) for reliable operation [42]. The algorithm based on iterative least squares estimation proposed by Han and Kim is suitable only for compensating linear first order errors in PZT and requires the initial guess of the phase step to be as close as possible to its true value [43]. The concept of identifying the reference phase steps proposed by Jüptner et al is suitable for only linear miscalibration error and can produce large number of outliers for noisy interferograms [44]. Although the five frame algorithm proposed by Stoilov and Dragostinov gives the flexibility of using arbitrary phase steps, its use is limited only to linear errors encountered in the PZT [45]. The stochastic approach [46] to detect the PZT nonlinearity in PZT proposed by Patil et al [47] is sensitive to random sources of error.

To the best of authors knowledge no generalized phase shifting technique has been reported so far which at a time caters to the presence of multiple PZTs in an optical configuration, is insensitive to errors due to nonsinusoidal waveforms, to PZT miscalibration and to additive white Gaussian intensity noise, and has the ability to work with diverging as well as converging beam. The use of multiple PZTs has been introduced in an optical configuration by Rastogi [48–49] for automating the extraction of information conveyed by both carrier and moiré in a holographic moiré setup. With the field of phase shifting interferometry advancing rapidly, the possibility of incorporating multiple PZTs in an interferometric configuration is a significant conceptual progress which should be useful for developing multiple channel measurement capabilities in moiré based interferometers.

The aim of this paper is to propose an integral generalized approach which first caters to the presence of multiple PZTs in an optical configuration, second is capable of handling harmonics, third is insensitive to errors arising from PZT miscalibration, fourth is applicable to spherical beams, and fifth is a robust performer even in the presence of white Gaussian intensity noise. We intend to apply the super-resolution frequency estimation approach which basically takes the advantage of the fact that a polynomial can be associated to a periodic intensity fringe pattern recorded by a CCD. Another polynomial termed as *annihilating filter* [50] is designed which has zeroes at the frequencies that are associated with the polynomial for the intensity fringes. The discrete convolution of the filter and the intensity fringe yields zero. The spectral information embedded in the signal corresponds to the phase steps imparted by the multiple PZTs in the interferometric setup. Hence, the parametric estimation of this annihilating filter yields the desired spectral information embedded in the signal, which in our case are the phase steps. Although the discrete Fourier transform is an efficient tool for the estimation of well separated frequencies, the separation of closely spaced frequencies in the presence of noise and less number of samples can be handled efficiently with high resolution techniques [51]. This method for frequency estimation is also sometimes referred to as super-resolution technique because of its ability to resolve spectral lines (frequencies) separated in frequency *f* = *ω*/2*π* by less than 1/*N* cycles per sampling interval (here, *ω* and *N* refers to angular frequency and number of samples, respectively), which is the resolution limit for the classical periodogram-based methods, such as Fourier transform. As a case study, we show the effectiveness of our algorithm to extract the dual phase steps in holographic moiré. Since, the proposed technique functions by extracting the phase step at each pixel location of the acquired frames, this method allows the use of diverging as well as converging beams. The advantage of using any arbitrary phase step between 0 and *π* radians overcomes the limitation exhibited by previously reported methods. We also show that the robustness of our algorithm to additive white Gaussian noise can be enhanced by incorporating a denoising step using the concept of truncated singular value decomposition [52–53].

In what follows, we first present a brief introduction to holographic moiré. Section 3 presents the design of annihilating filter for phase step determination. Section 4 deals with the identification of the harmonic content of the waveform. Section 5 shows the denoising procedure. Section 6 shows the simulation results for phase step extraction and retrieval of phases in holographic moiré in the presence of noise.

## 2. Holographic moiré

Although holographic interferometry has proved its effectiveness for measuring displacements along the line of sight, it does not provide means for obtaining whole field distribution of displacements orthogonal to the line of sight. This limitation has led to the development of more complex interferometers to obtain direct visualization of displacement components and their derivatives. The method referred to as holographic moiré has proved its usefulness in determining in-plane displacements [54], slopes and curvature [55], and difference displacements [56] of objects.

The application of phase-stepping to holographic moiré is a relatively complex procedure in which a dual-phase stepping is applied in the two arms of a holographic moiré setup [48–49]. Besides providing automation to the method, one of the significant advances of holographic moiré is that it offers retrieving simultaneously the phase information carried on moiré and carrier at any point on the object surface. The recorded fringe intensity at a point (*x*,*y*) for *m ^{th}* phase step is given by

$$\sum _{k=-1}^{-\kappa}{b}_{k}\mathrm{exp}\left[\mathit{ik}\left({\phi}_{2}+\mathit{m\beta}\right)\right]+\sum _{k=1}^{\kappa}{b}_{k}\mathrm{exp}\left[\mathit{ik}\left({\phi}_{2}+\mathit{m\beta}\right)\right];\phantom{\rule{.2em}{0ex}}\mathrm{for}\phantom{\rule{.2em}{0ex}}m=\mathrm{1,2},\dots \dots ,N$$

where, *a _{k}* and

*b*is the complex Fourier coefficient of the

_{k}*k*order harmonic, $i=\sqrt{-1}$ and

^{th}*I*is the local average value for intensity; and pairs

_{dc}*φ*

_{1}and

*φ*

_{2}, and,

*α*and

*β*, represent phase differences and phase shifts, respectively, in the two arms of the holographic moiré setup. The coefficients

*a*and

_{k}*b*are in fact real and are related to the appropriate choice of phase origin at a point where the intensity reaches a maximum.

_{k}Phase shifted holographic moiré offers the advantage of extracting simultaneously the phase values *φ*
_{1} and *φ*
_{2}. Phase information corresponding to moiré, Φ, and carrier, Ψ, are consequently obtained by and which, for example, in a holographic moiré configuration are directly related to in-plane and out-of-plane displacement components, respectively.

## 3. Theory of estimation of phase steps using a super-resolution method

In order to explain the basis of the spectral estimation procedure, let us write Eq. (1) as

$$\phantom{\rule{10.2em}{0ex}}\mathrm{for}\phantom{\rule{.2em}{0ex}}m=n+1=\mathrm{1,2},\dots \dots ,N$$

where, ℓ* _{k}* =

*a*exp(

_{k}*ikφ*

_{1}),

*u*= exp(

_{k}*ikα*),

*℘*=

_{k}*b*exp(

_{k}*ikφ*

_{2}),

*v*= exp(

_{k}*ikβ*); superscript * denotes the complex conjugate;

*φ*

_{1}and

*φ*

_{2}, and,

*α*and

*β*, are the unknown parameters;

*I*(

_{n}*x*,

*y*;

*m*) corresponds to the

*n*frame, where

^{th}*n*= 0 refers to the data frame corresponding to the first phase shifted intensity pattern. Equation (4) represents the complex valued sinusoidal signals where

*α*,….,

*κα*; -

*α*,….,-

*κα*;

*β*,….,

*κβ*; -

*β*,….,-

*κβ*represent the frequencies embedded in the signal, and by estimating them, the phase steps

*α*and

*β*, imparted by the PZTs can be determined. Frequency estimation of sinusoids is a classical problem in spectral estimation and as mentioned in Section 1 can be better handled using a super-resolution method. In this method we first transform the discrete time domain signal

*I*in Eq. (4) into a complex frequency domain by taking its Z-transform. Let the Z-transform of

_{n}*I*be denoted by I(

_{n}*z*). The objective here is to design another polynomial P(

*z*) termed as

*annihilating filter*which has zeros at frequencies associated with I(

*z*), which in turn would result in I(

*z*)P(

*z*) = 0. Since the frequencies translate into zeros, spectral estimates feature high resolution. In what follows we will derive the expression for P(

*z*) and explore the condition for which the multiplication of P(

*z*) and I(

*z*) is zero.

The Z-transform of Eq. (4) is written as

$$\sum _{n=0}^{N-1}\sum _{k=1}^{\kappa}{\wp}_{k}{v}_{k}^{\left(n+1\right)}{z}^{-n}+\sum _{n=0}^{N-1}\sum _{k=1}^{\kappa}{\wp}_{k}^{*}{\left({v}_{k}^{*}\right)}^{\left(n+1\right)}{z}^{-n}$$

which after expansion and simplification can be written as

$$\phantom{\rule{2.5em}{0ex}}\sum _{k=1}^{\kappa}{\ell}_{k}^{*}{e}^{-\mathit{i\alpha k}}\left(1-{e}^{\mathit{-i\alpha kN}}{z}^{-N}\right)\left[\frac{{D}_{2}\left(z\right)}{{P}_{2}\left(z\right)}\right]+$$

$$\phantom{\rule{2.5em}{0ex}}\sum _{k=1}^{\kappa}{\wp}_{k}{e}^{\mathit{i\beta k}}\left(1-{e}^{\mathit{i\beta kN}}{z}^{-N}\right)\left[\frac{{D}_{3}\left(z\right)}{{P}_{3}\left(z\right)}\right]+$$

$$\phantom{\rule{2.5em}{0ex}}\sum _{k=1}^{\kappa}{\wp}_{k}^{*}{e}^{-\mathit{i\beta k}}\left(1-{e}^{\mathit{i\beta kN}}{z}^{-N}\right)\left[\frac{{D}_{4}\left(z\right)}{{P}_{4}\left(z\right)}\right]$$

where,

$${P}_{1}\left(z\right)=\prod _{k=1}^{\kappa}\left(1-{e}^{\mathit{i\alpha k}}{z}^{-1}\right),\phantom{\rule{.2em}{0ex}}{P}_{2}\left(z\right)=\prod _{k=1}^{\kappa}\left(1-{e}^{-\mathit{i\alpha k}}{z}^{-1}\right),\phantom{\rule{.2em}{0ex}}{P}_{3}\left(z\right)=\prod _{k=1}^{\kappa}\left(1-{e}^{\mathit{i\beta k}}{z}^{-1}\right)$$

$$\mathrm{and}\phantom{\rule{.2em}{0ex}}{P}_{4}\left(z\right)=\prod _{k=1}^{\kappa}\left(1-{e}^{-\mathit{i\beta k}}{z}^{-1}\right).$$

Taking into account the terms D’ s and P’ s and substituting them in Eq. (6), we obtain

$$\phantom{\rule{2.5em}{0ex}}\left\{\frac{{\displaystyle \sum _{k=1}^{\kappa}}{\ell}_{k}^{*}{e}^{-\mathit{i\alpha k}}\left[{\displaystyle \sum _{j=0}^{\kappa -1}}{d}_{2j}{z}^{-j}-{\displaystyle \sum _{j=0}^{\kappa -1}}{d}_{2j}{e}^{-\mathit{i\alpha kN}}{z}^{-\left(j+N\right)}\right]}{{\displaystyle \prod _{k=1}^{\kappa}}\left(1-{e}^{-\mathit{i\alpha k}}{z}^{-1}\right)}\right\}+$$

$$\phantom{\rule{2.5em}{0ex}}\left\{\frac{{\displaystyle \sum _{k=1}^{\kappa}}{\wp}_{k}{e}^{\mathit{i\beta k}}\left[{\displaystyle \sum _{j=0}^{\kappa -1}}{d}_{3j}{z}^{-j}-{\displaystyle \sum _{j=0}^{\kappa -1}}{d}_{3j}{e}^{-\mathit{i\beta kN}}{z}^{-\left(j+N\right)}\right]}{{\displaystyle \prod _{k=1}^{\kappa}}\left(1-{e}^{\mathit{i\beta k}}{z}^{-1}\right)}\right\}+$$

$$\phantom{\rule{2.5em}{0ex}}\left\{\frac{{\displaystyle \sum _{k=1}^{\kappa}}{\wp}_{k}^{*}{e}^{-\mathit{i\beta k}}\left[{\displaystyle \sum _{j=0}^{\kappa -1}}{d}_{4j}{z}^{-j}-{\displaystyle \sum _{j=0}^{\kappa -1}}{d}_{4j}{e}^{-\mathit{i\beta kN}}{z}^{-\left(j+N\right)}\right]}{{\displaystyle \prod _{k=1}^{\kappa}}\left(1-{e}^{-\mathit{i\beta k}}{z}^{-1}\right)}\right\}$$

Equation (7) can be represented as

where, D(*z*) represents the numerator and P(*z*) the denominator in the right hand side of Eq. (7), with the latter being given by

$$\phantom{\rule{1.6em}{0ex}}={\displaystyle \sum _{k=0}^{4\kappa +1}}{p}_{k}{z}^{-k}$$

Given the complex form of D(*z*), let us consider for the sake of simplicity writing only its first term,

$$\phantom{\rule{2.2em}{0ex}}={\displaystyle \sum _{j=0}^{4\kappa}}C{1}_{j}{z}^{-j}-{\displaystyle \sum _{j=0}^{4\kappa}}C{1}_{j+N}{z}^{-\left(j+N\right)}$$

From Eq. (10) we observe that some coefficients of C_{1}(*z*), specifically, *C*1_{0}, *C*1_{1}, *C*1_{2}……*C*1_{4κ}, and *C*1* _{N}*,

*C*1

_{N+1},

*C*1

_{N+2}……

*C*1

_{N+4κ}depend upon the unknown amplitude and frequency. We observe that the coefficients

*C*1

_{4κ+1},

*C*1

_{4κ+2},

*C*1

_{4κ+3}……,

*C*1

_{N-1}are identically zero for

*N*≥ 4

*κ*+ 1. This condition is at the essence of the design of annihilating filter P(

*z*). In what follows we show that the (4

*κ*+1)

*degree polynomial P(*

^{th}*z*), when discretely convolved with I(

*z*) yields zero.

Since the multiplication of two signals in frequency domain corresponds to its discrete convolution in time domain, inverse Z-transform of Eq. (8) gives

where, ⊗ represents the convolution operator. In Eq. (11), *p _{n}* is the inverse Z-transform of P(

*z*) defined by

where, *δ*(*n*) is the unit impulse signal. Equation (11) can be further written as

Using the fact,

$$\mathrm{and}\phantom{\rule{.2em}{0ex}}\delta \left(n\right)=\{\begin{array}{c}0\\ 1\end{array}\begin{array}{c}\phantom{\rule{5.5em}{0ex}}\mathrm{for}\phantom{\rule{.2em}{0ex}}n=0\\ \phantom{\rule{5.2em}{0ex}}\mathrm{else}\phantom{\rule{.2em}{0ex}}\mathrm{where}\end{array}$$

we can express Eq. (13) in a matrix form as

We partition the first matrix on the left hand side of Eq. (15), select the middle *N*-4*κ*-2 rows equations corresponding to zero row values in the matrix on the right hand side, and form a new matrix as follows:

We can infer from Eq. (16) that the polynomial P(*z*) is an annihilating filter which when convolved with the moiré intensity signal yields zero. It can further be deduced that at least *N* ≥ 8*κ*+2 samples are required to extract the roots of the polynomial P(*z*). This enables us to find the unknown values *u _{k}* and

*v*. The phase steps

_{k}*α*and

*β*, can hence be computed using

*α*= ℜ(ln

*u*

_{1}/

*i*) and

*β*= ℜ(ln

*v*

_{1}/

*i*).

## 4. Detection of nonsinusoidal waveform and the corresponding harmonic content

The occurrence of nonsinusoidal waveform - a consequence of detector nonlinearity or multiple reflections in laser cavity - results in an error in the computation of phase. This section explains how to deduce the number of harmonics *κ* present in the signal, which can subsequently be applied in the design of the Vandermonde system of equations in Section 6 for the determination of phase values *φ*
_{1} and *φ*
_{2}. Let us rewrite Eq. (1), for the case *κ* = 1 (i.e. pure sinusoidal wave) and noiseless signal, as

$$\phantom{\rule{5.2em}{0ex}}\underset{T2}{\underbrace{{b}_{-1}\mathrm{exp}\left[-i\left({\phi}_{2}+\mathit{m\beta}\right)\right]+{b}_{1}\mathrm{exp}\left[i\left({\phi}_{2}+\mathit{m\beta}\right)\right]};\phantom{\rule{.2em}{0ex}}\mathrm{for}\phantom{\rule{.2em}{0ex}}m}=\mathrm{1,2},\dots \dots ,N.$$

Initially we assume T1≠0 and T2 = 0, a case in which the intensity *I*(*x*,*y*;*m*) in Eq. (17) corresponds to two wave interferometry. Fourier transform of *I*(*x*,*y*;*m*) for *m* = 1,2,……,*N* should result in three peaks in frequency domain corresponding to *I _{dc}*,

*α*and -

*α*. Similarly for T1≠0 and T2≠0, a case in which the intensity

*I*(

*x*,

*y*;

**m**) in Eq. (17) corresponds to holographic moiré interferometry, Fourier transform of

*I*(

*x*,

*y*;

*m*) for

*m*= 1,2,……,

*N*should result in five peaks in frequency domain corresponding to

*I*,

_{dc}*α*, -

*α*,

*β*, and -

*β*. However, with limited data samples, the resolution of closely spaced frequencies is troublesome in the presence of noise. Also because of the “leakage” phenomenon, energy in the main lobe can leak into the side lobes obscuring and distorting other spectral responses.

Interestingly, the singular value decomposition of **R** in Eq. (16), represented as **R** = **USV**
^{T}, yields more reliable information than the Fourier transform method regarding the number of harmonics present in the signal, where **S** is a diagonal matrix with *M* nonzero and *N* - *M* zero singular values; **U** and **V** are unitary matrices [57]. When T1 ≠ 0 and T2 = 0, in Eq. (17), and for a noiseless case, the number of nonzero diagonal entries in **S** is *M* = 3 (corresponding to *I _{dc}*,

*α*, and -

*α*). On the other hand, for T1 ≠ 0 and T2 ≠ 0, the number of nonzero diagonal entries is

*M*= 5 (corresponding to

*I*,

_{dc}*α*, -

*α*,

*β*, and -

*β*). Hence, if H is the number of PZTs used in the optical setup, the number of harmonics

*κ*in the signal can be determined using

*M*= 2H

*κ*+1. In the presence of noise,

*M*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

**R**

^{T}

**R**are less susceptible to noise perturbations in comparison to the remaining

*N*-

*M*eigenvectors.

Figure 1 illustrates typical singular values for **S** when the number of samples (data frames) used are *N* = 25. In this example we assume the presence of second order harmonics in the signal, *κ* = 2, and the presence of two PZTs in the optical configuration, H = 2. In this case and for a noiseless signal we should expect *M* = 2H*κ* + 1 = 9 diagonal entries of **S** in Fig. 1 to be significantly larger in magnitude; the *N* - 9 diagonal entries being zero. Even in presence of noisy data for the case with (SNR = 10 dB), *M* = 9 principal values are still larger in magnitude than the *N* - 9 diagonal values. Hence by selecting only those significant values of **S** (which in this example, corresponds to *M* = 9), the number of harmonics *κ* is determined to be *κ* = 2. The plot thus allows for a reliable estimation of the number of harmonics present in the signal.

## 5. Denoising the signal

There are various sources of noise in optical interferometry and the relative contribution of each noise source depends on the particular system and its application [58]. Denoising is thus an important step towards reducing the effect of noise. We consider white Gaussian noise to contaminate the signal; the reason being that many distributions can be approximated by a Gaussian function. It should also be noted that many noise sources are known to follow a Gaussian noise source in accordance to the central limit theorem [21]. To enhance the signal in the presence of white Gaussian noise, we apply the concept of truncated singular value decomposition. In the present case, number of data frames ≥ 8*κ*+2 are needed for applying the denoising procedure.

The concept of TSVD involves the following steps [52–53]. The matrix **R** is first written in Hankel matrix form, say **R̂**, and its singular value decomposition **R̂** = **ÛS̄V̂**
^{T} shows the nonsingular principal values of **S̄** which are significantly different from zero. After setting the non significant *N*-*M* singular values of **S̄** to zero, a matrix **Ŝ** is formed. A denoised matrix **Z**
* _{M}* =

**ÛŜV̂**

^{T}which approximates

**R̂**in the least squares sense, is then obtained by using the first

*M*columns of

**Û**,

**Ŝ**and

**V̂**;

**Û**and

**V̂**being unitary matrices. Finally, a denoised signal

*Ī*is retrieved by arithmetic averaging along the anti-diagonals (or diagonals) of

_{n}**Z**

*using*

_{M}where, *r* = max[1, *n* - number of rows (**R̂**)+1) and *q* = min[number of rows (**R̂**), *n*]. The denoised signal *Ī _{n}* is subsequently applied in Eq. (16).

## 6. Evaluation of phase steps and phase distributions in presence of noise

The next step after denoising the signal is to extract the phase steps. We present the application of super-resolution technique for the extraction of phase steps by simulating moiré fringes given by

$$\phantom{\rule{5.6em}{0ex}}{\displaystyle \sum _{k=-2}^{2}}{b}_{k}\mathrm{exp}\left(i\left\{\frac{2\pi \left[{\left(x-{p}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}\right]}{{\lambda}_{2}\xi}+\mathit{m\beta}\right\}\right)$$

$$\phantom{\rule{15.em}{0ex}};\phantom{\rule{.2em}{0ex}}\mathrm{for}\phantom{\rule{.2em}{0ex}}m=n+1=\mathrm{1,2},\dots \dots N$$

where, (*x*
_{0},*y*
_{0}) is the origin of the *X*×*Y* pixels for fringe pattern with pitch *λ*
_{1}; (*p*
_{0},*y*
_{0}) is the origin of *X*×*Y* pixels for fringe pattern with pitch *λ*
_{2}; *ξ* is some arbitrary constant; *κ* = 2 and the phase steps are selected as *α* = *π*/4 and *β* = *π*/3. A white Gaussian noise with SNR from 0 to 100 dB is added to test the robustness of the proposed concept. Fringes shown in Fig. 2(a)–(c) correspond to Eq. (19) for *κ* = 1 and for three different noise levels. On the other hand, fringes shown in Figs. 3(a)–(c) corresponds to Eq. (19) for second order harmonic, *κ* = 2, but for same noise levels as in Figs. 2(a)–(c), respectively. These fringes have been generated under the assumption *a*
_{0, ± 1, ± 2} = *b*
_{0, ± 1, ± 2} = 0.5.

Figure 4(a) shows the plot for phase step extraction when the noise level is considered from 0 to 100 dB; eighteen data frames are used to obtain this plot. As can be seen from the plot, the phase steps cannot be accurately determined for SNR level below 40 dB. In order to extract the phase steps reliably at a lower SNR, additional data frames are mandatory so that denoising procedure explained in Section 5 can be successfully applied. Figures 4(b) and 4(c) show a plot for the phase step extraction when twenty seven and thirty six data frames are used in the computation of phase steps and denoising procedure is applied. The plot shows an increase in reliability in the computation of phase steps with a larger number of frames. Large data frames can induce unwanted noise in the measurement, if controlled environment is not used. Our denoising technique works well when random noise follows a Gaussian profile, which is the case most of time. For noise which does not follow a Gaussian profile (air turbulence, or vibrations), reference 30 provides an insight into the measures which can be taken to minimize their influence on the measurement of phase distribution.

Once the phase steps are estimated the parameters ℓ* _{k}* and

*℘*can be solved using the linear Vandermonde system of equations obtained from Eq. (4). The Vandermonde system of equations always has full rank as long as

_{k}*α*’

_{i}*s*and

*β*’

_{i}*s*are distinct. The matrix thus obtained can be written as

where (*α*
_{1},*β*
_{1}), (*α*
_{2},*β*
_{2}), (*α*
_{3},*β*
_{3}), .. and (*α _{N}*,

*β*) are the phase steps for frames

_{N}*I*

_{0},

*I*

_{1},

*I*

_{2}, .., and

*I*

_{N-1}, respectively. The phases

*φ*

_{1}and

*φ*

_{2}are subsequently computed from the argument of ℓ

_{1}and

*℘*

_{1}. Figures 5(a) and 5(b) show typical errors in the computation of phase

*φ*

_{1}without and with the application of the denoising procedure, respectively.

The noise level is assumed to be 30 dB. In this simulation phase steps obtained in Fig. 4(c) are used for the computation of phase. Similarly Figs. 6(a) and 6(b) show errors in the computation of phase *φ*
_{2} without and with the application of the denoising step, respectively, all the other parameters remaining the same. Figure 7 shows the wrapped phases *φ*
_{1} and *φ*
_{2} obtained with the denoising step.

## 7. Conclusion

To conclude, we have presented a novel integral approach - based on the use of a super-resolution frequency estimation method - to phase shifting interferometry, starting from phase step estimation to phase evaluation at each point on the object surface. A significant advantage of this approach is that it also offers the possibility to identify the harmonic content of the signal. The described method holds promise for determining simultaneously multiple phase distributions in the presence of higher order harmonics and in optical configurations with multiple PZTs. The proposed technique works well with the diverging as well as converging beams since the phase steps are retrieved point by point before being applied to the Vandermonde system of equation. The robustness of the algorithm to white Gaussian noise is enhanced by the incorporation of denoising step before the extraction of phase step. Further research will focus on estimating the lower bounds of the distance between phase steps that can be resolved using the high resolution method.

## Acknowledgments

This research is funded by the Swiss National Science Foundation.

## References and Links

**1. **P. Carré, “Installation et utilisation du comparateur photoélectrique et interférentiel du Bureau International des Poids et Mesures,” Metrologia **2**, 13–23 (1966). [CrossRef]

**2. **J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White, and D. J. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses, App. Opt. **13**, 2693–2703 (1974). [CrossRef]

**3. **K. Creath, “Phase-shifting holographic interferometry,” *Holographic Interferometry*, P. K. Rastogi, ed. (Springer Series in Optical Sciences, Berlin1994), Vol.68, pp. 109–150.

**4. **T. Kreis, “Holographic interferometry Principles and Methods,” Akademie Verlag, 1996, pp. 101–170.

**5. **J. E. Greivenkamp and J. H. Bruning, Phase shifting interferometry *Optical Shop Testing* ed D. Malacara (New York: Wiley) 501–598 (1992).

**6. **J. Schwider, R. Burow, K. E. Elssner, J. Grzanna, R. Spolaczyk, and K. Merkel, “Digital wave-front measuring interferometry: some systematic error sources,” Appl. Opt. **22**, 3421–3432 (1983). [CrossRef] [PubMed]

**7. **Y. Zhu and T. Gemma, “Method for designing error-compensating phase-calculation algorithms for phase-shifting interferometry,” App. Opt. **40**, 4540–4546 (2001). [CrossRef]

**8. **P. Hariharan, B. F. Oreb, and T. Eiju, “Digital phase-shifting interferometry: a simple error-compensating phase calculation algorithm,” App. Opt. **26**, 2504–2506 (1987). [CrossRef]

**9. **J. Schwider, O. Falkenstorfer, H. Schreiber, and A. Zoller, “New compensating four-phase algorithm for phase-shift interferometry,” Opt. Eng. **32**, 1883–1885 (1993). [CrossRef]

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

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

**12. **J. van Wingerden, H. J. Frankena, and C. Smorenburg, “Linear approximation for measurement errors in phase shifting interferometry,” Appl. Opt. **30**, 2718–2729 (1991). [CrossRef] [PubMed]

**13. **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]

**14. **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]

**15. **Y. -Y. Cheng and J. C. Wyant, “Phase-shifter calibration in phase-shifting interferometry,” App. Opt. **24**, 3049–3052 (1985). [CrossRef]

**16. **B. Zhao, “A statistical method for fringe intensity-correlated error in phase-shifting measurement: the effect of quantization error on the N-bucket algorithm,” Meas. Sci. Technol. **8**, 147–153 (1997). [CrossRef]

**17. **B. Zhao and Y. Surrel, “Effect of quantization error on the computed phase of phase-shifting measurements,” Appl. Opt. **36**, 2070–2075 (1997). [CrossRef] [PubMed]

**18. **R. Józwicki, M. Kujawinska, and M. Salbut, “New contra old wavefront measurement concepts for interferometric optical testing, Opt. Engg. **31**, 422–433 (1992). [CrossRef]

**19. **P. de Groot, “Vibration in phase-shifting interferometry,” J. Opt. Soc. Am. A **12**, 354–365 (1995). [CrossRef]

**20. **P. de Groot and L. L. deck, “Numerical simulations of vibration in phase-shifting interferometry,“ App. Opt. **35**, 2172–2178 (1996). [CrossRef]

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

**22. **P. L. Wizinowich, “Phase-shifting interferometry in the presence of vibration: a new algorithm and system,” Appl. Opt. **29**, 3271–3279 (1990). [CrossRef] [PubMed]

**23. **C. Joenathan and B. M. Khorana, “Phase measurement by differentiating interferometric fringes,” J. Mod. Opt. **39**, 2075–2087 (1992). [CrossRef]

**24. **K. Kinnnstaetter, A. W. Lohmann, J Schwider, and N. Streibl, “Accuracy of phase shifting interferometry,” Appl. Opt. **27**, 5082–5089 (1988). [CrossRef]

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

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

**27. **Y. Ishii and R. Onodera, “Phase-extraction algorithm in laser diode phase-shifting interferometry,” Opt. Lett. **20**, 1883–1885 (1995). [CrossRef] [PubMed]

**28. **Ch. Ai and J. C. Wyant, “Effect of spurious reflection on phase shift interferometry,” Appl. Opt. **27**, 3039–3045 (1988). [CrossRef] [PubMed]

**29. **J. Schmit and K. Creath, “Extended averaging technique for derivation of error-compensating algorithms in phase-shifting interferometry,” Appl. Opt. **34**, 3610–3619 (1995). [CrossRef] [PubMed]

**30. **B. V. Dorrío and J. L. Fernández, “Phase-evaluation methods in whole-field optical measurement techniques,” Meas. Sci. Technol. **10**, R33–R55 (1999). [CrossRef]

**31. **C. P. Brophy, “Effect of intensity error correlation on the computed phase of phase-shifting interferometry,” J. Opt. Soc. Am. A **7**, 537–541 (1990). [CrossRef]

**32. **S. Ellingsrud and G. O. Rosvold, “Analysis of data-based TV-holography system used to measure small vibration amplitudes,” J. Opt. Soc. Am. A **9**, 237–251 (1992). [CrossRef]

**33. **G. T. Reid, “Automatic fringe pattern analysis: a review,” Opt. and Lasers Engg. **7**, 37–68 (1986). [CrossRef]

**34. **Q. Kemao, S. Fangjun, and W. Xiaoping, “Determination of the best phase step of the Carré algorithm in phase shifting interferometry,” Meas. Sci. Technol. **11**, 1220–1223 (2000). [CrossRef]

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

**36. **K. G. Larkin, “A self-calibrating phase-shifting algorithm based on the natural demodulation of two-dimensional fringe patterns,” Opt. Exp. **9**, 236–253 (2001). [CrossRef]

**37. **C. T. Farrell and M. A. Player, “Phase step measurement and variable step algorithms in phase-shifting interferometry,” Meas. Sci. Technol. **3**, 953–958 (1992). [CrossRef]

**38. **C. T. Farrell and M. A. Player, “Phase-step insensitive algorithms for phase-shifting interferometry,” Meas. Sci. Technol. **5**, 648–652 (1994). [CrossRef]

**39. **F. L. Bookstein, “Fitting conic sections to scattered data,” Com. Graphics and Image Process. **9**, 56–71 (1979). [CrossRef]

**40. **A. Fitzgibbon, M. Pilu, and R. B. Fisher, “Direct least square fitting of ellipses,” Pattern Anal. and Machine Intelligence. **21**, 474–480 (1999).

**41. **L. Z. Cai, Q. Liu, and X. L. Yang, “Phase-shift extraction and wave-front reconstruction in phase-shifting interferometry with arbitrary phase steps,” Opt. Let. **28**, 1808–1810 (2003). [CrossRef]

**42. **X Chen, M. Gramaglia, and J. A. Yeazell, “Phase-shifting interferometry with uncalibrated phase shifts,” Appl. Opt. **39**, 585–591 (2000). [CrossRef]

**43. **G. -S. Han and S. -W. Kim, “Numerical correction of reference phases in phase-shifting interferometry by iterative least-squares fitting,” Appl. Opt. **33**, 7321–7325 (1994). [CrossRef] [PubMed]

**44. **W. Jüptner, T. Kreis, and H. Kreitlow, “Automatic evaluation of holographic interferograms by reference beam phase shifting,” In W. F. Fagan, ed., *Industrial Applications of Laser Technology*, Proc. Of Soc. Photo-Opt. Instr. Eng. , **1553**, 569–582 (1991).

**45. **G. Stoilov and T. Dragostinov, “Phase-stepping interferometry: five-frame algorithm with an arbitrary step,” Opt. and Lasers in Eng. **28**, 61–69 (1997). [CrossRef]

**46. **B. Raphael and I. F. C. Smith, “A direct stochastic algorithm for global search,” App. Mathematics & Computation **146**/2–3, 729–758 (2003). [CrossRef]

**47. **A. Patil, B. Raphael, and P. Rastogi, “Generalized phase-shifting interferometry by use of a direct stochastic algorithm for global search,” Opt. Lett. **29**, 1381–1383 (2004). [CrossRef] [PubMed]

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

**49. **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é,” App. Opt. **32**, 3669–3675 (1993). [CrossRef]

**50. **P. Stoica and R. Moses, Introduction to Spectral Analysis (Prentice Hall, New Jersey, 1997).

**51. **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]

**52. **M. Dendrinos, S. Bakamidis, and G. Carayannis, “Speech enhancement from noise: A regenerative approach,” Speech Communication **10**, 45–57 (1991). [CrossRef]

**53. **P. C. Hansen and S. H. Jensen, “FIR filter representations of reduced-rank noise reduction,” IEEE Transactions on Signal Processing **46**, 1737–1741 (1998). [CrossRef]

**54. **P. K. Rastogi and E. Denarié, “Visualization of in-plane displacement fields using phase shifting holographic moiré: application to crack detection and propagation,” Appl. Opt. **31**, 2402–2404 (1992). [CrossRef] [PubMed]

**55. **P. K. Rastogi, M. Barillot, and G. Kaufmann, “Comparative phase shifting holographic interferometry,” Appl. Opt. **30**, 722–728 (1991). [CrossRef] [PubMed]

**56. **P. K. Rastogi, “Visualization and measurement of slope and curvature fields using holographic interferometry: an application to flaw detection,” J. Mod. Opt. **38**, 1251–1263 (1991). [CrossRef]

**57. **R. Kumaresan and D. W. Tufts, “Estimating the angles of arrival of multiple plane waves,” IEEE Transactions on Aerospace and Electronic Systems **AES-19**, 134–139 (1983). [CrossRef]

**58. **K. B. Hill, S. A. Basinger, R. A. Stack, and D. J. Brady, “Noise and information in interferometric cross correlators,” Appl. Opt. **36**, 3948–3958 (1997). [CrossRef] [PubMed]