## Subspace-based method for phase retrieval in interferometry

Optics Express, Vol. 13, Issue 4, pp. 1240-1248 (2005)

http://dx.doi.org/10.1364/OPEX.13.001240

Acrobat PDF (492 KB)

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

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]

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]

*a priori*[2

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

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]

*π*radians. The paper also investigates the influence of white Gaussian noise on the performance of the algorithm [8

8. C. Rathjen, “Statistical properties of phase-shift algorithms,” J. Opt. Soc. Am. A **12**, 1997–2008 (1995). [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]

*subspace based method*.

*forward-backward*approach [9

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]

## 2. Subspace-based method

*x, y*) on the

*t*

^{th}frame is given by

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

*N*recorded phase shifted sequences. The covariance of signal

*I(t)*in Eq. (1) is defined as [10]

*E*[·] represents the expectation operator which averages over the ensemble of realizations; the terms

*σ*

^{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]

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

**can be shown to have the form**

*R*_{I}**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

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

**APAcG**=0 and since

**AP**has full column rank we have

**A**

^{c}

**G**=0. This means that sinusoidals

*noise subspace*. This can be stated as (

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

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

*n*is determined. The details on selection of

*m*and

*N*will be explained in next Section.

*κ*is unknown and can be determined [12

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]

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

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

*α*=

*π*/4,

*κ*=2, and phase

*φ*given by

*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

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]

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

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]

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

*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

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

*α*

_{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

*π*. 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

*t*=0,1,2,…

*m*,..,

*N*-1

*κ*=1 and rewrite Eq. (A1) as

*I**(

*t-p*) for

*κ*=1 as

*c*

_{1}=

*I*

_{dc}

*a1e-*

^{iφ}

*e-*

_{iαt}+

*I*

_{dc}

*a1eiφe*

^{iαt},

*c*

_{2}=

*I*

_{dc}

*a1e-iφe-*

^{iαt}

*e*

^{−iφt}

*e*

^{−iαt}, and

*c*

_{3}=

*I*

_{dc}

*a1e-iφe*

^{iαt}+

*e*

^{2iφ}

*e*

^{2iαt}.

**E**{

*c*

_{1}}=

*E*{

*c*2}=

*E*{

*c*

_{3}}=

*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

*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

## References and Links

1. | Y. Surrel, “Phase stepping: a new self-calibrating algorithm,” Appl. Opt. |

2. | Y. Surrel, “Design of algorithms for phase measurements by the use of phase stepping,” Appl. Opt. |

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 |

4. | K. G. Larkin and B. F. Oreb, “Design and assessment of symmetrical phase-shifting algorithms,” J. Opt. Soc. Am. A |

5. | A. Patil, R. Langoju, and P Rastogi, “An integral approach to phase shifting interferometry using a super-resolution frequency estimation method,” Opt. Express |

6. | R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” in |

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 |

9. | B. D. Rao and K. V. S. Hari, “Weighted subspace methods and spatial smoothing: analysis and comparison,” IEEE Transactions on Signal Processing |

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 |

11. | A. J. Barabell, “Improving the resolution performance of eigenstructure-based direction-finding algorithms,” in |

12. | J. J. Fuchs, “Estimating the number of sinusoids in additive white noise,” IEEE Transactions on Acoustics, Speech, and Signal Processing |

**OCIS Codes**

(120.3180) Instrumentation, measurement, and metrology : Interferometry

(120.5050) Instrumentation, measurement, and metrology : Phase measurement

**ToC Category:**

Research Papers

**History**

Original Manuscript: December 24, 2004

Revised Manuscript: December 24, 2004

Published: February 21, 2005

**Citation**

Abhijit Patil and Pramod Rastogi, "Subspace-based method for phase retrieval in interferometry," Opt. Express **13**, 1240-1248 (2005)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-4-1240

Sort: Journal | Reset

### References

- Y. Surrel, “Phase stepping: a new self-calibrating algorithm,” Appl. Opt. 32, 3598-3600 (1993). [CrossRef] [PubMed]
- Y. Surrel, “Design of algorithms for phase measurements by the use of phase stepping,” Appl. Opt. 35, 51-60 (1996). [CrossRef] [PubMed]
- 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]
- 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]
- 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), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-20-4681.">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-20-4681</a> [CrossRef] [PubMed]
- R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” in Proceedings RADC, Spectral Estimation Workshop, Rome, NY, (243-258) 1979.
- 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)
- C. Rathjen, “Statistical properties of phase-shift algorithms,” J. Opt. Soc. Am. A 12, 1997-2008 (1995). [CrossRef]
- 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]
- T. Söderström and P. Stoica, “Accuracy of higher-order Yule-Walker methods for frequency estimation of complex sine waves,” IEEE Proceedings-F 140, 71-80 (1993).
- 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).
- 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]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.