1. Introduction

One can usually evaluate the phase map of interferograms by means of two-beam phase-stepping algorithms (TBPSA) [1–4]. These methods can be fully automated and allow one to achieve great accuracy if a sinusoidal dependence of the interferogram intensity with the optical phase is ensured. However, some interferometers of great practical interest do not produce sinusoidal profiles [5, 6]. In the case of a Fizeau interferometer, because of multiple reflections between reference surface and test surface, fringes with an Airy profile are obtained [7]. Thus an error in the phase arises if two-beam phase-stepping algorithms are employed [8], especially for measuring a highly reflective test surface.

To deal with this problem, Clapham and Dew [9] use a lossy reflecting coating on the reference surface to reduce multiple reflections. However, the coating can modify the figure of the reference surface in two aspects: one is deformation of the substrate by stresses in the coating and another one is variations in the phase changes on reflection and transmission due to local variations in the thicknesses of the layers. Hariharan [10] presents a method by introducing an absorbing filter in the test path between the reference surface and the test surface, but the method needs two sets of measurements to reduce the error caused by absorbing filter. In addition, although these two methods can reduce multiple reflections, it cannot eliminate them completely. On the other hand, Hariharan [5] points out that if a fourframe phase calculation algorithm is used, the phase error caused by multiple reflections of a mirror is eliminated to the first-order approximation. Bonsch and Bohme [7] present an algorithm can completely eliminate the phase error (mainly quadruple frequency) introduced by the multiple reflections of a mirror. But both of the two algorithms require accurate phase shifts.

However, it is hard to obtain accurate calibration of the PZT because of the aging of the PZT, hysteresis and environmental changes [11]. Hibino et al [12] and Surrel [13] offer the possibility to reduce the calibration error in the presence of high-order harmonics. It shows that jth-order harmonic can be minimized with phase steps of 2π (j+2) between acquired data frames. Picart et al [14], Patil et al [15] and Langoju et al [16–17] proposed methods for estimating phase information from interferograms with harmonics and unknown phase steps. These methods mainly deal with the phase step, linear and nonlinear miscalibration errors in the PZT response. In practical experiments, however, the phase shifts errors contain not only the phase step, linear and nonlinear miscalibration errors but also the random phase-shift errors caused by environmental vibration. These methods [12–17] seldom consider the problem of random phase-shift errors. On the other hand, many other authors [18–20] have developed iterative techniques to deal with the problem of random phase shifts for two-beam interferograms. But these techniques fail to consider the problem of high-order harmonics. Therefore, if the test surface is filmed with high reflectivity coefficient, such as 0.7 or 0.9, and if the phase shifts is random and unknown, none of the mentioned algorithms [5, 7, 12–20] can deal with them effectively.

In this paper, we present an iterative algorithm to deal with multiple-beam Fizeau interferograms with random phase shifts. We estimate the phase shifts by FFT method and substitute them into the iterative algorithm to reduce the number of iterations. It can extract phase distribution and phase shifts accurately from seven multiple-beam interferograms with random phase shifts. The proposed method has the potential to work efficiently with measuring a highly reflective test surface.

2. Multiple-beam Fizeau interferograms

The Fizeau intensity profile in reflection for two nonabsorbing nearly parallel surfaces is given by the Airy formula that can be written in compact form as [7]

where A is the local mean intensity that depends on the experimental conditions such as intensity of laser source, detector sensitivity etc, and B and C are constants that depend on the reflectivity coefficients of reference surface r ₁ and test surface r ₂

In the Eq. (1), φ denotes the phase function, θ_n means the relative phase shifts for the nth frame and η_n is additive white Gaussian noise with a mean of zero and variance σ². Although Eq. (1) is clearly a nonsinusoidal profile, the term in brackets is a periodic and even function of φ+θ_n. Then it can be conveniently developed in Fourier series form, and Eq. (1) can be rewritten as [8]

It can be seen that the amplitude of the jth harmonic decreases as j becomes higher and also that this behavior is faster when both reflectivity coefficients are lower, since the amplitude is a function of (r ₁ r ₂)^j. Moreover, the rate of amplitude between consecutive harmonics is equal to r ₁ r ₂. According to Eq. (3), the problem of multiple-beam Fizeau interferograms with random phase shifts is, in fact, the problem of interferograms with high-order harmonics and random phase shifts.

3. Phase-shifts estimation by FFT method

Goldberg and Bokor [21] present a method estimating the phase shifts by FFT method without the necessity of iteration if the two-beam interferograms have spatial carrier frequency. Similar to Goldberg and Bokor’s method, we can estimate the phase shifts of multiple-beam Fizeau interferograms by FFT method. Following the derivation of the widely known Fourier-transform method of interferogram analysis [21], in the presence of a spatial carrier frequency k ₀, the piston constant and tilt linear phase terms can be separated from the phase function of interest. We define φ to include the tilt as

where Φ₀(r) denotes the phase to be tested, with no net piston or tilt components. Equation (3) then becomes

The Fourier transform of I_n(r) is G_n(k), which is given by

with A₀(k), C_j(k) and ζ_n(k) being the Fourier transforms of a₀A(r)/2, a_jA(r)exp[ij(Φ₀(r)]/2 and η_n(r), respectively. In the equation,* indicates the complex conjugate. Figures 1(a) and 1(b) show the schematic drawings of G_n(k) for r ₁=0.2, r ₂=0.7 and r ₁=0.7, r ₂=0.9,respectively (here G_n(0) has been weakened). It shows that the high-order harmonics become evident as the product of the reflectivity coefficients (r ₁ r ₂) increases. Typically dominated by low-spatial-frequency components, A₀(k), C_j(k) and ζ_n(k) are all strongly peaked at the zero frequency. So if k ₀ is large enough, the absolute values of $A_0(k),{\displaystyle \sum _{j=2}^{\infty }\exp (ij{\theta }_n)C_j}((1-j)k_0),{\displaystyle \sum _{j=1}^{\infty }\exp (-ij{\theta }_n)C_j^{*}((1+j)k_0)}$ and ζ_n(k ₀) are much smaller than that of C ₁(0), then

where C ₁(0) is a complex constant. Therefore, the phase shifts can be estimated as [20]

where θ ₀=tan^-1[Im(C ₁(0))/Re(C ₁(0))]. θ ₀ denotes a constant phase offset and is the same for each interferogram.

 

Fig. 1. Schematic drawing of frequency spectrum of multiple-beam Fizeau interferogram with spatial carrier frequency k ₀.(a) r ₁=0.2, r ₂=0.7 and (b) r ₁=0.7, r ₂=0.9

Download Full Size | PDF

4. Phase determination by Least-squares iteration

The multiple-beam Fizeau intensity data containing harmonic components up to the p th-order for n=1 to n=N data frames (assuming that Aa_j≪1 for j>P+1), is given by

where, for the nth data frame, b ₀=a ₀ A/2 represents the background intensity and b_j=Aa_j represents the amplitude of the jth-order harmonic components.

Step 1. Determination of phase distribution from known phase shifts.

As in the conventional phase-shifting algorithms [18], defining a new set of variables as X₀=b₀,X_2j-1=b_j cos(jφ),X _2j=-b_j sin(jφ); S₀=1, S _2j-1(n)=cos(jθ_n), S _2j(n)=sin(jθ_n) for j=1, 2, 3 and n=1, 2⋯N; Eq. (10) is rewritten as

where

For the first iteration cycle, θ_n is estimated by FFT method as described in sec.3 and g_n is set to zero. For the other iteration cycles, we can estimate g_n by Eq. (12) according to the iterative results of previous cycle. Therefore, if phase shift θ_n is known and N≥7 is satisfied, the unknown X can be solved by utilizing the overdetermined least-squares method. An expression of the least-squares error can be written as

The minimum of E occurs when the derivative of E with respect to all X_j vanishes. These conditions yield to the following equations:

From Eqs. (14–16), the unknown X can be solved. Then b_j and φ can be determined from

According to Eqs. (4) and (18), we can determine b_j for j>3

We also modify g_n according to given phase shift θ_n and obtained phase φ and b_j for j>3. Consequently, we substitute φ and the modified g_n into step 2 of this iterative cycle.

Step 2. Determination of phase shifts by known phase distribution

Determination of the global phase-step values does not require use of the entire interferogram, thus a subdomain of the data with high signal-to-noise ratio is chosen. In our algorithm, we choose the center part of the interferograms. In this subdomain, it is assumed that the background intensity b ₀ and modulation of harmonics b_j are only functions of frames but constants of pixels. We can define another set of variables as X ^′ ₀=1, X ^′ _2j-1=cos(jφ), X ^′ _2j=sin(jφ); S ^′ ₀(n)=b ₀,S ^′ _2j-1(n)=b_jcos(jθ_n),S ^′ _2j(n)=-b_jsin(jθ_n). If φ is known (as obtained from step 1) and the number of pixels M is larger than 7, the unknowns of S ^′ _2j(n) can be solved again by utilizing the overdetermined least-squares method. An expression for least-squares error can be rewritten as

where m denotes the mth pixel. The minimum of E occurs when the derivative of E with respect to all components of S ^′ vanishes. These conditions yield to the following equations:

From Eqs. (21)–(23), the unknown S ^′(n) can be solved and the phase shift θ_n can be determined from

Here θ_n has been subtracted by θ ₁ because only the relative phase shifts are meaningful. To converge fast, we modify g_n again according to φ and b_j obtained in step 1 and the phase shift θ_n obtained in step 2 of this iterative cycle and then substitute g_n and θ_n into step 1 during the next iterative cycle.

Step 3. Convergence criterion

The algorithm repeats step1 and step 2 until the phase-shift values converge. The convergence criteria for relative phase shifts can be expressed as

where q represents the number of iterations and ε is a predefined accuracy requirement(e.g. ε=10^-3). When the convergence criterion is satisfied, the phase determined from step 1 in the last iterative cycle is unwrapped to obtain the exact phase distribution [22].

5. Numerical simulation and discussion

Numerical simulations are carried out to test the performance of the proposed algorithm. In all the simulations below, we assume that A=exp[-0.02(x ²+y ²)] and Φ₀=π(x ²+y ²), where -1≤x≤1, -1≤y≤1. In addition, we set phase shift θ_n=(n-1)π/4+δ(n) where δ(n) is a random number and set spatial carrier frequency k ₀=√2k_0x=√2k_0y where k_0x and k_0y denote the components of k ₀ along x and y axes. r ₁ is set as 0.2 while r ₂ is in range from 0.1 to 1. When r ₂=0.7, k_0x=k_0y=2π and the white Gaussian noise level is 50 dB, we generate seven multiple-beam Fizeau interferograms with arbitrary given phase shifts. Figure 2(a) shows one of the multiple-beam fringes and Fig. 2(b) shows the phase calculated by our iterative algorithm with 2 cycles. Table 1 presents the phase shifts for FFT-estimation (no iteration) and two iterative cycles. The average errors of phase shifts for them are 0.3304 and 0.0087 radian respectively. Figures 1(c) and 1(d) present the phase errors for no iteration and two iterative cycles, respectively. The root-mean-square (RMS) values of them are 0.0494 and 0.0034 radian respectively. Figure 2 and Table 1 show that the proposed algorithm requires only two iterative cycles and accurately extracts phase from seven multiple-beam Fizeau interferograms with arbitrary phase shifts even if the initial estimated phase shifts have error bigger than 0.3 radian.

Table 1. Phase shifts for FFT-estimation (no iteration) and two iterative cycles. Unit: radian.

View Table

 

Fig. 2. (a). multiple-beam Fizeau interferogram, (b). phase obtained by the proposed iterative algorithm, and (c) and (d) the phase errors for no iteration and two iterative cycles respectively.

Download Full Size | PDF

We find that some factors will affect the performance of the proposed algorithm. Equations (1) and (9) show that the main factors are the spatial carrier frequency, the reflectivity coefficient and the random noise. Here these factors are analyzed and discussed as follows.

5.1 Influence of spatial carrier frequency

Equations (8) and (9) show the phase-shift error for FFT-estimation depends primarily on the spatial carrier frequency k ₀ and the functions of A₀(k), C_j(k) and ζ_n(k). We keep all the parameters of the multiple-beam fringes the same as the above assumptions except the amplitude of spatial carrier frequency. Thus the functions of A₀(k), C_j(k) and ζ_n(k) are fixed. We calculate the average phase-shift errors and the RMS values of phase errors for different spatial carrier frequencies. Figures 3(a) and 3(b) present the average phase-shift error and the RMS value of phase error for no iteration and 2 iterative cycles, respectively. Figure 3(a) shows that the phase-shift error for FFT-estimation scales inversely with the number of fringes in a window. When the number of fringes is more than 8, the average phase-shift error and the RMS value of phase error for no iteration are less than 0.0135 and 0.0034 radian, respectively. In addition, Figs. 3(a) and 3(b) also show that the phase-shift error and the phase error are reduced effectively by only two iterative cycles if the number of fringes is more than 2. For the number of fringes less than 2, the big estimated phase-shift error will lead to big error in obtained phase. Consequently, the estimated g_n by Eq. (12) may produce considerable errors which may result in instability in the iteration. This problem can be solved if the estimated phase-shift error is small enough. Thus we set k_0x=k_0y=4π (that is, the number of fringes is 8) in the simulations below to avoid instability in the iteration.

 

Fig. 3. (a). phase-shift error and (b) phase error for no iteration and two iterative cycles at different spatial carrier frequencies.

Download Full Size | PDF

5.2 Influence of reflectivity coefficient

Here we set r ₁ as 0.2 and let r ₂ vary from 0.1 to 1. The reflectivity coefficient of test surface will affect the amplitude of functions A₀(k) and C_j(k), and thus it will affect the estimated phase-shift error. Figures 4(a) and 4(b) show the average phase-shift errors and the RMS values of phase errors for different r ₂, respectively. It shows that the average phase-shift error and the RMS value of phase error increase with r ₂. The errors are reduced partially after two iterative cycles. For large product of the reflectivity coefficients (r ₁ r ₂), however, the effect of high-order harmonics will be more evident. Therefore, more iterative cycles are needed to compensate the effect of high-order harmonics. For example, when r ₁=0.2 and r ₂=0.9, the average phase-shift error and the RMS values of phase error are 2.2×10^-3 and 7×10^-4 radian for 2 iterative cycles, but they are reduced to 2×10^-4 and 4×10^-4 radian for 4 iterative cycles.

 

Fig. 4. (a). phase-shift error and (b) phase error for no iteration and two iterative cycles at different reflectivity coefficients of test surface.

Download Full Size | PDF

5.3 Influence of random noise

The measurement is invariably sensitive to noise, thus it is important to study the robustness of our algorithm in the presence of noise. We performed simulations at different signal-to-noise ratios (SNR) to compute the average phase-shift errors and the RMS values of phase errors, and the results of them are shown in Fig. 5. Figure 5(a) shows that the average phase-shift error for FFT-estimation and two iterative cycles keep nearly constants when the SNR is larger than 20db. For the SNR near 10db, however, ζ_n(k ₀) is comparable with C ₁(0) in Eq. (12), thus the estimated phase shifts have big errors. Figure 5(b) shows that the RMS values of phase errors decrease with the increasing of SNR. When the SNR becomes larger than 50, the effect of noise on phase error can be ignored (comparing with other factors such as phase-shifts miscalibration, high-order harmonics and so on). Thus the phase errors keep nearly constants when the SNR is larger than 50db. We can conclude from Fig. 5 that the phase-shift error and the phase error can be reduced after two iterative cycles if the SNR is higher than 20db.

In the above simulation, one important random noise, the quantization error of CCD [23], has not been considered. However, the effect of quantization error is more serious for multiple-beam fringe than that for two-beam fringes. Thus it is necessary to analyze the influence of quantization error on the proposed algorithm. The fringe-intensity error caused by quantization [23] is ΔI_n=I_n-INT(I_n), where INT indicates the nearest integer representation. We assume that the maximum of INT(I_n) to be 2^t-1, where t means the quantization bit of CCD. Figures 6(a) and 6(b) show the average phase-shift error and the RMS values of phase errors for different bits of CCD when r2=0.7, respectively. It shows that the phase errors decrease as the bit of CCD increases. As the product of the reflectivity coefficients (r ₁ r ₂) increases, the contrast of the multiple-beam fringe will decrease and the influence of the quantization error will be more serious. Thus CCD with higher quantization bit is required to test optics with high reflectivity coefficient.

 

Fig. 5. (a). phase-shift error and (b) phase error for no iteration and two iterative cycles at different noise levels.

Download Full Size | PDF

 

Fig. 6. (a). phase-shift error and (b) phase error for no iteration and two iterative cycles at different quantization bits.

Download Full Size | PDF

6. Experiments

To illustrate the application of the proposed algorithm, we apply it to the practical interferograms with random phase shifts and high-order harmonics. In the experiment, the reflectivity coefficients of the test and reference surfaces in Fizeau interferometer are 0.7 and 0.2 respectively. We capture 7 frames of multiple-beam interferograms generated in dynamic environments. Figure 7(a) shows one of the interferograms with the contrast about 0.25 and Fig. 7(b) is the intensity gray at the middle row of Fig. 7(a). Because the spatial carrier frequency is not large enough, the estimated phase shifts have big error. After 8 iterative cycles, our algorithm meets the predefined accuracy ε=10^-3. According to the results of our algorithm: the actual relative phase shifts among the 7 frames are 1.1755, 1.9347, 2.9187, 3.8574, 4.7360 and 5.2126 radian; and the obtained phase is shown in Fig. 7(c), with the peak-to-vale (PV) value of 2.22 and RMS value of 0.38 radian. We also measure the test surface by Zygo’s interferometer with a vibration-isolating platform, calibrated PZT, and a proper optical attenuator in Fizeau cavity. Then traditional phase shifting algorithm [13] is utilized to extract the phase shown in Fig. 7(d), with PV value of 2.11 and RMS value of 0.36 radian. Evidently, it can be concluded from Figs. 7(c) and 7(d) that our method eliminates the effects of high-order harmonics and random phase shifts. In addition, Fig. 7(b) shows that the maximum intensity is about 200. If we change the intensity of laser source to make the maximum intensity gray level close to 255(the quantization bit of CCD is 8), the quantization error will be reduced and the accuracy of the obtained phase will be improved.

 

Fig. 7. Experimental results. (a) multiple-beam Fizeau interferogram, (b) the intensity gray at the middle row of Fig. 7(a), (c) phase by the proposed algorithm and (d) phase by Zygo’s standard interferometer with traditional algorithm.

Download Full Size | PDF

7. Conclusion

To conclude, we have proposed a novel method for phase extraction from multiple-beam Fizeau interferograms with random phase shifts. The proposed method requires only two iterative cycles and seven multiple-beam fringes with spatial carrier frequency. It can accurately extract phase information by dealing with the problem of random phase shifts and high-order harmonics which exists in practical multiple-beam Fizeau interferograms. The numerical simulations and experiments show the feasibility of our method. It also shows that (1) the average phase-shift error for FFT-estimation decreases with the number of fringes increasing and with the product of the reflectivity coefficients (r ₁ r ₂) decreasing; (2) the RMS values of phase errors decreases with the number of fringes increasing, with the product of the reflectivity coefficients (r ₁ r ₂) decreasing, and with the increasing of SNR; and (3) the average phase-shift error and the RMS value of phase error are reduced effectively after two iterative cycles. The method has applications in high precision interferometry, especially for measuring a highly reflective test surface.