The general theory of phase shifting algorithms

We have been reporting several new techniques of analysis and synthesis applied to Phase Shifting Interferometry (PSI). These works are based upon the Frequency Transfer Function (FTF) and how this new tool of analysis and synthesis in PSI may be applied to obtain very general results, among them; rotational invariant spectrum; complex PSI algorithms synthesis based on simpler first and second order quadrature filters; more accurate formulae for estimating the detuning error; output-power phase noise estimation. We have made our cases exposing these aspects of PSI separately. Now in the light of a better understanding provided by our past works we present and expand in a more coherent and holistic way the general theory of PSI algorithms. We are also providing herein new material not reported before. These new results are on; a well defined way to combine PSI algorithms and recursive linear PSI algorithms to obtain resonant quadrature filters. ©2009 Optical Society of America OCIS codes: (120.3180) Interferometry; (120.2650) Fringe Analysis. References and links 1. K. Freischlad, and C. L. Koliopoulos, “Fourier description of digital phase-measuring interferometry,” J. Opt. Soc. Am. A 7(4), 542–551 (1990). 2. D. W. Phillion, “General methods for generating phase-shifting interferometry algorithms,” Appl. Opt. 36(31), 8098–8115 (1997). 3. Y. Surrel, “Design of algorithms for phase measurements by the use of phase stepping,” Appl. Opt. 35(1), 51–60 (1996). 4. D. Malacara, M. Servin, and Z. Malacara, Interferogram analysis for Optical Testing, 2th ed., (Marcel Deker, 2003). 5. M. Servin, J. C. Estrada, and J. A. Quiroga, “Spectral analysis of phase shifting algorithms,” Opt. Express 17(19), 16423–16428 (2009). 6. J. G. Proakis, and D. G. Manolakis, Digital Signal Processing, 4th-ed., (Prentice Hall, 2007). 7. J. Schmit, and K. Creath, “Extended averaging technique for derivation of error-compensating algorithms in phase-shifting interferometry,” Appl. Opt. 34(19), 3610–3619 (1995). 8. J. F. Mosiño, M. Servin, J. C. Estrada, and J. A. Quiroga, “Phasorial analysis of detuning error in temporal phase shifting algorithms,” Opt. Express 17(7), 5618–5623 (2009). 9. M. Servin, J. C. Estrada, J. A. Quiroga, J. F. Mosiño, and M. Cywiak, “Noise in phase shifting interferometry,” Opt. Express 17(11), 8789–8794 (2009). 10. J. C. Estrada, M. Servin, and J. A. Quiroga, “Easy and straightforward construction of wideband phase-shifting algorithms for interferometry,” Opt. Lett. 34(4), 413–415 (2009). 11. K. G. Larkin, and B. F. Oreb, “Propagation of errors in different phase-shifting algorithms: a special property of the arctangent function,” presented at the SPIE International Symposium on Optical Applied Science and Engineering, San Diego, California, SPIE, 1755, 219–227 (1992). 12. F. G. Stremler, Introduction to Communications Systems, 3rd ed., (Addison-Wesley, 1990). 13. 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,” Appl. Opt. 13(11), 2693–2703 (1974). 14. K. Hibino, “Susceptibility of systematic error-compensating algorithms to random noise in phase-shifting interferometry,” Appl. Opt. 36(10), 2084–2093 (1997). 15. C. J. Morgan, “Least-squares estimation in phase-measurement interferometry,” Opt. Lett. 7(8), 368–370 (1982). 16. V. K. Madisetti, and D. B. Williams, eds., Digital Signal Processing Handbook, (CRC Press, IEEE Press, 1998). #116619 $15.00 USD Received 1 Sep 2009; revised 28 Oct 2009; accepted 29 Oct 2009; published 13 Nov 2009 (C) 2009 OSA 23 November 2009 / Vol. 17, No. 24 / OPTICS EXPRESS 21867


Introduction
Spectral analysis in PSI algorithm theory is systematically made since 1990, and it is based on the work of Freishlad and Koliopoulos (F&K) [1].Application of the F&K spectral analysis to PSI algorithm synthesis may be seen in [2].In 1996 Surrel [3] developed an algebraic approach to analyze PSI algorithms based on what he called the characteristic polynomial associated to the quadrature filter.Even though the properties of any PSI algorithms may be deduced from the roots of this polynomial, no spectral plot (similar to the one proposed in [1]) where one may gauge in a glance the frequency response of PSI algorithms was proposed.Over the years, the F&K spectral plot emerged as the standard way of gauging the frequency response of PSI algorithms [4].However the F&K spectral analysis has a major drawback, namely: the spectrum varies when the PSI algorithm's reference signal is time-shifted.However it is well known that the estimated phase of a reference time-shifted PSI algorithm remains unchanged, except for an irrelevant piston [4, page 239].In reference [5] we have recently proposed a new way to analyze the spectra of PSI algorithms based on the Frequency Transfer Function (FTF).This new way of spectral analysis of PSI algorithms is however not new in engineering.In particular in electrical engineering the use of the FTF has been the standard way of analyzing the spectra of digital filters for decades now [6].As we show in this paper, this new way of analyzing the spectra in PSI interferometry is very useful and productive for the analysis, synthesis and extensions of PSI algorithms.
In this paper we start by reviewing the basic problem faced in few steps PSI.We continue by describing this new (in PSI) way of spectral analysis based on the FTF [5].We then apply this new spectral tool to illustrate how this powerful analysis gives more straightforward formulae for different analyzing aspects of PSI algorithms.In particular we analyze in a more general way than before several aspects of PSI filtering.These are 1) reference time-shift invariant spectral analysis of PSI algorithms; 2) close form formulae to estimate the PSI algorithm's detuning error; 3) general formulae to estimate the output-power phase noise of PSI algorithms; 4) a well defined way to combine several PSI algorithms in such a way that the composed PSI algorithm inherits the good properties of its components, 5) synthesis of arbitrarily complex PSI algorithms based on much simpler building blocks; first and second order filters; and finally 6) how recursive linear quadrature filters may be used as "fast" PSI algorithms.We are not only reviewing and expanding on these subjects (themes 1, 2, 3, and 5).But we are also offering original research on; a well defined way to combine PSI algorithms (theme 4), and (for the first time in the PSI literature) the use of recursive quadrature filters as PSI demodulators (theme 6).
In the past, several methods to "combine" new algorithms have been reported (see [7] and its references).The intention of combining two PSI algorithms to obtain a third one is that in some way the offspring algorithm is expected to inherit the good features of its parents.This is a useful technique to obtain more complex PSI algorithms from simpler ones.However the techniques used are not systematic, nor well defined [7].In particular it is not easy to see which and how these good features were inherited.With our new approach herein presented, based on the FTF the combining problem becomes clear and much better defined.
Up to this date all PSI algorithms [4] are quadrature filters having a Finite Impulse Response (FIR).That is, their quadrature filter's impulse response is composed only of a finite number of samples.The "size" of these quadrature filters is normally the size of our interferometric sampled data.In this way we colloquially talk about a 5-step PSI algorithm when we have five temporal interferograms.Implicitly in this way we assume that we are using a quadrature digital filter having an impulse response composed of only five samples.On the other hand, an Infinite Impulse Response (IIR) quadrature filter (or PSI algorithm) is composed by an infinite number of samples.This might be surprising in the PSI field, but it has a very close analog in Optics, which is (as we analyze in this paper), the Fabry-Perot resonant interferometer.We will analyze this resonant PSI algorithm and how this recursive linear system may be used at high temporal sampling rates, extracting the phase in line as the experiment proceeds.

The aim of a quadrature filter in phase shifting interferometry
Let us begin by showing the usual mathematical model of a temporal interferometric signal as, 0 ( , , ) ( , ) ( , ) cos[ ( , ) ] .

I x y t a x y b x y
x y t Where a(x,y), b(x,y), and φ(x,y) are respectively the background, contrast and the searched phase.Finally ω 0 is the temporal carrier (in radians) of the phase shifted interferograms.
Rewriting it as, The spatial dependence of the functions a(x,y), b(x,y) and φ(x,y) were omitted for clarity.
The aim of a PSI algorithm is to filter out (using a quadrature filter h(t)) a(x,y), and one of the two complex exponentials.Therefore the desired output complex signal Where h(t) is the filter's impulse response, and ] its Fourier transform.The symbol * denotes a one-dimensional temporal convolution.The function h(t) is the PSI filter's complex impulse response given as,

( ) ( ) ( ) . h t hr t i hi t
The minimal conditions on the frequency transfer function H(ω) to obtain the complex signal Ic(t) in Eq. ( 3) from the interferogram are: This equation states that in the frequency domain we need to filter out (at least) the background signal at ω = 0, and one of the two complex signal, in this case the one centered at ω = ω 0 .There exist an infinite number of possible filters which comply with the conditions stated by Eq. (5).For example in Fig. 1 we graphically illustrate two such possibilities for H(ω).The input of these linear systems is the real-valued interferogram I(t) and its output is the searched complex signal Ic(t).5).The input to these filters is the interferogram I(t) and the output the searched analytical signal Ic(t).These two FTFs are valid PSI algorithms with temporal carrier ω0.Both have the minimum required zeroes for H(ω), these are H(0) = 0, and H(ω0) = 0.
If the detector has non-linear behaviour or it is saturated (above or below) then higher order harmonics are generated in the interferogram data.In this case one may extend the minimal zeroes in Eq. (5) to include any desired higher n-order harmonic in the condition H(nω 0 ) = H(nω 0 ) = 0, where n>1.These extended conditions will of course require a higher order PSI algorithm to meet the additional required zeroes.In section 9 in this paper we show how to synthesize higher order PSI filters to meet with these and/or other design criteria.

The general mathematical form of a PSI algorithm (quadrature filter)
In this section we will assume that the interferometric data and the filter h(t) are centered around t = 0. Also we are assuming a set of N phase shifted interferograms as data, and that the digital filter h(t) also has N samples.
The most general impulse response of a N-step, linear-quadrature digital-filter (PSI algorithm) is Where T is the digital filter's sampling rate, the coefficients a n , are N weighting (possible complex) constants, and the angular frequency of the reference signal (or local oscillator) is ω 0 .This equation is a product of two temporal functions: one is the complex periodic reference exp[iω 0 t].The other one is a set of N temporal deltas weighted by a n .By itself the weighted finite set of deltas constitutes a digital filter centered at ω = 0 (a low pass filter).The constants a n shape along with the deltas δ(t-nT), the frequency response of this filter.The complex reference is the responsible of transforming the low pass filter at ω = 0 into a quadrature (one sided) filter centered at ω = ω 0 (a PSI algorithm).Quadrature filtering a set of N temporal interferograms using the digital complex filter h(t), may be represented as, Where I s (t) represents our N sampled data.The temporal sequence I s (t)*h(t) has therefore a support of 2N samples.However one normally needs only the complex signal at the temporal origin (t = 0).Therefore the output analytical signal at t = 0 is, 0 The demodulated phase of the interferogram is the angle of this complex signal at t = 0, and it is given by This is called the PSI algorithm associated with the filter's impulse response h(t) in Eq. ( 6) [5].

The spectral analysis based on the Frequency Transfer Function (FTF)
In this section we briefly review the spectral analysis proposed by F&K [1] and in the next one we show a major drawback of it.The spectral analysis in [1] take the Fourier transform of the real hr(t) and the imaginary hi(t) components of h(t) separately,

Hi F hi t Hr
F hr t Having these Fourier transforms, F&K recommend to analyze this spectrum without constant or common phase factors [1], and in this way obtain two real functions to plot.On the other hand, we propose to analyze the spectra of the complex sum h(t) = hr(t) + i hi(t), The function H(ω) is the Frequency Transfer Function (FTF).As shown in the coming sections of this paper, the two most relevant virtues of this new PSI spectrum representation (Eq.( 11)) is that it is invariant to PSI filter's rotation (reference time-shift) and that its rejected frequencies are clearly shown as zeroes over the spectral line [5].

Reference time-shift in PSI algorithms
Our first example on the use of Eq. ( 11) (the FTF) instead of Eq. (10) to spectrally analyze PSI algorithms is the time-shift of the reference signal of an algorithm.Let us see what happens if the reference signal in h(t) is time-shifted t 0 , one obtains, From this we can see that time-shifting the reference implies filter's rotation by an angle of ∆ 0 = -ω 0 t 0 radians.This was seen by Malacara et.al. [4], page 239.Rotating h(t) by ∆ 0 one obtains, Filter's rotation linearly mix-up the real and imaginary parts of the original PSI filter h(t).These give apparently "new" PSI algorithms that look quite different albeit being the same one [7].The family of PSI algorithms obtained from the rotation ∆ 0 = -ω 0 t 0 are, According to the spectral analysis proposed by F&K [1], one must plot two real functions associated to the following Fourier transforms,

Hi F hr t hi t
Hr F hr t hi t The spectrum obtained from Eq. ( 15) is clearly different from the one in Eq. ( 10) for the same (albeit rotated) PSI algorithm!.This difference makes us wrongly think that we are dealing with two different PSI algorithms.
Most PSI algorithms introduce zero or constant phase shifts on their own.In particular, rotated PSI algorithms introduce a constant phase shift ∆(ω) = ∆ 0 [4].Therefore in most cases |H(ω)| may be used as the main analyzing function given that it is rotation invariant, this function is given by Where H(ω) = Re(ω) + i Im(ω), being Re(ω) and Im(ω) real functions of ω.Moreover, in this equation we also show the well known fact that [6]; the rotated h(t)exp[i∆ 0 ] and the original PSI filter h(t) have the same magnitude |H(ω)|.In conclusion, while the F&K spectral analysis give different spectra for the same (rotated) PSI algorithm, our analysis based on the magnitude of the FTF remains invariant to rotations.
As an example of this let us write down two (apparently different) 5-steps PSI algorithms, One may wrongly think that these two algorithms are different.According to the spectral analysis proposed by F&K they also have a different spectral plot.However it is known long ago [4,7] that these two algorithms are the same albeit rotated.In other words, the basic difference between them is a phase shift of the reference equal to π/4 radians.In Fig. 2 we show this as an example of the invariance to rotation of the FTF.The magnitude of H(ω) of these two (rotated) versions of the 5-step Schwider-Hariharan [4] algorithm is of course identical.This means that both 5-step algorithms have identical phase demodulation properties.
Fig. 2. The magnitude of the spectrum of the complex quadrature filter associated with the "two" 5-step Schwieder-Hariharan algorithms in Eq. ( 17) is identical.The spectral components of the real-valued interferogram that are rejected are the background H(0) = 0 and the complex signal at ω0, i.e.H(ω0) = 0 This spectrum remains invariant to the reference (local oscillator) time-shift.

Evaluation of the detuning error in PSI algorithms
The most common systematic source of error in PSI arises when the PSI algorithm has a reference signal at frequency ω 0 , but the actual interferometric data has a carrier of ω 0 + ∆.This frequency mismatch ∆ between the reference and the data is called detuning.An ideal one-sided quadrature filter h ideal (t) do not have such a problem.Because this ideal quadrature filter would have a FTF given by: [ ] This ideal filter will extract the analytical signal of any sequence of interferograms no matter what temporal carrier they had.However, such a sharp filter would require an infinite size PSI algorithm.
In few-step PSI, the more humble expectations over the H(ω), stated in Eq. ( 5) are searched.If our PSI algorithm is sensitive to detuning the output signal at t = 0 is obtained as, This is formed not only by the desired complex signal exp[-iφ], but also by a spurious one exp[iφ] at the other (mirror) side of the spectrum (see Fig. 3).A phasorial representation of Eq. ( 19) is given in Fig. 4. The angle of this signal is no longer φ but ψ, and this erroneous phase is given by, This expression may be simplified assuming that the spurious signal H(ω 0 + ∆)exp[iφ] is small compared to the desired one H(-ω 0 -∆)exp[-iφ], and for small detuning tan(φ−ψ) = φ−ψ.Using these assumptions one obtains the erroneous demodulated phase ψ as [8], We have used the modulus of the FTF provided that this function may be complex.Equation (21) contains the well known fact that the demodulated phase ψ obtained from the analytical signal in Eq. ( 19) equals the desired phase φ plus a spurious signal proportional to the original fringe pattern doubling its fringe number [4,8,11].What is new in Eq. ( 21) is the amplitude of this spurious signal sin(2φ) as a function of the spectrum of the PSI filter.This amplitude |H(ω 0 + ∆)|/|H(-ω 0 -∆)| is given as a function of H(ω) evaluated at the detuned frequencies at both sides of the spectrum.
To deal with higher n-order (n>1) harmonics miscalibration or detuning, the condition in Eq. ( 21) must be extended to these harmonics.In section 9 in this paper we propose a general synthesis technique that may solve this problem by generating zeroes in the FTF of the PSI filter everywhere you want in the spectral line.

Output phase noise estimation in PSI quadrature filtering
Now let us turn our attention to noisy interferograms corrupted by additive noise, and estimate the noise power of the demodulated phase.Let us start with the following model for an interferogram corrupted by additive noise at pixel (x,y), 0 ( ) cos( ) ( ) .In t a b t n t Where n(t) is a stationary zero-mean white gaussian random process with a flat spectral power density given by F[E{n(t)n(t + τ)}] = η/2 [12].Where E{.} stands for the ensemble average, and F[.] its Fourier transform.After this signal passes through a linear quadrature filter h(t) one obtains the output noisy analytic signal as, Where Φ(t) is a random process uniformly distributed in [0,2π] [12].We now have transformed our real noise n(t) into a complex-stationary band-pass gaussian-additive noise nc(t)exp[-iΦ(t)-iω 0 t].The power E{nc 2 (t)} of the output noise is given by [12] 2 2 Where H(ω) is the FTF of our PSI algorithm H(ω) = F[h(t)], and H*(ω) denotes its complex conjugate.As usual we need to know the angle of the following sum of complex signals, Before applying a PSI algorithm, one normally low-pass the fringes by means of averaging spatial convolution filters.Depending upon the amount of corrupting noise one normally passes a 3x3 averaging window several times until clearer fringes are obtained.In practice, this strong spatial low-pass filtering makes the condition (b/2)H(-ω 0 )>>nc(0) to hold.One may demonstrate [9] that with the aid of the phasorial diagram in Fig. 5, and the just mentioned inequality, the standard deviation of the output noisy phase φ n is, The circle on this diagram reefers to the fact that the complex band-pass noise's phasor nc(0) may point anywhere due to its uniformly distributed random phase Φ(0).Using the approximation, tan(x) = x and the FTF of the quadrature filter, the noisy-phase's expected value and its deviation are better displayed as, As we intuitively expected [9,12], the ensemble average of our noisy phase E{φ n } is the searched modulating phase φ.And its standard deviation E{φ n 2 } is proportional to the output noise-energy σ 2 nc and inversely proportional to the output signal energy H 2 (-ω 0 )(b 2 /4).We finally mention that a similar result of that in Eq. ( 27) was obtained by Hibino [14], applying a kind of demodulated-phase sensitivity variation to a random variation in the intensity of the interferograms.However in [14], the demodulated phase-noise analysis is not given in terms of the frequency response of the PSI algorithm as we did in this section.Fig. 6.Here we show the difference of noise rejection between the standard 3-step PSI algorithm [13] and the 5-step Schwider-Hariharan algorithm [4].For the same output signal power the integral under |H(ω)| 2 is greater for the case of a 3-step algorithm than for the 5-step one.As a consequence, the 5-step algorithm rejects more noise than the 3-step one.
As example, in Fig. 6 we show the magnitude of the FTF of two popular PSI algorithms; the 3-step one [13] and the 5-step Schwider-Hariharan [4].In this figure we show the area under |H(ω)| 2 and see how the 5-step algorithm has less output noise (area) for the same output signal amplitude.Finally in reference [9] we sketch a proof that although we have assumed additive noise, the fact of preprocessing the fringe patterns with averaging (low-pass) spatial filters makes that multiplicative noise be converted into additive noise for every situation that occurs in practice.

Combining several PSI algorithms
In the past several papers have been published (see [7] and its referenced papers) on techniques to combine two or more PSI algorithms to obtain more complex algorithms inheriting some desired properties of their components.Despite these papers' claims, these combining techniques are not systematic, nor intuitive, neither well defined.On the other hand working within the FTF framework the combination of lower order PSI algorithms to obtain a more complex PSI algorithm, follows the standard "combining by convolution" of any sequence of linear filters.
Let us start by considering for example the combination of the following two PSI algorithms, Of course properly tuned (which it is assumed) these two algorithms will estimate the same phase φ.However the algorithm properties may be different, such as detuning robustness, harmonic rejection, or noise filtering.These two PSI algorithms have their corresponding impulse response as,

( ) . h t hr t i hi t h t hr t i hi t
Well, the best way to combine them is, as most engineers do; by convolution!In linear system engineering the natural way to combine filters to obtain a "combined" filter h(t) is,

] . h t h t h t hr hr hi hi i hr hi hi hr
The t variable has been omitted for clarity.The "combined" PSI algorithm now reads, Of course the generalization of the above scheme is straightforward.As we see using the FTF one may easily "combine" in an intuitive and in a well defined way any number of previously known PSI algorithms, and obtain more complex PSI algorithms inheriting the good features of its components.The obtained algorithms have as temporal support the sum of its parents.That is combining by convolution a M-step algorithm with a N-step one, one obtains a (M + N)-step PSI algorithm.
In the next section we offer a much better technique to synthesize sophisticated PSI algorithms than the one just presented by combining existing PSI algorithms.

Synthesis of sophisticated PSI algorithms from basic building blocks
The synthesis of new PSI algorithms is a recurrent theme reported in countless number of papers [4].Researches in this field have proposed a wide range of interesting and useful PSI algorithms since the basic and best known 3-step algorithm was first published [13].Early attempts to synthesize higher order PSI algorithms were made more or less ad hoc [4].However the first systematic way to design arbitrarily long N-step PSI algorithms was the algebraic, least-squares approach [15].Although the least-squares technique gives all the required coefficients of an N-step PSI algorithm straightforwardly, its spectral properties remains unknown and uncontrollable to the PSI filter's designer.As it occurs with the leastsquares technique, all designs strategies known to the authors, have synthesized N-step algorithms as a whole.There is however an interesting alternative to this and it comes from linear system engineering [6,16].In this field, synthesis of sophisticated linear filters is typically made from much simpler subunits or building blocks.The two subunits normally used are first and second order filters.These two basic building blocks have zeroes and/or poles that can be positioned anywhere in the spectrum by the designer.In this way more interesting filters may be synthesized as a convolution product of many of these two basic components.This procedure is similar to the process presented in the preceding "combining PSI algorithms" section.In standard digital filter design [16], the required zeroes and poles are first positioned in specialized spectral-charts as the design goal.From this careful collocation of zeroes and poles one may synthesize in a straightforward way arbitrarily sophisticated filters having the desired frequency and temporal behavior.Well known examples of this design philosophy are the Butterworth (maximally flat), or the Chebyshev approximations, among many others [16].
In PSI interferometry one may use the same strategy [10] with the exception that in our case (up to this date) our PSI filters only generate zeroes.That is, there are no PSI algorithms having (real or complex) poles.This is a fortunate situation for the PSI subfield, provided that zeroes-only systems are always stable.The reason for having zeroes-only PSI filters is that we do not use recursive PSI algorithms, where poles in their FTF H(ω) may arise.In spite of this we may manage to build sophisticated PSI algorithms with any desired spectrum shape from much simpler building blocks.
Our design strategy consists of the following three simple steps: first, to construct a FTF of the desired quadrature filter by carefully collocating its zeroes anywhere we want in the spectral line; second, find the inverse Fourier transform of this FTF; finally separate the real and imaginary parts of the resulting filter and construct from it the formula for the desired PSI algorithm.We want to remark once again that, these steps follows very closely the standard way of filter design in linear digital systems engineering [6,16].
Let us start by showing our simplest building block to construct higher order PSI algorithms; this is the first order difference equation, Where (without loose of generality) we have normalized the time step to T = 1.This first order filter has the following FTF, This FTF correspond to a tunable (through ω 0 ) first-order digital-filter that "punches a zero" (rejects the signal at ω = ω 0 ) at any desired frequency ω 0 (and its harmonics ω 0 + 2mπ) in the spectral line.To have a viable PSI algorithm using Eq.(34) and according to Eq. ( 5), we would need at least two first order filters.One to reject the background at ω = 0, and another one to filter-out the unwanted analytical signal at ω = + ω 0 .Convolving these two filters one obtains the composed FTF given by, 0 ( ) 4 sin( ) sin( ) .
This FTF "punches" at least two zeroes in the spectrum; one at ω = 0, and the other at ω = + ω 0 .The explicit collocation of the desired zeroes is the first step of our synthesis process.The impulse response of this filter is the inverse Fourier transform of this equation, that is, This impulse response determination is second step in our synthesis technique.From this we find the real hr(t) and the imaginary hi(t) parts, and write down its associated PSI algorithm, ] cos( ) This is our last step in our PSI algorithm synthesis technique.Although this algorithm is not particularly interesting; it was only presented to illustrate our synthesis method.Now we present the most interesting building block in phase shifting interferometry.This is the following second order quadrature filter, 2 0 With Fourier transform (its FTF) given by, This second order system along with the first order one in Eq. (33) are our main workhorse in our synthesis procedure.The advantage of using this second order filter is its robustness to detuning at ω 0 .That is H 2 (ω 0 ) = 0 touches zero tangentially; i.e. [dH 2 (ω)/dω] = 0 at ω = ω 0 .A simple, yet useful PSI algorithm that may be constructed with this second order building block along with the first order one h 1 (t) is the following, In particular, it is interesting to note that substituting ω 0 = π/2 as the reference frequency in this equation, one recovers the well known Schwider-Hariharan 5-step PSI algorithm [4].We may in a straightforward way continue assigning first or second order zeroes over the spectral line to synthesize PSI algorithms having rejection at every desired frequency (or harmonics).Apart from punching zeroes anywhere we want using these two subunits; we know the spectral behaviour in the neighborhood of these zeroes.For example an interesting and sophisticated filter built from a first order, and three second order components is, Choosing the three tuning frequencies at ω 1 = π/4, ω 2 = π/2 and ω 3 = 3π/4 one obtains the ultra-wide-band 9-step quadrature filter [10] shown in Fig. 7.  42).This is an ultra-wide-band filter that corresponds to a 9-step PSI algorithm.This algorithm would demodulate any 9-step sequence of interferograms no matter what carrier they may have.However to obtain the maximum signal to noise ratio, according to Eq. ( 27), the best carrier is π/2 radians/frame.The "flat zone" is strictly zero only at the tuning frequencies: ω1 = π/4, ω2 = = π/2 and ω3 = 3π/4.But compared to the amount of measuring noise, the non-zero ripples in the right-hand "flat-zone" are small enough to be negligible..This quadrature filter has the nice properties (detuning robustness) of the Schwider-Hariharan 5-step algorithm repeated at frequencies ω 1 , ω 2 and ω 3 .This filter does not need to know the temporal carrier of the interferograms to demodulate them.However it has its maximum signal to noise ratio (Eq.( 27)) for a signal tuned at π/2 radians per temporal interferogram.The impulse response h(t) of the 9-step quadrature filter is calculated by the inverse Fourier transform of Eq. ( 42).Once the impulse response is obtained, its real part hr(t) and its imaginary part hi(t) are found, and finally we get our searched 9-step PSI algorithm As mentioned in section 2 if we want to reject higher n-order harmonic (n>1), one needs to generate zeroes at H(-nω 0 ) = H(nω 0 ) = 0, and probably increase the number of building blocks of the PSI filter.

High Q resonant PSI algorithms
We now turn to our last theme on this paper which is the use of recursive linear quadrature filters in PSI.Let us start by showing (in our view) the simplest of such recursive systems, 0 ( , , ) ( , , 1) exp[ ] ( , , ) .

Ic x y t
Ic x y t i I x y t Where we have used unit sampling rate T = 1, and as later we will show, the parameter η<1.0 controls the bandwidth of this quadrature filter.The function Ic(x,y,t) is the output analytical signal associated to our real-valued I(x,y,t) interferogram.The demodulated (estimated) phase of the temporal interferogram in Eq. ( 1) tuned at ω 0 is, Where Re{} and Im{} take the real and imaginary parts of Ic(x,y,t).A block diagram (without the spatial coordinates) of this system is given in Fig. 8.As we can see the impulse response of this filter has infinite time support.However, choosing η = 0.8, only few interferogram samples (about 10) are required to reach the stable regime of the system.On the other side of the fence, choosing η = 0.99 one would need about 200 interferograms to reduce the contribution of the first arrived one by a factor of (0.99) 200 = 0.134, this filter would have a pass-band roughly equivalent to a 200 step PSI algorithm.Intuitively this is why the parameter η controls the band-pass of this resonant filter.The advantage of having a very narrow band-pass response is its extremely high noise immunity; the disadvantage is that we would need to wait many interferogram samples to collect reliable demodulated phase.
Why does one would claim that the system in Eq. ( 44) is a quadrature filter?, because its Frequency Transfer Function H(x,y,ω) says so.Let us take the Fourier transform with respect to t of this equation, one obtains, 0 ( , , ) The FTF is also defined as the ratio between the Fourier transform of the output signal Ic(x,y,ω) to the input signal I(x,y,ω) [6].Which in our current case this ratio equals,

Ic x y H x y F h x y t I x y i
This spectrum represents a resonator centered at -ω 0 , and having a band-width (or Q factor) controlled by η.This FTF may also be seen in phasor representation as, The magnitude and phase delay of this resonant quadrature filter are plotted in Fig. 9 for η = 0.8, and ω 0 = π/2.48).The magnitude response |H(ω)| is resonant at our temporal carrier ω0 = -π/2.Here we have used η = 0.8 which produces a broad resonant peak.The angle introduced by this recursive filter is however "new" in PSI.This is because standard PSI algorithms have all a constant phase delay.This delay in this case is always less than π/2 radians.This quadrature filter poorly rejects the background at ω = 0, and the unwanted analytical signal at ω = π/2.To obtain better rejection at these two frequencies, one may proceed in two ways.One strategy consists in increasing the η parameter to 0.99 to obtain a very high Q resonator.Another strategy is to add to our basic recursive filter Eq. ( 48) with the minimum required zeroes (see Eq. ( 5)) to obtain a bona-fide PSI algorithm.These zeroes may come from our main work-horse system; the tunable 5-step algorithm in Eq. (40).In this way one obtains a FTF given by, With these two zeroes at ω = 0, and at ω = + ω 0 added, our new composed resonator now qualifies as a valid PSI algorithm.This modified FTF spectrum looks as in Fig. 10.This spectral magnitude is almost as good as the one shown for the 9-step algorithm.However the main difference between these two is the non-constant phase delay introduced by the recursive filter.This phase delay may be seen as a drawback by some people, but to others this may be an opportunity to extend the possibilities of PSI algorithms.Another singularity of all recursive filters is their transient response.That is, we need to process (filter) several interferogram samples before the recursive system reaches its stable regime, and thereafter start taking valid phase demodulated data.Finally the recurrence relation that gives rise to the FTF spectrum of our 5-step recursive PSI algorithm is, { } We have changed the notation to the discrete variable n, to obtain a less cumbersome equation.
In our view, the best application for this recursive PSI algorithm would be continuous temporal-interferogram demodulation.For example, using a highly peaked PSI resonator with η = 0.99, and a sampling rate of 200 interferograms per second one would only need to wait one second in order to collect the demodulated phase with extremely low noise, without the overhead computation time of strong spatial filtering.

Conclusions
We have presented the most comprehensive and coherent theory behind few-steps Phase Shifting Interferometry (PSI) to this date.We have presented our theory based on the Frequency Transfer Function (FTF) of the quadrature filter or H(ω) = F[h(t)].We have also presented several applications and analyzed different aspects of PSI such as; 1) reference time-shift invariant spectrum; 2) closed-form formulae for the detuning error in PSI algorithms; 3) noise estimation of the demodulated phase as a function of the spectrum H(ω) of the PSI algorithm used; 4) a well defined strategy to combine existing PSI algorithms; 5) a systematic synthesis technique to design sophisticated PSI algorithms based on two simple units, namely first and second order quadrature filters; 6) finally we have studied some advantages and drawbacks of using recursive digital filtering in PSI, analyzing a simple resonant PSI algorithm.

Fig. 1 .
Fig. 1.Two different spectral magnitudes for |H(ω)| that comply with Eq. (5).The input to these filters is the interferogram I(t) and the output the searched analytical signal Ic(t).These two FTFs are valid PSI algorithms with temporal carrier ω0.Both have the minimum required zeroes for H(ω), these are H(0) = 0, and H(ω0) = 0.

Fig. 3 .
Fig. 3.The interferogram data has a detuning error equal to ∆ radians/interferogram.The spurious signal that contributes to the erroneous phase demodulation is H(ω0 + ∆)exp[iφ].In a well tuned PSI algorithm (∆ = 0), this signal is zero as shown in Fig. 1 and Fig.2.

Fig. 7 .
Fig. 7. Spectral plot of the FTF in Eq. (42).This is an ultra-wide-band filter that corresponds to a 9-step PSI algorithm.This algorithm would demodulate any 9-step sequence of interferograms no matter what carrier they may have.However to obtain the maximum signal to noise ratio, according to Eq. (27), the best carrier is π/2 radians/frame.The "flat zone" is strictly zero only at the tuning frequencies: ω1 = π/4, ω2 = = π/2 and ω3 = 3π/4.But compared to the amount of measuring noise, the non-zero ripples in the right-hand "flat-zone" are small enough to be negligible..

Fig. 8 .
Fig. 8. Block diagram of the resonant PSI algorithm.This is a very simple quadrature filter which requires only a single memory frame to keep the last output complex signal Ic(t).If the η parameter equals 0.99, one would obtain a high Q resonator with an extremely narrow bandwidth roughly equivalent to a 200-step PSI algorithm.The time response of this recursive filter is,

Fig. 9 .
Fig. 9.This is the magnitude |H(ω)| and phase delay introduced by the recursive quadrature filter shown in Eq. (48).The magnitude response |H(ω)| is resonant at our temporal carrier ω0 = -π/2.Here we have used η = 0.8 which produces a broad resonant peak.The angle introduced by this recursive filter is however "new" in PSI.This is because standard PSI algorithms have all a constant phase delay.This delay in this case is always less than π/2 radians.

Fig. 10 .
Fig.10.This is the spectrum of the 5-step recursive resonator in Eq. (50).This PSI filter has the same pole at ω = −π/2 as the basic resonator (with η = 0.8), but now we have added two zeroes at ω = 0, and at the unwanted conjugate signal at ω = + π/2.We have now obtained an almost flat-zero right hand side frequency response. .