## Abstract

A method combining the principal component analysis (PCA) and the least squares method (LSM) is proposed to extract the phase from interferograms with random phase shifts. The method estimates the initial phase by PCA, and then determines the correct global phase sign and reduces the residual phase error by LSM. Some factors that may influence the performance of the proposed method are analyzed and discussed, such as the number of frames used, the number of fringes in interferogram and the amplitude of random phase shifts. Numerical simulations and optical experiments are implemented to verify the effectiveness of this method. The proposed method is suitable for randomly phase-shifted interferograms.

©2011 Optical Society of America

## 1. Introduction

Phase-shifting interferometry (PSI) has been widely used in surface testing of optical elements, especially the large-aperture mirrors in astronomical telescopes or satellite cameras. Owing to the large size of the measured surface, the length of optical path is set to be several meters or even longer to complete the testing. However, the mechanical vibration and the air turbulence will change the preset phase shifts and cause inevitable errors to the measurement [1]. The traditional PSI algorithms [2–9] require a special constant phase shift in each step (usually 2π/*K*, where *K* is a positive integer). Thus it imposes a strict requirement on the piezo transducers (PZT)’s accuracy and the environmental stability. However, it is often difficult to meet the requirement exactly in practice, especially for the measurement of large-aperture optics.

To deal with the problem of random phase shifts, on the one hand, many numerical iterative methods [10–16] have been proposed to determine the phase-shift amounts and the phase distribution simultaneously. These iterative methods need three or more frames of interferograms and require an iteration process that will lead to substantial computation loads. On the other hand, some non-iterative methods are also presented. Dobroiu *et al* proposed the statistical self-calibrating algorithms to extract the phase distribution on the assumption of a constant fringe contrast [17], a quasi-uniformly distributed phase in the range of [0, 2π] [18] or an uniform illumination [19]. Cai *et al* [20], Xu *et al* [21,22], Meng *et al* [23] and Gao *et al* [24] introduced some algorithms to directly extract the unknown phase shifts on the assumption that the measured surface has an equal distribution in [0,2π] over the whole interferogram. However, these assumptions are too strong to be met exactly in practical experiments. For example, the measured phase does not distribute equally in [0, 2π] when the peak-to-valley value of the measured phase is less than 2π.

To relax the requirement of the measured surface, Goldberg and Bokor [25] used Fourier transform to determine the phase shift between two frames by introducing a large carrier frequency into the interferograms. The phase-shift estimation error cannot be ignored if the carrier frequency is not large enough. Chen *et al* [26,27] proposed a max-min algorithm based on two-point correlation and Xu *et al* [28] presented a non-iterative approach to extract the unknown phase shift according to the histogram of phase difference without the assumption of equal distribution of measured phase in [0,2π]. However, both Chen’s and Xu’s methods need to determine the background and the modulation amplitude of the fringe from many frames of interferograms by finding the maximum and the minimum intensity in each pixel [26–28].

Recently, Vargas *et al* [29,30] proposed a new method based on principal component analysis (PCA), which can fast extract the phase distribution from randomly phase-shifted interferograms with a weaker assumption of the measured surface. The PCA is an efficient technique for phase extraction by converting a set of observations of possibly correlated variables into a set of values of uncorrelated variables [29].The PCA method does not need to determine the background and the modulation amplitude of the fringe before the phase measurement. We have extended the PCA technique to multiple beam interferometry [31] and found that the PCA method is very fast and easy to implement in practical experiments. However, we also found that the PCA method cannot determine the correct global sign of the measured phase and the obtained phase still has an unignorable residual phase error [29–31]. These two problems are very important for the accurate phase measurement of practical optics. Therefore, in this paper, we propose a method combining the principal component analysis (PCA) and the least squares method (LSM) to solve the two problems existing in PCA method. Firstly, the estimated phase by PCA is worked as the initial phase, which may have a wrong global phase sign and an unignorable residual phase error. Secondly, the correct global phase sign is determined and the residual phase error is reduced by least squares method (LSM), which may include one or a few iterative cycles. In this paper, we refer the proposed method as the “PCA+LSM” method for brevity. We will discuss the principle of the “PCA+LSM” method, and then give its verification by computer simulations and experiments.

## 2. Principal component analysis

Suppose *N* frames of randomly phase-shifted interferograms are collected and each image is reshaped into one column with size of *M*, where *M* means the number of pixels in each image. The intensity of the *m*th pixel in the *n*th interferogram is expressed as

*${B}_{m}$*and ${\phi}_{m}$represent the background intensity, the modulation amplitude and the measured phase of the fringe in the

*m*th pixel. ${\theta}_{n}$is the unknown phase shift for the

*n*th frame. Then these

*N*images can be expressed in a matrix form aswhere${[]}^{T}$denotes the transposing operation. Supposing that the mean of $I$is ${\mu}_{I}={\displaystyle \sum _{n=1}^{n=N}{I}_{n}/N}\approx A$, then $I-{\mu}_{I}$means the intensity of

*N*frames without background intensity and can be approximately expressed as

In probability theory, two variables are said to be uncorrelated if their covariance is zero. So Eqs. (3)–(6) show that $[I-{\mu}_{I}]$can be decomposed into two approximately uncorrelated quadrature signals $u$and$v$. So now the question is how to obtain the uncorrelated quadrature signals $u$and $v$ from$[I-{\mu}_{I}]$. We solve it by principal component analysis (PCA), which 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 [32].

According to the probability theory, the covariance matrix of *N* images can be expressed as

*D*is a diagonal matrix consisting of ${\lambda}_{i}$. Once the covariance matrix is diagonalized and the orthogonal matrix $\Phi $is obtained, the principal components can be calculated by principal component transform [29,34]where

*Y*consists of the principal components of

*$I-{\mu}_{I}$,$Y={[{y}_{1},{y}_{2},\cdot \cdot \cdot ,{y}_{n},\cdot \cdot \cdot ,{y}_{N}]}^{T}$*. The

*i*th component of

*Y*is${y}_{i}={\varphi}_{i}{}^{T}(I-{\mu}_{I})$, which means the projection of $I-{\mu}_{I}$on ${\varphi}_{i}$. The first two principal components with the biggest and second biggest eigenvalues (${y}_{1}$and${y}_{2}$) correspond to $u$and $v$, or $v$and $u$, respectively [29]. Therefore, the measured phase is estimated as

Here, we cannot determine the correct global phase sign as we arbitrarily assign the cosine and sine signals ($u$and$v$) to the first and second principal components. In addition, the accuracy of PCA depends on some approximations [30], such as the background intensity estimation and the approximations of Eqs. (5)–(6). However, in a real experiment, neither of these approximations will be fully accomplished. In this case, the estimated phase will have an unignorable residual phase error. Although the estimated phase by PCA method may have a wrong global phase sign and an unignorable phase error, the PCA method works well and obtains an approximating phase distribution even if the phase shifts are complete random. Thus, the estimated phase by PCA method works well as the initial phase distribution of the least-squares based iterative method.

## 3. Least squares method

Once$\phi $ is estimated, the phase shift ${\theta}_{n}$can be easily determined by using of the least squares method [13]. We assume that the background intensity and the modulation amplitude do not have pixel-to-pixel variation, i.e. ${A}_{1}={A}_{m}={A}_{M}$ and ${B}_{1}={B}_{m}={B}_{M}$. Then ${a}_{n}={A}_{m}$, ${b}_{n}={B}_{m}\mathrm{co}s{\theta}_{n}$ and ${c}_{n}=-{B}_{m}\mathrm{sin}{\theta}_{n}$ can be determined by

*N*and it is wrapped in the range from $-\pi $to$\pi $. We unwrap the estimated $\stackrel{~}{{\theta}_{n}}$ by one-dimension phase unwrapping, $\stackrel{\wedge}{\theta}=unwrap(\stackrel{~}{\theta})$. In practical experiments, we can control the first real phase step to ensure$0<{\theta}_{2}-{\theta}_{1}<\pi $. Therefore, if $\stackrel{\wedge}{{\theta}_{2}}-\stackrel{\wedge}{{\theta}_{1}}$ is less than zero, then ${\theta}_{n}=-{\stackrel{~}{\theta}}_{n}$; otherwise, ${\theta}_{n}={\stackrel{~}{\theta}}_{n}$. In this way, we can determine the correct direction of phase shift.

Once the phase shifts are estimated with correct direction, the phase distribution can be calculated with correct global sign by using least squares method [10–13]. We assume that ${A}_{1}={A}_{n}={A}_{N}$and ${B}_{1}={B}_{n}={B}_{N}$.Then ${a}_{m}={A}_{n}$, ${b}_{m}={B}_{n}\mathrm{co}s{\phi}_{m}$ and ${c}_{m}=-{B}_{n}\mathrm{sin}{\phi}_{m}$ can be can be determined by

## 4. Numerical simulation and discussion

Numerical simulations are carried out to test the performance of the proposed method. *N* frames are generated according to Eq. (1) by setting the parameters as follows:$A=130$,*$B=120$*, and $\phi =D\xb7\pi [({x}^{2}+{y}^{2})]$, where *$-1\le x\le 1$*, $-1\le y\le 1$ and *D* is a parameter that denotes the number of fringes in the filed of interferogram. ${\theta}_{1}=0$and ${\theta}_{n}=(n-1)\pi /2+E\xb7r(n)$ for *n*>1, where $r(n)$ is random rational number ranging from −1 to 1 and E denotes the amplitude of random phase shifts. In addition, a Gaussian white noise with mean zero and standard deviation 1 is added to each fringe. Then the intensities of fringes are assumed to be gathered by a CCD with 8-bit quantization.

We generate three interferograms with *D* = 1.5 and *E* = 1. One typical interferogram and the given measured phase are shown in Figs. 1(a)
and 1(b). The three interferograms are analyzed by the PCA method and the “PCA+LSM” method, and the extracted phases of them are shown in Figs. 1(c) and 1(e) respectively. By comparing Figs. 1(c) and 1(e) with the given measured phase, the residual phase errors of the PCA method and the “PCA+LSM” method are obtained and shown in Figs. 1(d) and 1(f), respectively. The RMS values of Figs. 1(d) and 1(f) are 0.0382 rad and 0.0157 rad. It can be concluded from Figs. 1(b)–(d) that, the phase extracted by PCA method has an incorrect global sign while the “PCA+LSM” method can correct the global sign. In addition, Figs. 1(d) and 1(f) indicate that the “PCA+LSM” method can reduce the residual phase error and thus has a smaller phase error than the PCA method.

In the simulations, we find that some factors may affect the performance of the proposed algorithm, such as the number of frames used (*N*), the number of fringes in interferogram (*D*), and the amplitude of random phase shifts (*E*). These factors are analyzed and discussed as follows. Here we define the phase-shift extraction error as the average of the absolute values of the difference between the extracted phase shifts and the given phase shifts over *N* frames.

Firstly, we assume that *D* = 1.5 and *E* = 1, and let *N* vary from 3 to 20. Then we use the “PCA+LSM” method to analyze these interferograms. The phase-shift extraction error and the phase error for different *N* are obtained and shown in Figs. 2(a)
and 2(b), respectively. Meanwhile, the PCA method is also used for comparison, and the phase error is given in Fig. 2(b). Figure 2(a) shows that the phase-shift extraction error oscillates between 0.026 rad and 0.031 rad as *N* increases from 3 to 20. Figure 2(b) shows that the phase errors of the PCA method and “PCA+LSM” method decreases slowly with the increasing of *N,* but the latter is less than one third of the former for a given *N*. It can be also concluded from Figs. 2(a) and 2(b) that, for *D* = 1.5, *E* = 1 and *N*≥3, the phase-shift extraction error and the RMS value of phase error of the proposed method are less than 0.031 rad and 0.012 rad, respectively.

Secondly, we assume that *N* = 3 and *E* = 1 and let *D* vary from 1 to 20. The PCA method and the “PCA+LSM” method are applied and the results are shown in Fig. 3
. Since the approximations of Eqs. (5)–(6) are more close to the truth for a larger *D*, the phase-shift extraction error and the phase error decrease with the increasing of *D*, just as shown in Figs. 3 (a) and 3(b). Figure 3(b) indicates that the “PCA+LSM” method is more accurate than the PCA method, but the results of them will converge slowly as *D* increases. Figures 3(a) and 3(b) also show that the phase-shift extraction error and the RMS value of phase error of the “PCA+LSM” method are less than 0.035 rad and 0.03 rad for *D*>1.

Thirdly, we assume that *N* = 3 and *D* = 2.5, and let *E* vary from 0.1 to 3. The PCA method and the “PCA+LSM” method are used and the phase errors of them are shown in Fig. 4(a)
. Figure 4(a) shows that the phase errors are kept under a small value for *E*<1. For *N* = 3 and *E*<1, the relative phase shifts meet the condition that **$\text{0}\text{.57}<{\theta}_{2}-{\theta}_{1}<2.57$**and$\text{2}\text{.14}<{\theta}_{3}-{\theta}_{1}<4.14$, so it is easy to determine the global signs of phase shift and phase distribution. However, for *N* = 3 and *E*>1.5, the relative phase shifts ${\theta}_{2}-{\theta}_{1}$ and ${\theta}_{3}-{\theta}_{1}$ may close to 2*kπ,* where *k* = 0, 1, 2….In this case, the three interferograms become more interdependent and finally it results in a big phase error, just as shown in Fig. 4 (a). This problem can be solved by increasing the number of frames used. Figure 4(b) shows the phase errors of the two methods at different *E* for *N* = 5 and *D* = 2.5. For *N*≥5, it is easy to ensure that there are two or more frames whose relative phase shifts are not close to 2*kπ.* Thus there is no big phase error in Fig. 4(b) even if *E*>1.5. It can be concluded from Figs. 4(a) and 4(b) that we should collect more frames to improve the performance of the proposed method if the phase shifts are completely random.

Finally, to demonstrate the performance of our method, we make a comparison between our method and Gao’s method [24] published in 2009. Both of the two methods can extract the unknown phase shifts from three interferograms, so we assume that *N* = 3 and *E* = 1 and let *D* vary from 0.25 to 5. The phase-shift extraction errors of our method and Gao’s method are obtained and shown in Fig. 5
. It shows that the phase-shift extraction error of our method decreases with the increasing of *D* and the error is less than 0.12 for$D\ge 0.5$, while the phase-shift extraction error of Gao’s method is a damping oscillation of *D* and the error is less than 0.12 for$D\ge 1.\text{75}$. Figure 5 also indicates that phase-shift extraction error of our method is less than that of Gao’s method in most general case.

## 5. Experiment

Optical experiments have also been carried out to investigate the performance of our method. Three randomly phase-shifted interferograms, with size of 512 × 512, are collected. One of them is shown in Fig. 6(a) . Firstly, the advanced iterative algorithm (AIA) [13], with the initial phase shifts of 0, π/2 and π, is used to extract the phase distribution and the phase shifts by 20 iterative cycles. The extracted phase shifts are 0, 0.9460, and 3.2272 rad and the extracted phase is shown in Fig. 6(b), which works as the reference phase. Next, we process the three interferograms with the PCA method and the “PCA+LSM” method, and the results of them are shown in Figs. 6(c) and 6(d) respectively. By comparing with the reference phase, the phase errors of the PCA method and the “PCA+LSM” method are obtained and shown in Figs. 6(e) and 6(f). The RMS values of Figs. 6(e) and 6(f) are 0.0137 rad and 0.0058 rad. In addition, the fringe patterns are processed using MATLAB, and the processing times of the AIA method, the PCA method and the “PCA+LSM” method are 23.063s, 0.329s and 1.498s, respectively. Finally, some conclusions are drawn from the experimental results mentioned above. (1) The AIA method can accurately extract the phase distribution and the phase shifts by iterative cycles, but it is an extremely time-consuming process for randomly phase-shifted interferograms. (2) The PCA method is fast, but it cannot determine the correct global phase sign and has a visible residual phase error. (3) The “PCA+LSM” method can determine the correct global phase sign and has a smaller phase error than the PCA method. In addition, the processing time of the “PCA+LSM” method is far less than that of AIA method. Thus, the “PCA+LSM” method is more suitable than the AIA and PCA methods for processing of randomly phase-shifted interferogram in practical experiments.

## 6. Conclusion

In conclusion, we have presented a method to extract phase from randomly phase-shifted interferogram by combining PCA and LSM. The method is fast, correct and accurate. Simulated and experimental results demonstrate the effectiveness of the proposed method. The phase-shift extraction error and residual phase error of the proposed method decrease with the increasing of the number of fringes in interferogram. The proposed method works well if the number of fringes in interferogram is larger than 1, the number of frames used is no less than 3 and the amplitude of random phase shifts is less than 1. The simulation shows that the proposed method works better than Gao’s method in phase-shift extraction. Experiment shows that the proposed method can determine the correct global phase sign and reduce the residual phase error. Thus, the proposed method is more suitable than AIA and PCA methods for randomly phase-shifted interferogram in practical experiments.

## Acknowledgments

This research is supported by the Natural Science Foundation of Zhejiang Province, China (Y1110125). We would like to express our gratitude to Javier Vargas for invaluable help and useful discuss. Furthermore, we are grateful to the reviewers for their helpful comments and suggestions.

## References and links

**1. **L. L. Deck, “Suppressing phase errors from vibration in phase-shifting interferometry,” Appl. Opt. **48**(20), 3948–3960 (2009). [CrossRef] [PubMed]

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

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

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

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

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

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

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

**9. **K. Hibino, B. F. Oreb, D. I. Farrant, and K. G. Larkin, “Phase-shifting algorithms for nonlinear and spatially nonuniform phase shifts,” J. Opt. Soc. Am. A **14**(4), 918–930 (1997). [CrossRef]

**10. **K. Okada, A. Sato, and J. Tsujiuchi, “Simultaneous calculation of phase distribution and scanning phase shift in phase shifting interferometry,” Opt. Commun. **84**(3-4), 118–124 (1991). [CrossRef]

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

**12. **B. Kong and S.-W. Kim, “General algorithm of phase-shifting interferometry by iterative least-squares fitting,” Opt. Eng. **34**(1), 183–188 (1995). [CrossRef]

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

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

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

**16. **H. Guo, Y. Yu, and M. Chen, “Blind phase shift estimation in phase-shifting interferometry,” J. Opt. Soc. Am. A **24**(1), 25–33 (2007). [CrossRef] [PubMed]

**17. **A. Dobroiu, P. C. Logofatu, D. Apostol, and V. Damian, “Statistical self-calibrating algorithm for three-sample phase-shift interferometry,” Meas. Sci. Technol. **8**(7), 738–745 (1997). [CrossRef]

**18. **A. Dobroiu, D. Apostol, V. Nascov, and V. Damian, “Self-calibrating algorithm for three-sample phase-shift interferometry by contrast leveling,” Meas. Sci. Technol. **9**(5), 744–750 (1998). [CrossRef]

**19. **A. Dobroiu, D. Apostol, V. Nascov, and V. Damian, “Statistical self-calibrating algorithm for phase-shift interferometry based on a smoothness assessment of the intensity offset map,” Meas. Sci. Technol. **9**(9), 1451–1455 (1998). [CrossRef]

**20. **L. Z. Cai, Q. Liu, and X. L. Yang, “Generalized phase-shifting interferometry with arbitrary unknown phase steps for diffraction objects,” Opt. Lett. **29**(2), 183–185 (2004). [CrossRef] [PubMed]

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

**22. **X. Xian-Feng, C. Lu-Zhong, W. Yu-Rong, and L. Dai-Lin, “Accurate phase shift extraction for generalized phase-shifting interferometry,” Chin. Phys. Lett. **27**(2), 024215 (2010). [CrossRef]

**23. **X. F. Meng, X. Peng, L. Z. Cai, A. M. Li, J. P. Guo, and Y. R. Wang, “Wavefront reconstruction and three-dimensional shape measurement by two-step dc-term-suppressed phase-shifted intensities,” Opt. Lett. **34**(8), 1210–1212 (2009). [CrossRef] [PubMed]

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

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

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

**27. **X. Chen, M. Gramaglia, and J. A. Yeazell, “Phase-shift calibration algorithm for phase-shifting interferometry,” J. Opt. Soc. Am. A **17**(11), 2061–2066 (2000). [CrossRef] [PubMed]

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

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

**30. **J. Vargas, J. A. Quiroga, and T. Belenguer, “Analysis of the principal component algorithm in phase-shifting interferometry,” Opt. Lett. **36**(12), 2215–2217 (2011). [CrossRef] [PubMed]

**31. **J. Xu, L. Sun, Y. Li, and Y. Li, “Principal component analysis of multiple-beam Fizeau interferograms with random phase shifts,” Opt. Express **19**(15), 14464–14472 (2011). [CrossRef] [PubMed]

**32. **http://en.wikipedia.org/wiki/Principal_component_analysis.

**33. **http://en.wikipedia.org/wiki/Singular_value_decomposition.

**34. **http://en.wikipedia.org/wiki/Karhunen%E2%80%93Lo%C3%A8ve_theorem.