Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Noise in phase shifting interferometry

Open Access Open Access

Abstract

We present a theoretical analysis to estimate the amount of phase noise due to noisy interferograms in Phase Shifting Interferometry (PSI). We also analyze the fact that linear filtering transforms corrupting multiplicative noise in Electronic Speckle Pattern Interferometry (ESPI) into fringes corrupted by additive gaussian noise. This fact allow us to obtain a formula to estimate the standard deviation of the noisy demodulated phase as a function of the spectral response of the preprocessing spatial filtering combined with the PSI algorithm used. This phase noise power formula is the main result of this contribution.

©2009 Optical Society of America

1. Introduction

Performance of Phase Shifting Interferometry (PSI) in the presence of noise is important to estimate the average noise-power of the demodulated phase. In the past several attempts have been made to analyze the effect of noisy interferograms over the estimated phase in PSI [1-5]. Despite of these efforts to our view, no general enough study has been given; in one way or another these authors [1-5] assume particular PSI algorithms to estimate their signal to noise ratios and no mention is made on the influence of the interferograms’ preprocessing. The aim of this paper is far more general; our purpose is to analyze the amount of phase noise that results from applying a spatial preprocessing combined with a PSI algorithm (as actually occurs in practice). We finally obtain a simple yet useful general formula to estimate the signal to noise ratio of the estimated phase as a function of both; the preprocessing spatial-filtering combined with the PSI algorithm (temporal-filtering) used.

The presentation of this paper is as follows: in section 2 we discuss the fact that linear filtering an interferogram having multiplicative noise gives as a result an interferogram corrupted by additive gaussian noise. Section 3 reviews temporal PSI interferometry and its spectral response. Section 4 analyzes how the combination of spatial and temporal PSI filtering determines the average noise power of the output analytical signal. Section 5 analyzes the output analytical signal as phasor to study the influence of the quadrature noise in the estimated phase and derive the main result of this paper. Finally in section 6 we draw some conclusions.

2. Spatial linear filtering of ESPI interferograms

In interferometry we must mainly deal with two kinds of noise [6]. One is the additive noise where the observed irradiance is corrupted by additive uncorrelated noise. The other one is the multiplicative phase (speckle) noise. However all interferograms contain a mixture of these two extreme cases. Interferograms corrupted mainly by additive noise arises in optical metrology of smooth surfaces such as mirrors where most of the noise comes from ambient and the electronic equipment used. On the other hand, interferograms corrupted mainly by speckle or phase noise normally arises in TV speckle metrology or ESPI [6].

As mentioned, the interferogram intensity I 0(x,y) is normally corrupted by a mixture of speckle noise ns(x,y) and additive noise na(x,y) according to the following model,

I0(x,y)=a(x,y)+b(x,y)cos[ns(x,y)+ϕ(x,y)]+na(x,y).

Where a(x,y) is the two dimensional (2D) background illumination and b(x,y) is the contrast of the fringes and finally ϕ(x,y) is the spatial modulating phase that we desire to estimate. The spatial domain (x,y) of the interferograms will be a regular 2D lattice of size K×L pixels. When dealing with smooth-surface interferometric metrology the main source of noise is additive electronic noise. On the other hand in ESPI metrology the main source of noise is multiplicative. The reason is that, although the speckle is a phase noise, we normally require two correlated speckle patterns to obtain fringes modulated by the phase being measured. These correlated fringe patterns may be modeled as,

I1(x,y)=a(x,y)+b(x,y)cos[ns(x,y)]
I2(x,y)=a(x,y)+b(x,y)cos[ns(x,y)+ϕ(x,y)].

The speckle noise ns(x,y) is a random phase uniformly distributed within [0,2π] and ϕ(x,y) is the deterministic smooth phase being estimated. The most usual form to obtain fringes phase modulated by ϕ(x,y) is to find the absolute subtraction (correlation by subtraction) of the two preceding speckle patterns as,

I(x,y)=I1I2=2bcos[ns]cos[ns+ϕ]=2bsin(ns+ϕ/2)sin[ϕ2].

Where the spatial dependence (x,y) has been omitted for clarity purposes. The fundamental harmonic in Eq. (3) of the signal ∣sin[f(x,y)/2]∣ is something proportional to a+cos[ϕ(x,y)] (a>1) and this signal is corrupted by the multiplicative speckle noise ∣sin(ns+ϕ/2)∣.

Given that our interferograms (without carrier) are low frequency signals, the first step in processing any set of phase shifted fringe patterns (whether having additive and/or multiplicative noise) is to low pass filter them to improve our signal to noise ratio, specially in ESPI fringes. This low pass filtering is typically done using (quite often repeatedly) a small 3×3 averaging convolution mask. This process however simple gives us a highly improved interferometric data. Even an ESPI interferogram convolved with a single 3×3 moving-average spatial-filter will remove pretty much multiplicative noise. Additionally to this, there is also a transformation of the noise statistics (a very welcome bonus) which transforms the multiplicative noise into additive gaussian noise. The process of changing the statistical properties of the noise has its theoretical foundations in the Central Limit Theorem (CLT) [7]. For our purposes this theorem says that a random process obtained as a linear filtering of another random process with finite mean and variance will give a filtered random signal having gaussian statistics no matter which statistical distribution the input process had. In practical ESPI metrology convolving I(x,y) (Eq. (3)) with a single 3×3 averaging window may be considered a minimum (for all practical purposes), to obtain a filtered interferogram corrupted by gaussian additive noise. Additional convolutions with this 3×3 averaging kernel will turn the output noise “more gaussian and additive”. In mathematical terms,

If(x,y)=(ξ,η)Rh1(xξ,yη)I(ξ,η)=a(x,y)+b(x,y)cos[ϕ(x,y)]+n(x,y).

Where If(x,y) is the spatial filtered ESPI image, the coefficients h 1(x,y) are the low-pass filter impulse response used to preprocess the fringe pattern and the two dimensional region R is the spatial support of h 1(x,y). The subscript a in the additive noise n(x,y) is dropped because from now on only additive noise will be considered. Passing your 3×3 averaging filter an even number of times the fringe’s modulating phase do not change. That is because the 3×3 averaging kernel has real-valued frequency response, and even number of “passes” of this filter have a real-valued and positive frequency response, so no additional phase will be introduced by the filter at any spatial frequency.

3. Symmetric quadrature filtering of phase shifted interferograms

Now that we have already filtered out most of the noise of our interferograms through spatial filtering with the “bonus” of transforming the multiplicative noise into additive gaussian noise, the next step is to estimate its modulating phase ϕ(x,y). The most useful and widely used way to estimate the modulating phase of these interferograms is Phase Shifting Interferometry or PSI. A symmetric set of N+1 phase shifted interferograms with additive corrupting noise may be modeled by

I(x,y,t)=k=N/2N/2{a+bcos[ϕ+t]+n}δ(t),α=2π/N.

Where the spatial 2D dependence of a, b, ϕ, and n has been omitted for clarity purposes. In Eq. (5) we have N+1 phase steps α uniformly distributed over the unit circle. This set of interferograms is filtered by a quadrature filter having the following symmetrical form,

h2(t)=k=N/2N/2hr(k)δ(t)+ik=N/2N/2hi(k)δ(t),α=2π/N.

In order to obtain a temporal quadrature filters from this model the coefficients of the real part hr(k) must be symmetric and the ones of the imaginary part hi(k) must be asymmetric, that is,

hr(k)=hr(k),hi(k)=hi(k)andhi(0)=0.

Just as an example, a 3-step algorithm with a phase shift equal to α=π/2 may be written as:

h3(t)=[2δ(t)δ(tπ/2)δ(t+π/2)]+i[δ(tπ/2)δ(t+π/2)].

Continuing with the general form of our complex-temporal filter model given in Eq. (6) we now calculate its Fourier transform obtaining its frequency response (spectrum) as,

H2(ωt)=2k=N/2N/2hr(k)cos(ωt)2k=N/2N/2hi(k)sin(ωt),α=2π/N.

It is important to note that the frequency spectrum H 2(ωt) is real even though h 2(t) is complex. That is, symmetrical PSI algorithms are complex in time but real in the frequency domain. Spectral analysis of PSI algorithms was presented in [11] however in a less compact way.

Just for the sake of completeness before leaving this section, we show how to calculate the phase from our symmetric temporal quadrature filter Eq. (6) [8]. Our set of N+1 interferograms convolved with Eq. (6) has a support of 2N+1 pixel. The complex output is

Ic(x,y,t)=If(x,y,t)*hr(t)+iIf(x,y,t)*hi(t)=[I(x,y,t)**h1(x,y)]*h2(t).

As we see, the analytic signal Ic(x,y,t) is the result of two steps, one is the convolution with a 2D filter h 1(x,y) followed by h 2(t). From these 2N+1 temporal signals, we are only interested with the one at t=0, which is the one that have the best phase estimation, which is

ϕn(x,y,0)=If(x,y,t)*hi(t)If(x,y,t)*hr(t)t=0=[I(x,y,t)**h1(x,y)]*hi(t)[I(x,y,t)**h1(x,y)]*hr(t)t=0.

4. Noise power of the output complex space-temporal signal Ic(x, y, t)

Let us apply the basics of random process theory in linear systems dealing with additive noise [7] to our interferometry case. A noisy set of space-temporal interferograms In(x,y,t) corrupted by additive gaussian noise n(x,y,t) may be modeled by

In(x,y,t)=I(x,y,t)+n(x,y,t).

Passing this signal through a space-temporal linear filter h(x,y,t) we obtain the filtered signal

Ic(x,y,t)=In(x,y,t)***h(x,y,t).

Where the symbol *** means a triple convolution and the analytic signal Ic(x,y,t) is obtained passing our original (noisy) interferograms by the low-pass spatial filter h 1(x,y) followed by our quadrature temporal filter h 2(t). Taking the Fourier transform this equation one obtains

Ic(ωx,ωy,t)=In(ωx,ωy,ωt)H(ωx,ωy,ωt).

Where the standard result that a convolution of two functions equals its Fourier product was used. We also assume that our additive noise n(x,y,t) is white and has a power density of η/2 watts/hertz. Now the average noise power of the complex signal Ic(x,y,t) equals [9],

σn2=E(n2)=Rnn(0,0,0)=η2∫∫∫(π,π)×(π,π)×(π,π)H2(ωx,ωy,ωt)xyt.

Where σ 2 n=E(n 2) is the square of the noise deviation (average noise power) which also equals the autocorrelation function of the random process Rnn(0,0,0) at the origin [7]. Given that our space-temporal filter h(x,y,t) may be written as a convolution between the spatial filter h 1(x,y)δ(t) and the PSI algorithm h 2(t) as h(x,y,t)= h 1(x,y)δ(t)*h 2(t), the above integral may be expressed as the following product

σn2=E(n2)Rnn(0,0,0)=η2{∫∫(π,π)×(π,π)H12(ωx,ωy)dωxy}{(π,π)H22(ωt)dωt}.

The objective to separate Eq. (14) as a product of two integrals is to clearly point out the contribution to the final noise power of the spatial filtering combined with the temporal filtering or PSI algorithm used. The first integral in Eq. (15) depends only on the spectrum of h 1(x,y) and the second integral depends only on the spectrum of h 2(t) (the PSI algorithm). Therefore in the noise reduction process, both filters the spatial and the temporal are equally important in reducing the noise power of the output analytic signal Ic(x,y,t).

 figure: Fig. 1.

Fig. 1. Frequency response of a 3 and 5-step [10] PSI algorithms both having a phase step of π/2 radians. H3(ωt) represents the 3-step PSI and H5(ωt) the 5-step one, having the same response be at ωt, =-1. We also show the analytical signal (ωt+1.0)e which is the output of these filters to a+bcos[ϕ+(π/2) t]. Finally the value of the square integral of their frequency response is given to obtain a relative noise rejection.

Download Full Size | PDF

To compare the relative figure of noise attenuation among several PSI algorithms having the same output signal power b 2, we must compare the value of the average power E(n 2) among them according to Eq. (15). An easier way to compare two PSI algorithms A and B, keeping h 2(x,y) constant is the use of the following inequality,

{(π,π)HA2(ωt)dωt}<{(π,π)HB2(ωt)dωt}.

In this case the PSI algorithm A rejects more noise than algorithm B. As an example of this, Fig. 1 shows two PSI algorithms having the same response at ωt=-1 which is be . On the left side we have the frequency response of a 3-step PSI Eq. (7), and on the right side the response of a 5-step one [10], both have a phase step α=π/2. Finally the value of their square modulus is given. We can formally see the intuitive result that a 5-step PSI algorithm rejects more noise than a 3-step one; both having the same preprocessing filter h 2(x,y) and output signal.

Before leaving this section, it is very important to note that although we have focused in symmetrical N+1 PSI algorithms giving a real spectrum for H 2(ωt); in general however H 2(ωt) is complex (non-symmetrical PSI). In which case you must use [H 2(ωt)]2= H(ωt) H *(ωt).

5. Estimated phase-noise as a function of the interferograms’ processing and noise

As we saw in the preceding section the combination of the spatial low-pass filtering along with the quadrature temporal filter (the PSI algorithm) are equally responsible for the average noise power σ 2 n=E(n 2) at the output complex signal Ic(x,y,t). The noisy estimated phase of the interferograms is given as the angle ϕn(x,y,0) of the following phasor,

Ic(x,y,0)=Re{Ic(x,y,0)}+iIm{Ic(x,y,0)}.

Where the operators Re{.} and Im{.} take the real and imaginary parts of their argument. The signal in the last equation can also be represented as function of the “original” (noiseless) complex interferogram at ϕ(x,y,0) corrupted by additive-gaussian quadrature noise as,

Ic(x,y,0)=b(x,y)cos[ϕ(x,y,0)]+nr(x,y,0)+ib(x,y)sin[ϕ(x,y,0)]+ni(x,y,0).

Where nr(x,y,0) and ni(x,y,0) are the quadrature components of the noise n(x,y,0) [9]. Figure 2 shows a phasorial representation of all the interesting parameters involved in this paper. Having a fringe contrast b and looking at Fig. 2 we see that from the triangle having as base σn and “height” b the standard deviation of the phase noise is,

 figure: Fig. 2.

Fig. 2. Phasor representation (spatial and time dependences were omitted) of the estimated analytical signal Ic, and its corrupting noise n. The noiseless amplitude of Ic is represented by the phasor b which has a noiseless angle ϕ. The noise n has an average power σn 2 given by Eq. (15) and its angle is uniformly distributed within [0.2π], so the phasor n may point anywhere within the circle shown. The noisy estimated phase ϕn is obtained by Eq. (10).

Download Full Size | PDF

σϕn=2tan1(σn2b).

Where σn 2 is the total average noise power given in Eq. (15). Using the approximation tan(x)~x, and Eq. (15), equation nineteen may be re-written as,

σϕn2σn2b2=η2b2{∫∫(0,2π)×(0,2π)H12(ωx,ωy)dωxy}{(0,2π)H22(ωt)dωt}

This equation is the main result of this paper, and confirms the intuitive facts that; 1) the corrupting noise at the output of the spatial-PSI filtering (Eq. (15)) is gaussian due to the linear filtering; 2) as a consequence the output phase noise is also a gaussian random process with the noiseless phase ϕ as its mean; 3) the noise power σϕn 2 of the noisy estimated phase ϕn is proportional to the output noise power σn 2 of our linear filters and inversely proportional to the interferogram’s contrast b; and 4) the amount of noise on the estimated phase depends equally on the spatial low-pass filter and on the PSI algorithm used. The condition bn, used to derive Eq. (19) is for all practical purposes always fulfilled; otherwise the recovered phase would be too noisy and the physical quantity under measurement would be highly unreliable.

6. Conclusion

We have analyzed the influence of the noise in temporal interferometry. We have shown the shared weight between the preprocessing spatial filtering and the temporal PSI algorithm in estimating the fringe’s phase of a set of noisy phase shifted interferograms. Equation (20) states the main result of this paper which is an estimate of the standard deviation of the noisy demodulated phase due to a particular combination of a spatial and a PSI algorithm. We have also formally shown in Eq. (20 ) the long standing intuitive idea that having a set of good spatially filtered interferograms along with a PSI algorithm having many phase shifted interferograms (narrow-bandwidth) one would necessarily obtain a good estimated phase.

Acknowledgments

We acknowledge the Mexican National Science and Technology Council (Consejo Nacional de Ciencia y Tecnologia, CONACYT) for its valuable support.

References and links

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

2. 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).

3. Y. Surrel, “Additive noise effect in digital phase detection,” Appl. Opt. 36, 271–276 (1997). [PubMed]  

4. J. Schmit and C. Katherine, “Window function influence on phase error in phase-shifting algorithms,” Appl. Opt. 35, 5642–5649 (1996). [PubMed]  

5. G. Paez and M. Strojnik, “Analysis and minimization of noise effects in phase shifting interferometry,” SPIE Vol. 3744, 295–305 (1999).

6. K. J. Gasvik, Optical Metrology (John Wiley & Sons Ltd); 2th ed., (1996).

7. A. Papoulis, Probability, Random Variables, and Stochastic Processes, (McGraw-Hill, 3th ed., 1991).

8. M. Servin and M. Kujawinska, Modern Fringe Pattern Analysis in Interferometry, in Handbook of Optical Engineering, D. Malacara and B. J. Thompson eds., (Marcel Dekker Inc., 2001) Chap. 12.

9. L. W. Couch, Digital & Analog Communication Systems, (Prentice Hall, 2006) 7th ed.

10. P. Hariharan, B. F. Oreb, and T. Eyui, “Digital phase shifting interferometry: a simple error-compensating phase calculation algorithm,” Appl. Opt. 26, 2504–2505 (1987). [PubMed]  

11. K. Freischland and C. L. Koliopolous, “Fourier description of digital phase-measuring interferometry,” J. Opt. Soc. Am. A 7, 542–552 (1990).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (2)

Fig. 1.
Fig. 1. Frequency response of a 3 and 5-step [10] PSI algorithms both having a phase step of π/2 radians. H3(ωt ) represents the 3-step PSI and H5(ωt ) the 5-step one, having the same response be at ωt , =-1. We also show the analytical signal (ωt +1.0)e which is the output of these filters to a+bcos[ϕ+(π/2) t]. Finally the value of the square integral of their frequency response is given to obtain a relative noise rejection.
Fig. 2.
Fig. 2. Phasor representation (spatial and time dependences were omitted) of the estimated analytical signal Ic, and its corrupting noise n. The noiseless amplitude of Ic is represented by the phasor b which has a noiseless angle ϕ. The noise n has an average power σn 2 given by Eq. (15) and its angle is uniformly distributed within [0.2π], so the phasor n may point anywhere within the circle shown. The noisy estimated phase ϕn is obtained by Eq. (10).

Equations (22)

Equations on this page are rendered with MathJax. Learn more.

I0(x,y)=a(x,y)+b(x,y)cos[ns(x,y)+ϕ(x,y)]+na(x,y).
I1(x,y)=a(x,y)+b(x,y)cos[ns(x,y)]
I2(x,y)=a(x,y)+b(x,y)cos[ns(x,y)+ϕ(x,y)].
I(x,y)=I1I2=2bcos[ns]cos[ns+ϕ]=2bsin(ns+ϕ/2)sin[ϕ2].
If(x,y)=(ξ,η)Rh1(xξ,yη)I(ξ,η)=a(x,y)+b(x,y)cos[ϕ(x,y)]+n(x,y).
I(x,y,t)=k=N/2N/2{a+bcos[ϕ+t]+n}δ(t),α=2π/N.
h2(t)=k=N/2N/2hr (k) δ (t)+i k=N/2N/2hi(k)δ(t),α=2π/N .
hr(k)=hr(k),hi(k)=hi(k)andhi(0)=0.
h3(t)=[2δ(t)δ(tπ/2)δ(t+π/2)]+i[δ(tπ/2)δ(t+π/2)].
H2(ωt)=2k=N/2N/2hr (k) cos(ωt) 2k=N/2N/2hi(k)sin(ωt),α=2π/N.
Ic(x,y,t)=If(x,y,t)*hr(t)+iIf(x,y,t)*hi(t)=[I(x,y,t)**h1(x,y)]*h2 (t) .
ϕn(x,y,0)=If(x,y,t)*hi(t)If(x,y,t)*hr(t)t=0=[I(x,y,t)**h1(x,y)]*hi(t)[I(x,y,t)**h1(x,y)]*hr(t)t=0.
In(x,y,t)=I(x,y,t)+n(x,y,t).
Ic(x,y,t)=In(x,y,t)***h(x,y,t).
Ic(ωx,ωy,t)=In(ωx,ωy,ωt)H(ωx,ωy,ωt).
σn2=E(n2)=Rnn(0,0,0)=η2∫∫∫(π,π)×(π,π)×(π,π)H2(ωx,ωy,ωt)xyt.
σn2=E(n2) Rnn (0,0,0) = η2 {∫∫(π,π)×(π,π)H12(ωx,ωy)dωxy} {(π,π)H22(ωt)dωt}.
{(π,π)HA2(ωt)dωt}<{(π,π)HB2(ωt)dωt} .
Ic(x,y,0)=Re {Ic(x,y,0)}+iIm{Ic(x,y,0)} .
Ic(x,y,0)=b(x,y)cos[ϕ(x,y,0)]+nr(x,y,0)+ib(x,y)sin[ϕ(x,y,0)]+ni(x,y,0).
σϕn=2tan1(σn2b).
σϕn2σn2b2=η2b2{∫∫(0,2π)×(0,2π)H12(ωx,ωy)dωxy}{(0,2π)H22(ωt)dωt}
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.