Design of algorithms for phase shifting interferometry using self-convolution of the rectangle window

Zhenguang Shi; Jian Zhang; Yongxin Sui; Ji Peng; Feng Yan; Huaijiang Yang

doi:10.1364/OE.19.014671

1. Introduction

An effective means of profiling smooth surfaces interferometrically is to analyze a sequence of fringe patterns shifted in phase relative to each other [1]. The phase difference between an object wavefront and a reference wavefront is precisely controlled by a phase shift mechanism, and the resulting sequential irradiance distribution is exposed by a photoelectric detector. Then the object phase can be obtained from the arctangent of the ratio between two combinations of the observed irradiances, according to the chosen phase shifting algorithm.

The most commonly used phase shift mechanism and photoelectric detector are piezoelectric transducer (PZT) and charge coupled device (CCD). However, due to the imperfect performance of hardware, systematic phase-shift errors caused by PZT nonlinearity and high-order harmonics caused by CCD nonlinearity usually result in systematic errors in retrieving the object information. Furthermore, metrology environments including vibration, airflow and temperature stability et al, would also lead to serious random errors. Therefore the phase shifting algorithm should be designed to suppress the effects of these error sources.

Up to now, many studies have been reported on error-resistant algorithms by which systematic errors can be effectively minimized. The Fourier theory was introduced as a rule to evaluate the quality of algorithms and to derive new ones less sensitive to errors [2–5]. de Groot designed algorithms via integer approximation of well-known window functions such as the von Hanning window [6,7]. Hibino constructed algorithms to eliminate phase errors due to harmonics and phase-shift miscalibration error in terms of solving linear simultaneous equations [8]. Surrel used the characteristic polynomial theory to devise algorithms based on the one-to-one corresponding relation between phase shifting algorithms and characteristic polynomials [9]. Phillion used recursion rules to generate new algorithms that inherit good properties from the old ones [10]. Servin also used the Freischlad spectral analysis theory to combine new algorithms from the basic building blocks [11-12]. However, to our knowledge, the previous methods are either based on choosing the well-known window functions or based on somewhat complicated calculations.

In this paper, a theoretically simple and effective approach is presented to derive new algorithms based on the self-convolution of a rectangle window, which allows generating customized algorithms with very simple calculations just like convolution. In addition, customization can achieve the aim that the constructed algorithms are not only error-resistant to the harmonics of interference signal with any order, but also error-resistant to the phase-shift nonlinearity up to arbitrary order, which could evidently reduce the required performance of PZT and CCD. The error-resistant algorithms can be customized by only two procedures that are hardware calibration and data fitting.

The paper is organized as follows: Section 2 introduces the basic theory about the effects of error sources in reconstructing the object phase and the principle to suppress systematic errors; Section 3 presents the core approach on how to construct the desired window function based on self-convolution of the rectangle window; Section 4 evaluates the self-designed algorithm via theoretical analysis to prove its error-resistant ability and characterizes it with several parameters. Finally Section 5 offers some numerical simulations for checking the performance of the self-designed algorithm compared with the Zygo 13-frame algorithm [13].

2. Principle

The basic idea is to picture the phase shifting algorithms as a filtering process in the frequency domain and to analyze the performance of various algorithms according to their frequency response, following the Freischlad Fourier transform (FT) spectral analysis theory [2].

Most phase shifting algorithms can be characterized as follows,

θ^{'} = \tan^{- 1} (\sum_{n} s_{n} I_{n} / \sum_{n} c_{n} I_{n}) = \arg (\sum_{n} (c_{n} + i s_{n}) I_{n}) = \arg (\sum_{n} w_{n} I_{n}) .

Similarly, they can be expressed in a continuous form,

θ^{'} = \arg (\int w (t) I (t) d t),

where s_n and c_n are the coefficients of algorithms, w_n is the complex coefficients, I_n is the nth fringe,

w (t)

is the continuous window function and

I (t)

is the continuous fringe signal.

Using Parseval’s identity, Eq. (2) can be rewritten by,

θ^{'} = \arg (\int \tilde{W} (v) {\tilde{I}}^{*} (v) d v),

where

\tilde{W} (v)

and

\tilde{I} (v)

are the FT of the window function

w (t)

and the signal

I (t)

respectively.

In view of the nonlinearity of PZT ramping and CCD intensity response, the mathematical model of a distorted signal could be described as the following form [3,5,9-10],

I (t) = \sum_{k = 0}^{q} γ_{k} \cdot \cos (θ_{k} + k (v_{0} t + \sum_{j = 1}^{p} ε_{j} t^{j})),

where

γ_{k}

and

θ_{k}

are the amplitudes and phases of the kth-order harmonics;

θ_{0}

equals zero and

θ_{1}

is the object phase; v ₀ is the PZT ramping velocity and

ε_{j}

is the nonlinear coefficients.

Provided that the coefficients $ε_{j}$ are small enough, the signal $I (t)$ can be approximated by a first-order Taylor series expansion as follows,

I (t) \approx \sum_{k = - q}^{q} \frac{γ_{k}}{2} \exp (i θ_{k}) (\exp (i k v_{0} t) (1 + i k \sum_{j = 1}^{p} ε_{j} t^{j})),

where the Euler formula

\exp (i θ) = \cos (θ) + i \sin (θ)

is used and

γ_{- k} = γ_{k}, θ_{- k} = - θ_{k}

.

The Fourier transform of the signal can be written by,

\tilde{I} (v) = \sum_{k = - q}^{q} \frac{γ_{k}}{2} \exp (i θ_{k}) (δ (v - k v_{0}) + k \sum_{j = 1}^{p} i^{j + 1} ε_{j} δ^{(j)} (v - k v_{0})) .

As depicted in Fig. 1 , the harmonics of the signal spectrum are mainly the integer multiple of the PZT linear ramping velocity as shown in Fig. 1.

Fig. 1 The signal spectrum

Download Full Size | PDF

According to the property of delta function,

$\int x (t) δ^{(j)} (t - t_{0}) d t = x^{(j)} (t_{0}), j i s integer,$

Equation (3) can be evolved into the following form,

θ^{'} = \arg (\sum_{k = - q}^{q} \frac{γ_{k}}{2} \exp (i θ_{- k}) (\tilde{W} (k v_{0}) + k \sum_{j = 1}^{p} {(- i)}^{j + 1} ε_{j} {\tilde{W}}^{(j)} (k v_{0}))) .

So the necessary and sufficient condition for $θ^{'} = θ_{1}$ is as follows,

{\tilde{W}}^{(j)} (k v_{0}) = 0, j = 0, 1, \dots, p, k = - q, \dots, - 2, 0, 1, \dots, q .

From Eq. (8), the FT of the window function must have 2q equidistant multiple roots of order p+1 except at k = −1. Namely it can filter out the unwanted components with the first negative fundamental frequency passed. An example is shown in Fig. 2 . The greater the p and q values satisfy Eq. (8), the stronger the ability of the algorithm error-resistant to the PZT and CCD nonlinearity will be, but always at the cost of more needed frames. In what follows, a self-convolution construction method is presented.

Fig. 2 The amplitude spectrum of windows

Download Full Size | PDF

3. Design method based on self-convolution

As is known, the Fourier transform of a rectangle window is a sinc function which has equidistant simple roots [7]. If the rectangle window is convolved with itself once, the Fourier transform of the obtained window function will have multiple roots of order two, increased by one compared to the rectangle window. The theoretical basis is that the Fourier transform of convolution of two functions in the time domain is equal to the products of their Fourier transforms in the frequency domain [7]. When the rectangle window is convolved with itself repeatedly, the window functions which meet the requirements of Eq. (8) might be obtained.

Firstly, construct a time domain rectangle window function as follows,

x (t) = {\begin{cases} 1, | t | \leq T_{0} / 2 \\ 0, | t | > T_{0} / 2 \end{cases}, T_{0} = \frac{2 π}{v_{0}},

whose Fourier transform is a sinc function,

\tilde{X} (v) = \frac{2 π}{v_{0}} \sin c (\frac{v}{v_{0}} π),

here

\tilde{X} (v)

has simple roots at

v = k v_{0}

for all but

k = 0

.

According to the convolution property of Fourier transform, define a function as follows,

y (t) = x (t) \underset{p}{\underset{}{* \dots *}} x (t),

so its Fourier transform is the p + 1 power of a sinc function,

\tilde{Y} (v) = {\tilde{X}}^{p + 1} (v),

where

\tilde{Y} (v)

has multiple roots of order p+1 at

v = k v_{0}

for all but

k = 0

.

Based on the frequency-shifting property of Fourier transform, define a window function,

w (t) = \exp (- i v_{0} t) y (t),

so its Fourier transform is,

{\tilde{W}}_{c} (v) = \tilde{Y} (v + v_{0}),

Obviously ${\tilde{W}}_{c} (v)$ has multiple roots of order p+1 at $v = k v_{0}$ for all but $k = - 1$ .

In practice, a few fringe intensity values are discretely sampled in interferometry. Hence the window function w(t) should be sampled. Suppose the sampling function is, $p_{s} (t) = \sum_{m = - \infty}^{+ \infty} δ (t - m T_{s}),$ so the discrete sampling window is, $w_{d} (t) = w (t) p_{s} (t),$ and the corresponding Fourier transform is,

{\tilde{W}}_{d} (v) = \frac{1}{T_{s}} \sum_{m = - \infty}^{+ \infty} {\tilde{W}}_{c} (v - m v_{s}), v_{s} = \frac{2 π}{T_{s}} .

As seen from Eq. (15), ${\tilde{W}}_{d} (v)$ is a periodic function with period $v_{s}$ , and ${\tilde{W}}_{c} (v)$ has equidistant roots at the integral multiple of the fundamental frequency $v_{0}$ , so in order to make ${\tilde{W}}_{d} (v)$ meet the requirement in Eq. (8), the following condition must be satisfied,

r = \frac{v_{s}}{v_{0}} \geq q + 2 and r i s integer
.

Up to now, the continuous case of designing the window function of phase shifting algorithm is finished. Next it’s the turn to show the discrete realization of the algorithm.

The construction method is very simple: firstly construct a vector of r(= q+2) dimensions; secondly self-convolve the vector p times; at last multiply the self-convolved vector with a phase factor to shift the frequency components of the Fourier domain. Then the coefficients of the algorithm are the real and imaginary parts respectively, and the phase shift step is $Δ φ = 2 π / r$ . The number of fringe patterns needed is (p+1)r-p. Usually r = q+2 is chosen to design the algorithm with the least number of fringes.

For example, if the nonlinear orders of PZT and CCD are 3 and 2 respectively, that is p = 3, q = 2, the discrete approach can be done as follows: A four dimensional vector $a = [\begin{matrix} 1 & 1 & 1 & 1 \end{matrix}]$ is constructed first, then self-convolve it three times, so a vector is derived as, $b = [\begin{matrix} 1 & 4 & 10 & 20 & 31 & 40 & 44 & 40 & 31 & 20 & 10 & 4 & 1 \end{matrix}],$ after shifting the frequency components, the following vector is obtained, $w = [\begin{matrix} 1 & 4 i & - 10 & - 20 i & 31 & 40 i & - 44 & - 40 i & 31 & 20 i & - 10 & - 4 i & 1 \end{matrix}] .$

Then the parameters of our 13-frame phase shifting algorithm are

\begin{array}{l} s = imag (w) = [\begin{matrix} 0 & 4 & 0 & - 20 & 0 & 40 & 0 & - 40 & 0 & 20 & 0 & - 4 & 0 \end{matrix}] \\ c = real (w) = [\begin{matrix} 1 & 0 & - 10 & 0 & 31 & 0 & - 44 & 0 & 31 & 0 & - 10 & 0 & 1 \end{matrix}] \\ Δ φ = \frac{π}{2} \end{array}

4. Theoretical analysis

To verify the approach developed in Section 3, we analyze the algorithm in Eq. (17) in theory. The window function and its Fourier transform curves are shown in Fig. 3 .

Fig. 3 The correlative curve of the designed phase shifting algorithm in Eq. (17). (a) The window function; (b) its amplitude spectrum.

Download Full Size | PDF

In Fig. 3(b), it is obvious that the period of the amplitude spectrum is 4 and there are multiple roots of order 4 at integers except at −1, which means that the Fourier transforms of the associated window function satisfy the following equation,

{\tilde{W}}^{(j)} (k v_{0}) = 0, (j = 0, 1, 2, 3, k = - 2, 0, 1, 2) and \tilde{W} (v) is of period 4 .

So we can affirm that this algorithm should be error-resistant to as high as third-order phase-shift nonlinearity and second-order signal nonlinearity simultaneously.

To provide evidence for the claim, we calibrate two main interferometer hardware devices including PZT and CCD, and then obtain the following calibration curves shown in Fig. 4(a) and Fig. 4(b). The measured values are fitted with polynomials, in which the well fitted orders of PZT and CCD values are three and two respectively.

Fig. 4 The calibration curve of hardware devices. (a) PZT calibration curve; (b) CCD calibration curve. Inside each curve, the red marks mean measured points and the blue line mean fitted polynomial curve.

Download Full Size | PDF

The fitted equation of the PZT curve is,

\begin{matrix} P (x) = (p_{0} + x) + p_{1} x + p_{2} x^{2} + p_{3} x^{3} \\ = (- 15.79 + x) - 0.25 x + 0.00057 x^{2} - 0.00000031 x^{3} \end{matrix}

and the fitted equation of the CCD curve is,

C (x) = a_{0} + a_{1} x + a_{2} x^{2} = 8.93 + 180.30 x + 42.85 x^{2},

where p ₁, p ₂ and p ₃ are the PZT miscalibration error, second-order and third-order nonlinearity coefficients respectively. The coefficients a ₀, a ₁ and a ₂ of CCD curve will affect the amplitude of harmonics of an interference signal. Because the coefficients p ₁, p ₂ and p ₃ are close to zero, taking first-order Taylor series approximation like Eq. (5) is reasonable.

An ideal interference signal is $I_{o} = A_{0} + A_{1} \cos (θ),$ but according to the CCD fitted model, the nonlinearity will distort the signal as follows,

\begin{matrix} I_{r} = C (I_{o}) = a_{0} + a_{1} I_{o} + a_{2} I_{o}^{2} \\ = (a_{0} + a_{1} A_{0} + a_{2} A_{0}^{2} + \frac{a_{2} A_{1}^{2}}{2}) + (a_{1} A_{1} + 2 a_{2} A_{0} A_{1}) \cos (θ) + \frac{a_{2} A_{1}^{2}}{2} \cos (2 θ) \\ = B_{0} + B_{1} \cos (θ) + B_{2} \cos (2 θ) \end{matrix}

So the nonlinearity of CCD will result in harmonics of higher order. The higher the CCD nonlinearity order is, the higher order the distorted signal harmonics will have.

To take into account the nonlinearity of PZT, the sampled frame of the interference signal will be distorted further as follows, $\begin{matrix} I_{n}^{'} = B_{0} + B_{1} \cos (θ + \frac{4 π}{λ} P (n v_{0} \cdot \frac{λ}{4 π})) + B_{2} \cos (2 (θ + \frac{4 π}{λ} P (n v_{0} \cdot \frac{λ}{4 π}))) \\ \approx B_{0} \\ + \frac{B_{1}}{2} \exp (i (θ + n v_{0})) (1 + i p_{1} n v_{0} + i p_{2} {(n v_{0})}^{2} \frac{λ}{4 π} + i p_{3} {(n v_{0})}^{3} {(\frac{λ}{4 π})}^{2}) \\ + \frac{B_{1}}{2} \exp (- i (θ + n v_{0})) (1 - i p_{1} n v_{0} - i p_{2} {(n v_{0})}^{2} \frac{λ}{4 π} - i p_{3} {(n v_{0})}^{3} {(\frac{λ}{4 π})}^{2}) \\ + \frac{B_{2}}{2} \exp (2 i (θ + n v_{0})) (1 + 2 i p_{1} n v_{0} + 2 i p_{2} {(n v_{0})}^{2} \frac{λ}{4 π} + 2 i p_{3} {(n v_{0})}^{3} {(\frac{λ}{4 π})}^{2}) \\ + \frac{B_{2}}{2} \exp (- 2 i (θ + n v_{0})) (1 - 2 i p_{1} n v_{0} - 2 i p_{2} {(n v_{0})}^{2} \frac{λ}{4 π} - 2 i p_{3} {(n v_{0})}^{3} {(\frac{λ}{4 π})}^{2}) \end{matrix}$

\begin{matrix} = B_{0} \cdot I_{n}^{0} \\ + \frac{B_{1}}{2} \exp (i θ) \cdot I_{n}^{1} (1 + i p_{1} v_{0} \cdot n + i p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot n^{2} + i p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot n^{3}) \\ + \frac{B_{1}}{2} \exp (- i θ) \cdot I_{n}^{- 1} (1 - i p_{1} v_{0} \cdot n - i p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot n^{2} - i p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot n^{3}) \\ + \frac{B_{2}}{2} \exp (2 i θ) \cdot I_{n}^{2} (1 + 2 i p_{1} v_{0} \cdot n + 2 i p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot n^{2} + 2 i p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot n^{3}) \\ + \frac{B_{2}}{2} \exp (- 2 i θ) \cdot I_{n}^{- 2} (1 - 2 i p_{1} v_{0} \cdot n - 2 i p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot n^{2} - 2 i p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot n^{3}) \end{matrix}

where

I_{n}^{0} = 1, I_{n}^{1} = \exp (i n v_{0}), I_{n}^{- 1} = \exp (- i n v_{0}), I_{n}^{2} = \exp (2 i n v_{0}), I_{n}^{- 2} = \exp (- 2 i n v_{0}) .

According to the differential property of Fourier transform, $x (n) \leftrightarrow \frac{d X (v)}{d v} \Rightarrow n^{j} x (n) \leftrightarrow i^{j} \frac{d^{j} X (v)}{d v^{j}},$ the Fourier transform of $I_{n}^{'}$ can be expressed as,

\begin{matrix} I_{v}^{'} = B_{0} {\tilde{I}}_{v}^{0} \\ + \frac{B_{1}}{2} \exp (i θ) ({\tilde{I}}_{v}^{1} + i^{2} p_{1} v_{0} \cdot \frac{d {\tilde{I}}_{v}^{1}}{d v} + i^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \frac{d^{2} {\tilde{I}}_{v}^{1}}{d v^{2}} + i^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \frac{d^{3} {\tilde{I}}_{v}^{1}}{d v^{3}}) \\ + \frac{B_{1}}{2} \exp (- i θ) ({\tilde{I}}_{v}^{- 1} - i^{2} p_{1} v_{0} \cdot \frac{d {\tilde{I}}_{v}^{- 1}}{d v} - i^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \frac{d^{2} {\tilde{I}}_{v}^{- 1}}{d v^{2}} - i^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \frac{d^{3} {\tilde{I}}_{v}^{- 1}}{d v^{3}}) \\ + \frac{B_{2}}{2} \exp (2 i θ) ({\tilde{I}}_{v}^{2} + 2 i^{2} p_{1} v_{0} \cdot \frac{d {\tilde{I}}_{v}^{2}}{d v} + 2 i^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \frac{d^{2} {\tilde{I}}_{v}^{2}}{d v^{2}} + 2 i^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \frac{d^{3} {\tilde{I}}_{v}^{2}}{d v^{3}}) \\ + \frac{B_{2}}{2} \exp (- 2 i θ) ({\tilde{I}}_{v}^{2} - 2 i^{2} p_{1} v_{0} \cdot \frac{d {\tilde{I}}_{v}^{2}}{d v} - 2 i^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \frac{d^{2} {\tilde{I}}_{v}^{2}}{d v^{2}} - 2 i^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \frac{d^{3} {\tilde{I}}_{v}^{2}}{d v^{3}}) \end{matrix}

where

\begin{matrix} I_{n}^{0} \leftrightarrow {\tilde{I}}_{v}^{0} = 2 π \sum_{l = - \infty}^{+ \infty} δ (v - 2 π l) & I_{n}^{1} \leftrightarrow {\tilde{I}}_{v}^{1} = 2 π \sum_{l = - \infty}^{+ \infty} δ (v - v_{0} - 2 π l) \\ I_{n}^{- 1} \leftrightarrow {\tilde{I}}_{v}^{- 1} = 2 π \sum_{l = - \infty}^{+ \infty} δ (v + v_{0} - 2 π l) & I_{n}^{2} \leftrightarrow {\tilde{I}}_{v}^{2} = 2 π \sum_{l = - \infty}^{+ \infty} δ (v - 2 v_{0} - 2 π l) \\ I_{n}^{- 2} \leftrightarrow {\tilde{I}}_{v}^{- 2} = 2 π \sum_{l = - \infty}^{+ \infty} δ (v + 2 v_{0} - 2 π l) \end{matrix}

According to Eq. (3), one can get $\begin{matrix} θ^{'} = \arg (\int \tilde{W} (v) {\tilde{I}}^{*} (v) d v) \\ = \arg (\begin{array}{l} B_{0} \cdot 2 π \sum_{l = - \infty}^{+ \infty} \tilde{W} (2 π l) \\ + \frac{B_{1}}{2} \exp (- i θ) \cdot 2 π (\sum_{l = - \infty}^{+ \infty} \tilde{W} (v_{0} + 2 π l) + {(- i)}^{2} p_{1} v_{0} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(1)} (v_{0} + 2 π l) + {(- i)}^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(2)} (v_{0} + 2 π l) + {(- i)}^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(3)} (v_{0} + 2 π l)) \\ + \frac{B_{1}}{2} \exp (i θ) \cdot 2 π (\sum_{l = - \infty}^{+ \infty} \tilde{W} (- v_{0} + 2 π l) + {(- i)}^{2} p_{1} v_{0} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(1)} (- v_{0} + 2 π l) + {(- i)}^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(2)} (- v_{0} + 2 π l) + {(- i)}^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(3)} (- v_{0} + 2 π l)) \\ + \frac{B_{2}}{2} \exp (- 2 i θ) \cdot 2 π (\sum_{l = - \infty}^{+ \infty} \tilde{W} (2 v_{0} + 2 π l) + 2 {(- i)}^{2} p_{1} v_{0} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(1)} (2 v_{0} + 2 π l) + 2 {(- i)}^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(2)} (2 v_{0} + 2 π l) + 2 {(- i)}^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(3)} (2 v_{0} + 2 π l)) \\ + \frac{B_{2}}{2} \exp (2 i θ) \cdot 2 π (\sum_{l = - \infty}^{+ \infty} \tilde{W} (- 2 v_{0} + 2 π l) + 2 {(- i)}^{2} p_{1} v_{0} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(1)} (- 2 v_{0} + 2 π l) + 2 {(- i)}^{3} p_{2} v_{0}^{2} \frac{λ}{4 π} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(2)} (- 2 v_{0} + 2 π l) + 2 {(- i)}^{4} p_{3} v_{0}^{3} {(\frac{λ}{4 π})}^{2} \cdot \sum_{l = - \infty}^{+ \infty} {\tilde{W}}^{(3)} (- 2 v_{0} + 2 π l)) \end{array}) \end{matrix}$

To ensure $θ^{'} = θ$ , the Fourier transform $\tilde{W} (v)$ must satisfy the following equation, ${\tilde{W}}^{(j)} (k v_{0} + 2 π l) = 0, (j = 0, 1, 2, 3, k = - 2, 0, 1, 2, l = - \infty, \dots, + \infty),$ which is consistent with Eq. (18) as long as $2 π / v_{0}$ equals four. So it’s proved that the self-designed 13-frame algorithm is extremely error-resistant to three-order PZT nonlinearity and two-order CCD nonlinearity at the same time.

The suppression ability of a phase shifting algorithm to higher harmonics can be estimated by the highest side-lobe level [7]. The smaller level stands for better suppression ability to harmonics. The immunity of a phase shifting algorithm to phase-shift nonlinearity can be estimated by the order of zeros at the integer harmonics of the window function. The bigger the value of the order of zeros is, the stronger the error-resistant ability to phase-shift nonlinearity will be. The sensitivity of a phase shifting algorithm to random noise can be estimated by the equivalent noise bandwidth (ENBW for short [5,7,10]). Usually the smaller the ENBW value is, the stronger the error-resistant capability to the random noise will be.

Table 1 shows the error-resistant capability of algorithms with several window functions. When the rectangle window is self-convolved repeatedly, the side-lobe level will fall off quickly and the order of error-resistant capability for phase-shift nonlinearity will increase evidently, but at the cost of being more sensitive to random noise.

Table 1. Performance of several algorithms for higher harmonic suppression, compensation for phase-shift error and sensitivity to random noise.

View Table

5. Numerical simulations

Next we will give some numerical simulations to compare the Zygo 13-frame high resolution phase shifting algorithm [4] and the self-designed 13-frame algorithm expressed in Eq. (17).

The amplitude and log-amplitude spectrums of the Zygo and self-designed algorithms are shown in Fig. 5(a) and Fig. 5(b) respectively. Obviously the highest side-lobe level of the self-designed algorithm is smaller, and the spectrum curve of the self-designed algorithm touches zeros more smoothly at integer harmonics, which means that the spectrum of the self-designed algorithm has higher-order derivatives equal to zero than the Zygo algorithm. It indicates more error-resistant ability to harmonics and phase-shift nonlinearity.

Fig. 5 Comparison of Zygo 13-frame and self-designed algorithm. (a) The amplitude spectrum of two algorithms; (b) the log-amplitude spectrum of two algorithms; (c) the PV phase error versus the miscalibration coefficient; (d) the PV phase error versus the second-order nonlinearity coefficient; (e) the PV phase error versus the third-order nonlinearity coefficient; (f) the PV phase error versus the coefficient of second order signal harmonics.

Download Full Size | PDF

The PV phase error versus the PZT miscalibration coefficient, second-order nonlinearity coefficient and third-order nonlinearity coefficient curves of the two algorithms is shown in Fig. 5(c), Fig. 5(d) and Fig. 5(e) respectively. In all curves, the second harmonic component is with amplitude 20% of the fundamental, and the PZT miscalibration, second-order and third-order nonlinearity coefficients ranges from −0.2 to 0.2, from −0.04 to 0.04, from −0.008 to 0.008 respectively. As seen from the curves, the PV phase errors for the self-designed algorithm are much smaller than those for the Zygo algorithm on a large scale. So it proves that the self-designed algorithm is more error-resistant to PZT nonlinearity errors.

The PV phase error versus the coefficient of signal second order harmonics curves of two algorithms is shown in Fig. 5(f), where the PZT miscalibration, second-order nonlinearity and third-order nonlinearity coefficients are −0.2, 0.04 and 0.008 respectively and the second order harmonic coefficient ranges from 0 to 0.4. Obviously the PV phase error for the self-designed algorithm is much smaller than that for the Zygo algorithm on a large scale. Even if the second order harmonic coefficient reaches 0.4 of the fundamental one, the PV phase error of the self-designed algorithm would be just 0.025 radians, far smaller than 0.2 radians from the Zygo algorithm. So it proves that the self-designed algorithm is more error-resistant to high-order harmonics resulted from CCD nonlinearity.

In general, to reach the same accuracy, the self-designed 13-frame algorithm allows a larger range of the miscalibration or nonlinearity coefficients of PZT and CCD compared to the Zygo 13-frame algorithm, which means that the self-designed phase shifting algorithm can tolerate greater systematic errors in the interferometer hardware.

6. Summary

This paper presents a comprehensive theory to design a new effective phase shifting algorithm error-resistant to various error sources, such as the phase-shift nonlinearity, the existing high-order harmonic signal and random noise et al. Through analysis, we know that the merits of the phase shifting algorithms completely depend on the Fourier frequency domain characteristics of the corresponding window function. And we also use the self-convolution of the rectangle window to design the satisfying window function. According to the analysis and simulation section, the self-designed phase shifting algorithm is very straightforward and more error-resistant to error sources. Future research will focus on how to reduce the random errors simultaneously while reducing the systematic errors.

References and links

1. H. Schreiber and J. H. Bruning, “Phase shifting interferometry,” in Optical Shop Testing Third Edition, D. Malacara, (New York, 2007), pp. 547–666.

2. K. Freischlad and C. L. Koliopoulos, “Fourier description of digital phase-measuring interferometry,” J. Opt. Soc. Am. A 7(4), 542–551 (1990). [CrossRef]

3. H. Zhang, M. J. Lalor, and D. R. Burton, “Robust, accurate seven-sample phase-shifting algorithm insensitive to nonlinear phase-shift error and second-harmonic distortion: a comparative study,” Opt. Eng. 38(9), 1524–1533 (1999). [CrossRef]

4. P. de Groot, “Measurement of transparent plates with wavelength-tuned phase-shifting interferometry,” Appl. Opt. 39(16), 2658–2663 (2000). [CrossRef] [PubMed]

5. K. Hibino, R. Hanayama, J. Burke, and B. F. Oreb, “Tunable phase-extraction formulae for simultaneous shape measurement of multiple surfaces with wavelength-shifting interferometry,” Opt. Express 12(23), 5579–5594 (2004). [CrossRef] [PubMed]

6. P. Groot, “Derivation of algorithms for phase-shifting interferometry using the concept of a data-sampling window,” Appl. Opt. 34(22), 4723–4730 (1995). [CrossRef] [PubMed]

7. F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proc. IEEE 66(1), 51–83 (1978). [CrossRef]

8. K. Hibino, B. F. Oreb, D. I. Farrant, and K. G. Larkin, “Phase shifting for non-sinusoidal waveforms with phase shift errors,” J. Opt. Soc. Am. A 12(4), 761–768 (1995). [CrossRef]

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

10. D. W. Phillion, “General methods for generating phase-shifting interferometry algorithms,” Appl. Opt. 36(31), 8098–8115 (1997). [CrossRef] [PubMed]

11. M. Servin, J. C. Estrada, and J. A. Quiroga, “The general theory of phase shifting algorithms,” Opt. Express 17(24), 21867–21881 (2009). [CrossRef] [PubMed]

12. 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). [CrossRef] [PubMed]

13. Zygo is a registered trademark of Zygo Corporation.

Window functions	Highest side-lobe level (dB)	Capability for phase-shift error	Sensitivity to random noise	Number of Frames
Rectangle	−12	None	1.00	4
self-convolved once	−23	of first order	1.33	7
self-convolved twice	−34	of second order	1.65	10
self-convolved thrice	−45	of third order	1.92	13
von Hanning	−32	moderate	1.50	13
Hamming	−38	moderate	1.36	13

Design of algorithms for phase shifting interferometry using self-convolution of the rectangle window

Abstract

1. Introduction

2. Principle

3. Design method based on self-convolution

4. Theoretical analysis

5. Numerical simulations

6. Summary

References and links

Cited By

Figures (5)

Tables (1)

Equations (23)

Optics Express