## Abstract

A non-iterative method based on principal component analysis (PCA) is presented to directly extract the phase from multiple-beam Fizeau interferograms with random phase shifts. The PCA method is the approach that decomposes the multiple-beam Fizeau interferograms into many uncorrelated quadrature signals and then applies principal component analysis algorithm to extract the measured phase without any prior guess about the phase shifts. Some factors that affect the performance of the proposed method are analyzed and discussed. Numerical simulations and experiments demonstrate that the proposed method extracts phase fast and exhibits high precision. The method can be applied in high precision interferometry.

© 2011 OSA

## 1. Introduction

Phase-shifting interferometry (PSI) [1,2] has been widely used in surface testing of optical elements. One can usually evaluate the phase map of interferograms by means of two-beam phase-stepping algorithms [1–4], which allows one to achieve great accuracy if the interferogram intensity is a sinusoidal dependence on the phase map. However, some interferometers of great practical interest do not produce sinusoidal profiles [5,6]. For example, in a Fizeau interferometer, because of multiple reflections between reference surface and test surface, fringes with an Airy profile are obtained [7]. Thus significant phase error will occur if two-beam phase-stepping algorithms are employed [8], especially for the measurement of test surface with high reflectivity coefficient.

To deal with this problem, methods of reducing multiple reflections are proposed. Clapham used a lossy reflecting coating on the reference surface [9], and Hariharan introduced an absorbing filter in the test path between the reference surface and the test surface [10]. However, Clapham’s method will modify the figure of the reference surface and Hariharan’s method needs two sets of measurement to reduce the error caused by the absorbing filter. In addition, although the two methods can reduce multiple reflections, it cannot eliminate them completely.

Therefore, some special algorithms are presented to reduce the effect of multiple reflections.

If a four-frame phase calculation algorithm is used, the phase error caused by multiple reflections of a mirror is eliminated to the first-order approximation [5]. Bonsch and Bohme [7] presented an algorithm that can completely eliminate the phase error (mainly quadruple frequency) introduced by the multiple reflections of a mirror. Hibino et al. [11] and Surrel [12] offered some methods which show that the *j*th-order harmonic can be minimized with phase steps of$2\pi /(j+2)$ between acquired data frames. These special algorithms require accurate phase shifts; otherwise they will fail to reduce the effect of multiple reflections.

Unfortunately, it is hard to obtain accurate phase shifts in practical measurement because of miscalibration error of PZT and environmental vibration [13]. On the one hand, Picart et al. [14], Patil et al. [15] and Langoju et al. [16,17] proposed some methods to estimate phase information from interferograms with harmonics and unknown phase steps. These methods mainly deal with the phase step, linear and nonlinear miscalibration errors in the PZT response, but fail to consider the problem of random phase shifts caused by environmental vibration. On the other hand, some iterative algorithms [18–20] and non-iterative algorithms [21–25] are also developed to estimate the phase shifts from interferograms with random phase shifts, but these algorithms fail to consider the problem of high-order harmonics. Therefore, none of the mentioned algorithms [5,7,11–25] can effectively deal with the problem of high-order harmonics and random phase shifts.

Recently, the author has presented an iterative algorithm to analyze the multiple-beam Fizeau interferograms with random phase shifts [26]. It can extract phase distribution from seven multiple-beam interferograms with random phase shifts by iteration process. It requires the multiple-beam interferograms have large spatial carrier frequency, so that it can estimate the phase shifts by FFT method [21]. Otherwise, the iteration process is time-consuming and leads to substantial computation loads.

Here we propose a non-iterative approach to directly extract the phase distribution from multiple-beam Fizeau interferograms with random phase shifts. The method is based on principal component analysis (PCA), and it is similar as the method proposed by Vargas [27]. However, Vargas uses PCA to extract phase from sinusoidal interferograms, while we use PCA to extract phase from multiple-beam Fizeau interferograms. We will first discuss the principle of our method, and then give its verification by computer simulations and experiments.

## 2. Multiple-beam Fizeau interferograms

The Fizeau intensity profile in reflection for two nonabsorbing nearly parallel surfaces is given by the Airy formula that can be written in compact form as [7]

where *A* is the local mean intensity, and *B* and *C* are constants that depend on the reflectivity coefficients of reference surface *r*
_{1} and test surface *r*
_{2}

In the Eq. (1),*φ*, ${\theta}_{n}$and${\eta}_{n}$denote the measured phase, the phase shift and the noise, respectively. Although Eq. (1) is clearly a nonsinusoidal profile, the term in brackets is a periodic and even function of$\phi +{\theta}_{n}$. Then it can be conveniently developed in Fourier series form, and Eq. (1) can be rewritten as [8,26]

It can be seen that the amplitude of the *j*th harmonic decreases with increasing of *j*. Moreover, the rate of amplitude between consecutive harmonics is equal to${r}_{1}{r}_{2}$. According to Eq. (3), the problem of multiple-beam Fizeau interferograms with random phase shifts is, in fact, the problem of interferograms with high-order harmonics and random phase shifts. Equation (3) can be rewritten as

In Eq. (5), ${b}_{j}=\mathrm{cos}(j{\theta}_{n})$ and ${c}_{j}=-\mathrm{sin}(j{\theta}_{n})$. According to probability theory, If we have more than one fringe in the field of view, then the covariances of $\mathrm{cos}(j\phi )$ and $\mathrm{sin}(k\phi )$, can be approximately expressed as

where ${M}_{x}$and${M}_{y}$ mean the size of interferogram along x and y axes, respectively. Equation (6) shows that the covariance of $\mathrm{cos}(j\phi )$and $\mathrm{sin}(k\phi )$approximately equal to zero. Similarly,

Thus we can easily prove that the covariance of $\mathrm{cos}(j\phi )$ and$\mathrm{cos}(k\phi )$, and the covariance of $\mathrm{sin}(j\phi )$ and$\mathrm{sin}(k\phi )$approximately equal to zero when$j\ne k$. The more fringes in the field of interferograms have, the more accurate Eqs. (6) and (7) will be. In probability theory, two variables are said to be uncorrelated if their covariance is zero. So Eq. (5) shows that any multiple-beam Fizeau interferograms without background term $\frac{A{a}_{0}}{2}$ and noise term ${\eta}_{n}$ can be decomposed into many uncorrelated quadrature signals $\mathrm{cos}(j\phi )$ and $\mathrm{sin}(k\phi )$.

## 3. Principal component analysis

Principal component analysis (PCA) [28] is a mathematical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of uncorrelated variables called principal components. The principal components are linear combinations of the original variables and are the single best subspace of a given dimension in the least-square sense. In practice, Our PCA algorithm includes four steps.

#### 3.1 Calculating the covariance matrix

Suppose that we have *N* images of multiple-beam Fizeau interferograms and each image is reshaped into one column with size of *M*, (here$M={M}_{x}\times {M}_{y}$). Then this *N* images can be expressed in a matrix form as

where ${I}_{n}$means the *n*th image with size of *M*, and ${[]}^{T}$denotes the transposing operation. Equation (8) shows that the size of matrix of *I* is $N\times M$. Suppose that the mean of *I* is ${\mu}_{I}={\displaystyle \sum _{n=1}^{n=N}\frac{{I}_{n}}{N}}$, and then the covariance matrix of the data set can be expressed as

In Eq. (9), $I-{\mu}_{I}$corresponds to a background suppression operation. The accuracy of background suppression operation increases with the number of fringes used (*N*). The size of ${C}_{I}$is $N\times N$. The components of${C}_{I}$, denoted by${C}_{pq}$, represent the covariances between the variable components${I}_{p}-{\mu}_{I}$ and ${I}_{q}-{\mu}_{I}$ and represent the variance of the component ${I}_{p}-{\mu}_{I}$ when *p* = *q*. The covariance matrix of ${C}_{I}$is, by definition, always symmetric.

#### 3.2 Diagonalization process

From a symmetric matrix of${C}_{I}$, we can calculate an orthogonal basis by finding its eigenvalues and eigenvectors. The eigenvectors ${\varphi}_{i}$and the corresponding eigenvalues ${\lambda}_{i}$are the solutions of the equation

We can construct a matrix *Φ* consisting of${\varphi}_{i}$,

Since the matrix of ${C}_{I}$ is a square matrix, the eigenvectors ${\varphi}_{i}$of ${C}_{I}$ is an orthogonal basis and the matrix Φ is orthogonal (unitary). So we have

where *Λ* is a diagonal matrix consisting of ${\lambda}_{i}$ and ${\Phi}^{-1}$is the conjugate transpose of *Φ*. In fact, Eq. (12) is the diagonalization process of the symmetric matrix${C}_{I}$. In our paper the diagonalization process is performed by the singular value decomposition (SVD) algorithm. By ordering the eigenvectors in the order of descending eigenvalues (largest first), one can create an ordered orthogonal basis with the first eigenvector having the direction of largest variance of the data. In this way, we can find directions in which the data set has the most significant amounts of energy.

#### 3.3 Principal Component Transform

We can now define the orthogonal principal component transform of *I*. The forward transform [27,28] is

We are free to create new variables Y based on linear combinations of the existing ones. Components of ** Y** can be seen as the coordinates in the orthogonal base. We can reconstruct the original data vector $(I-{\mu}_{I})$from

**by**

*Y*Equation (14) shows that the original vector $(I-{\mu}_{I})$is projected on the coordinate axes defined by the orthogonal basis *Y*. This minimizes the mean-square error between the data and this representation with given number of eigenvectors. Instead of using all the eigenvectors of the covariance matrix, we may represent the data in terms of only a few basis vectors of the orthogonal basis. For the analysis of multiple beam Fizeau interferograms expressed by Eq. (1), it has only five unknows (A, B, C, *φ* and${\theta}_{n}$), thus more than four eigenvectors are required (i.e.$N\ge 5$).

#### 3.4 Calculating the measured phase

The first two principal components with the biggest and second biggest eigenvalues (${y}_{1}$and${y}_{2}$) correspond to ${a}_{1}\mathrm{cos}(\phi )$and${a}_{1}\mathrm{sin}(\phi )$, respectively. The third and fourth biggest eigenvalues (${y}_{3}$and${y}_{4}$) correspond to ${a}_{2}\mathrm{cos}(2\phi )$and${a}_{2}\mathrm{sin}(2\phi )$, respectively. Therefore, the measured phase is obtained as

In our case, we only use ${y}_{1}$and ${y}_{2}$ to calculate the measured phase because bigger eigenvalues leads to higher accuracy. Our method cannot determine the correct global phase sign as we arbitrarily assign $\mathrm{cos}(\phi )$ and $\mathrm{sin}(\phi )$ to the first and second principal component, respectively.

## 4. Numerical simulation and discussion

Numerical simulations are carried out to test the performance of the proposed algorithm. We assume that $A=\mathrm{exp}[-0.02({x}^{2}+{y}^{2})]$and$\phi =5({x}^{2}+{y}^{2})+\pi x+\pi y$, where$-1\le x\le 1$, $-1\le y\le 1$. In addition, ${r}_{1}$ = 0.2, ${r}_{2}$ = 0.7, $\delta (n)$is a random number ranging from 0 to 2*π*, and ${\eta}_{n}$is white Gaussian noise with signal-to-noise ratio (SNR) of 30 dB. According to Eq. (1), we generate 7 multiple-beam Fizeau interferograms with random phase shifts, and one of them is shown in Fig. 1 (a)
with the contrast about 0.29. Then we extract the measured phase from the multiple-beam Fizeau interferograms by the proposed PCA algorithm. The normalized eigenvalues of matrix ${C}_{I}$ are 1.0000, 0.8641, 0.0360, 0.0228, 0.0148 and 0.0147 and 0. We use the first two principal components to calculate the measured phase, and the wrapped and unwrapped phases are shown in Figs. 1 (b) and 1(c), respectively. By comparing the unwrapped phase with the given phase *φ*, the phase error of our algorithm is obtained and shown in Fig. 1 (d). The peak-to-valley (PV) and root-mean-square (RMS) values of the phase error are 0.670 and 0.103 rad, respectively. Figure 1 shows that the proposed algorithm can accurately extract the phase from seven multiple-beam Fizeau interferograms with random phase shifts.

In the above simulation, we find that some factors (such as reflectivity coefficient, random noise and frame number used) will affect the performance of the proposed algorithm. These factors are analyzed and discussed as follows.

First, we assume that the number of multiple-beam Fizeau interferograms used (*N*) in our simulation is 7 and each frame has signal-to-noise ratio (*SNR*) of 30 dB. We let *r*
_{1} = 0.2 and *r*
_{2} vary from 0.1 to 1. Then we calculate the RMS value of phase error for different *r*
_{2} and the result is shown in Fig. 2 (a)
. For$0.1\le {r}_{2}\le 0.3$, the contrast of fringes is higher than 0.78 and the rate of amplitude between consecutive harmonics is less than 0.06. The relative amplitude of high-order harmonics is small and its effect is not significant, thus the RMS value of phase error varies slowly with *r*
_{2} when$0.1\le {r}_{2}\le 0.3$. However, as the increasing of *r*
_{2}, the contrast of the multiple-beam fringe decreases, the relative amplitude of high-order harmonics increases, and the influence of noise becomes serious. In addition, the covariance of $\mathrm{cos}(\phi )$ and $\mathrm{cos}(2\phi )$or $\mathrm{cos}(3\phi )$is a small value but not absolutely equals to zero, which will affect the accuracy of phase extraction. Therefore, the RMS value of phase error increases evidently with the increasing of *r*
_{2} when${r}_{2}\ge 0.3$, as shown in Fig. 2 (a). Figure 2 (a) also shows that the RMS value of phase error is less than 0.15 rad even if *r*
_{2} approximates to 1.

Second, we assume that *r*
_{1} = 0.2, *r*
_{2} = 0.7, and *N* = 7. Then we perform simulation at different *SNRs* to compute the RMS value of phase error, and the result is shown in Fig. 2(b). Figure 2(b) shows that the RMS value of phase error decreases with the increasing of *SNR*. When SNR becomes larger than 40, the effect of noise on phase error becomes insignificant, so the RMS value of phase error keeps nearly constant of 0.025rad.

Third, we assume that *r*
_{1} = 0.2, *r*
_{2} = 0.7, and *SNR* = 30dB. Then we calculate the RMS value of phase error when we use different numbers of multiple-beam fringes, and the result is shown in Fig. 2(c). Figure 2(c) shows that the RMS value of phase error decreases as the used number of multiple-beam fringes increases. From Figs. 2 (a-c) we can conclude that the RMS value of phase error is less than 0.1rad when${r}_{2}\le 0.7$, $N\ge 7$and$SNR\ge 30dB$.

## 5. Experiments

To illustrate the application of the proposed PCA algorithm, we apply it to the practical multiple-beam Fizeau interferograms. In the experiment, the reflectivity coefficients of the reference and test surfaces in Fizeau interferometer are 0.2 and 0.7, respectively. We capture 7 frames of fringes in dynamic environments and one typical fringe is shown in Fig. 3 (a) with size of 200 × 200. We process these fringes by the proposed PCA algorithm and the least-squares based iterative algorithm [26] and the processing times of them are 0.33s and 35.82s (10 iterative cycles) respectively. Figures 3 (b) and 3(c) show the retrieved phase by the proposed PCA algorithm and the least-squares based iterative algorithm, respectively. In addition, we also measure the test surface by Zygo’s interferometer with a vibration-isolating platform, calibrated PZT, and a proper optical attenuator in Fizeau cavity [10]. The obtained phase is used as the reference phase of the test surface and shown in Fig. 3 (d). By comparing Fig. 3 (d) with Figs. 3 (b) and 3(c), the phase errors of the proposed PCA algorithm and the least-squares based iterative algorithm are obtained and shown in Figs. 3 (e) and 3(f), respectively. In Table 1 we show the peak-to-valley (PV) and RMS values of Figs. 3 (b-f).

As can be seen from Fig. 3 and Table 1, the retrieved phase by the proposed PCA algorithm exhibits high precision with very weak waviness, while the retrieved phase by the least-squares based iterative algorithm exhibits visible waviness because of the effect of high-order harmonics and random phase shifts. The RMS value of the phase error of the proposed PCA algorithm is 0.0511rad, which is less than that of the least-squares based iterative algorithm. If more fringes are processed, the proposed PCA algorithm can further improve the accuracy with a short processing time. For the least-squares based iterative algorithm, it can reduce the phase error by more iterative cycles and more fringes, but the processing time increases rapidly with the numbers of iterative cycles and fringes. Comparing with the least-squares based iterative algorithm, the proposed PCA algorithm does not need any prior guess about phase shifts, works faster and more efficiently, and exhibits higher precision. Therefore, the proposed PCA algorithm is more suitable for practical multiple-beam Fizeau interferograms than the least-squares based iterative algorithm.

## 6. Conclusion

To conclude, we have proposed a non-iterative method for phase extraction from multiple-beam Fizeau interferograms with random phase shifts. The method is based on principal component analysis and does not need any prior guess about the phase shifts. It decomposes the multiple-beam Fizeau interferograms into many uncorrelated quadrature signals and then extracts the phase from the quadrature signals. The numerical simulations and experiments demonstrate the validity of our method. Numerical simulations show that the RMS value of phase error of the proposed PCA algorithm is less than 0.1rad when${r}_{2}\le 0.7$, $SNR\ge 30dB$ and$N\ge 7$. Experiments show that the proposed PCA algorithm is very fast and easy to implement and exhibits high precision. Thus it is more suitable for practical multiple-beam Fizeau interferograms than the least-squares based iterative algorithm. The proposed method has applications in high precision interferometry, especially for measurement of surface whose reflectivity coefficient is greater than 0.3. Furthermore, it can also be used to analyze other form of nonsinusoidal fringes [29] with random phase shifts.

## Acknowledgments

This research is supported by the Natural Science Foundation of Zhejiang Province, China (Y1110125). Furthermore, we would like to express our gratitude to Javier Vargas for invaluable help and useful discussion.

## References and links

**1. **K. Creath, “Temporal phase measurement method, ” in *Interferogram Analysis,* D. W. Robinson and G. T. Reid, eds. (Institute of Physics, Bristol, UK, 1993), pp. 94–140.

**2. **D. Malacara, M. Servín, and Z. Malacara, *Interferogram Analysis for Optical Testing* (Marcel Dekker, New York, 1998).

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

**4. **Y. Surrel, “Fringe analysis,” in *Photomechanics*, P. K. Rastogi, ed., Vol. 77 of Topics in Applied Physics (Springer, 2000), pp. 55–102.

**5. **P. Hariharan, “Digital phase-stepping interferometry: effects of multiply reflected beams,” Appl. Opt. **26**(13), 2506–2507 (1987). [CrossRef] [PubMed]

**6. **C. Ai and J. C. Wyant, “Effect of retroreflection on a Fizeau phase-shifting interferometer,” Appl. Opt. **32**(19), 3470–3478 (1993). [CrossRef] [PubMed]

**7. **G. Bonsch and H. Bohme, “Phase-determination of Fizeau interferences by phase-shifting interferometry,” Optik (Stuttg.) **82**, 161–164 (1989).

**8. **B. V. Dorrío, J. Blanco-García, C. López, A. F. Doval, R. Soto, J. L. Fernández, and M. Pérez-Amor, “Phase error calculation in a Fizeau interferometer by Fourier expansion of the intensity profile,” Appl. Opt. **35**(1), 61–64 (1996). [CrossRef] [PubMed]

**9. **P. B. Clapham and G. D. Dew, “Surface-coated reference flats for testing fully aluminized surfaces by means of the Fizeau interferometer,” J. Sci. Instrum. **44**(11), 899–902 (1967). [CrossRef]

**10. **P. Hariharan, “Interferometric measurements of small-scale irregularities: highly reflecting surfaces,” Opt. Eng. **37**(10), 2751–2753 (1998). [CrossRef]

**11. **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**(4), 761–768 (1995). [CrossRef]

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

**13. **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**(21), 3421–3432 (1983). [CrossRef] [PubMed]

**14. **P. Picart, R. Mercier, and M. Lamare, “Influence of multiple-beam interferences in a phase-shifting Fizeau interferometer and error-reduced algorithms,” Pure Appl. Opt. **5**(2), 167–194 (1996). [CrossRef]

**15. **A. Patil and P. Rastogi, “Estimation of multiple phases in holographic moiré in presence of harmonics and noise using minimum-norm algorithm,” Opt. Express **13**(11), 4070–4084 (2005). [CrossRef] [PubMed]

**16. **R. Langoju, A. Patil, and P. Rastogi, “Resolution-enhanced Fourier transform method for the estimation of multiple phases in interferometry,” Opt. Lett. **30**(24), 3326–3328 (2005). [CrossRef] [PubMed]

**17. **R. Langoju, A. Patil, and P. Rastogi, “Phase-shifting interferometry in the presence of nonlinear phase steps, harmonics, and noise,” Opt. Lett. **31**(8), 1058–1060 (2006). [CrossRef] [PubMed]

**18. **Z. Wang and B. Han, “Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms,” Opt. Lett. **29**(14), 1671–1673 (2004). [CrossRef] [PubMed]

**19. **H. Guo and M. Chen, “Least-squares algorithm for phase-stepping interferometry with an unknown relative step,” Appl. Opt. **44**(23), 4854–4859 (2005). [CrossRef] [PubMed]

**20. **X. F. Xu, L. Z. Cai, X. F. Meng, G. Y. Dong, and X. X. Shen, “Fast blind extraction of arbitrary unknown phase shifts by an iterative tangent approach in generalized phase-shifting interferometry,” Opt. Lett. **31**(13), 1966–1968 (2006). [CrossRef] [PubMed]

**21. **K. A. Goldberg and J. Bokor, “Fourier-transform method of phase-shift determination,” Appl. Opt. **40**(17), 2886–2894 (2001). [CrossRef] [PubMed]

**22. **X. F. Xu, L. Z. Cai, Y. R. Wang, X. F. Meng, W. J. Sun, H. Zhang, X. C. Cheng, G. Y. Dong, and X. X. Shen, “Simple direct extraction of unknown phase shift and wavefront reconstruction in generalized phase-shifting interferometry: algorithm and experiments,” Opt. Lett. **33**(8), 776–778 (2008). [CrossRef] [PubMed]

**23. **P. Gao, B. Yao, N. Lindlein, J. Schwider, K. Mantel, I. Harder, and E. Geist, “Phase-shift extraction for generalized phase-shifting interferometry,” Opt. Lett. **34**(22), 3553–3555 (2009). [CrossRef] [PubMed]

**24. **T. E. Zander, V. Madyastha, A. Patil, P. Rastogi, and L. M. Reindl, “Phase-step estimation in interferometry via an unscented Kalman filter,” Opt. Lett. **34**(9), 1396–1398 (2009). [CrossRef] [PubMed]

**25. **J. Xu, Y. Li, H. Wang, L. Chai, and Q. Xu, “Phase-shift extraction for phase-shifting interferometry by histogram of phase difference,” Opt. Express **18**(23), 24368–24378 (2010). [CrossRef] [PubMed]

**26. **J. Xu, Q. Xu, L. Chai, and H. Peng, “Algorithm for multiple-beam Fizeau interferograms with arbitrary phase shifts,” Opt. Express **16**(23), 18922–18932 (2008). [CrossRef] [PubMed]

**27. **J. Vargas, J. A. Quiroga, and T. Belenguer, “Phase-shifting interferometry based on principal component analysis,” Opt. Lett. **36**(8), 1326–1328 (2011). [CrossRef] [PubMed]

**28. **R. C. Gonzalez and R. E. Woods, *Digital Image Processing*, 3rd ed. (Prentice-Hall, 2007).

**29. **B. Pan, Q. Kemao, L. Huang, and A. Asundi, “Phase error analysis and compensation for nonsinusoidal waveforms in phase-shifting digital fringe projection profilometry,” Opt. Lett. **34**(4), 416–418 (2009). [CrossRef] [PubMed]