1. Introduction

One of the major topics of research in interferometry has been the design of efficient algorithms for accurate measurement of wavefront phase. Two major approaches that have been followed in phase measurement can broadly be classified as temporal and spatial. In temporal techniques, three or more data frames are acquired while phase shifts are introduced in successive data frames [1–15]. The phase is subsequently computed at each pixel by using the intensity readings to solve for the phase. One of the devices commonly used for shifting the phase of the interference beams in interferometry is the piezo-actuator transducer (PZT). Unfortunately, the PZT itself is a source of error and one of the problems encountered during its use is the device miscalibration [5,8,12]. Other errors occurring due to the use of a PZT can be attributed to hysteresis, ageing, PZT tilt error, etc. Error in phase measurement can also crop up due to multiple harmonics; a consequence of detector nonlinearity, multiple reflections inside the laser cavity or the phase shifter itself. Few algorithms have been proposed which advocate selection of fixed values of phase steps in order to minimize these error sources [8,12]. These algorithms however restrict the user from applying arbitrary phase steps.

The use of spatial techniques such as Fourier and wavelet transforms has also become popular since the phase information can be retrieved from single data frame [16–18]. The Fourier transform technique has been widely applied in signal and image processing functions by transforming the intensity data into the Fourier domain. To apply the Fourier transform to interferometric signals, a carrier fringe is introduced to separate the first-order signal from zeroth-order components in the Fourier domain. The phase information is then computed by processing the first-order signal. However, the Fourier transform technique is basically a global frequency transformation technique and as such the signal in one position affects the signal in other positions in spectral analysis. To address this concern, wavelet transform was introduced which also allows for direct phase demodulation by doing away with the need of phase unwrapping [19]. However, this technique has met with limited success because of the need to add carrier fringes in both the horizontal and vertical directions. If the fringes are added only in the horizontal direction, only one dimensional wavelet transform is possible. Moreover, the need to introduce an optimization algorithm to minimize the wavelet ridge is one of the potential bottlenecks in accessing its optimal performance.

Since temporal techniques are known to offer advantage over the spatial techniques in terms of improved noise immunity, high-spatial-frequency resolution, and insensitivity to spatial variations in the detector response, these techniques are preferred in phase analysis. Given the advantages possessed by the temporal techniques, few algorithms have also been proposed in which the Fourier transform is applied temporally to estimate the global phase steps. Once the phase steps are estimated between the successive data frames, the least squares fit technique is applied for the estimation of phase. For instance, in the method proposed by Goldberg et al [17] for the estimation of global phase steps, the phase of the first-order maximum of one of the interferograms in the Fourier domain is compared with the first-order maximum of other interferograms. This method however, is highly dependent on the distribution of the zeroth- and the first-order spectrum. In such a case, the selection of higher carrier frequency to separate spectral components would seem to be a solution but its feasibility is limited by the spectral resolution of the detectors. The method is also sensitive to noise occurring during the acquisition of successive data frames. To overcome this problem, Guo et al [18] suggested an algorithm based on energy-minimum to estimate the global phase steps between frames. Although this algorithm is less sensitive to noise, it is iterative and requires an initial guess for minimizing the energy. Moreover, both these concepts assume that the phase step is uniform over the entire data frame. These methods also require the addition of carrier fringes.

Hence, there is a need to design a method which does not require the addition of carrier fringes, handles PZT miscalibration, computes phase step values pixel-wise, and compensates for non-sinusoidal fringe profile. For this, we propose a super-resolution Fourier transform method [20]. Conceptually, if the Fourier transform is applied temporally at (x, y) pixel for N data frames, the phase step information can be retrieved, however, the phase step resolution is limited to 2π/N radians. On the other hand, the proposed algorithm incorporates salient features of super-resolution technique and Fourier transform method. Super-resolution techniques based on the formation of covariance matrix from temporal data set offers the possibility to measure frequencies which are separated by less than 2π/N. In practical terms we first design a covariance matrix from the intensity fringes recorded at the (x, y) pixel for N data frames. The eigenvalue decomposition of this covariance matrix yields signal and noise subspaces. The Fourier transform is next applied to the noise subspace so that the phase step values can be estimated. Applying Fourier transform temporally overcomes the limitation mentioned earlier, compared to when applied spatially. The efficiency of the proposed method is studied in the presence of noise.

The following section introduces super-resolution Fourier transform method. Section 3 presents analysis of the proposed method and comparison with Discrete Fourier transform (DFT). Section 4 presents evaluation of the algorithm followed by phase estimation in Section 5.

2. Phase shifting interferometry

The method consists of acquiring data while voltages are applied to the PZT. The recorded fringe intensity at a point (x, y) of the t^th frame in presence of white Gaussian noise is given by

where, I_dc is the local average value of intensity, a_k is the coefficient, j = √-1, u_k = exp (jkα), superscript * denotes the complex conjugate, η is the white Gaussian noise with mean zero and variance σ ²; and φ, α, and κ represent phase distribution, phase step, and the number of harmonics, respectively. The first step is the computation of covariance for I(t) which is given by [21]

where, E[.] is the expectation operator and ω_k = 0,±α,..,±κα are the harmonics present in the intensity modulation, p is the time offset and term $A_k^2$ has been explained in Appendix. The covariance matrix can thus be written as

where, I(t) = [I(t-1),……,I(t-m)]^T, m is the covariance length, T is the transpose operator, and r* (-p) = r(p). This covariance matrix R _I can be written in a compact form as follows

where, H is the conjugate transpose operator and S = [s _κ ,s _-κ,…,s ₀]_m×(2κ±1). Here, s _k is the Vandermonde matrix of the form s _k =[1,exp (jω_k ),…..,exp(j(m-1)ω_k )]^T; R _s and R _ε are the signal and noise contributions; I _d is the m×m identity matrix; and P _{(2κ+1)×(2κ+1)} is a power matrix which is defined as

Thus, the eigenvalue decomposition of the matrix R _I can be written as

The eigenvalue decomposition of the matrix R _I yields a full set of orthonormal eigenvectors υ, ν and corresponding set of eigenvalues λ ₀ ≥ λ ₁ ≥,….,λ _m-1. The eigenvectors υ = υ ₀,υ ₁,…,υ _2κ span the signal subspace and the corresponding eigenvalues are such that λ ₀ ≥ λ ₁ ≥ …. ≥ λ _2κ > σ ². The remaining eigenvectors ν = ν ₀, ν ₁,…, ν _m-2κ-2 span the noise subspace and their corresponding eigenvalues are λ _2κ+1 =,……..,= λ _m-1≈σ ². Since the signal vectors are orthogonal to the noise subspace, the inner product of each of these subspaces is zero. In addition, since υ is the signal eigenvector, we can write

where, matrix γ is

The eigenvectors υ in Eq. (7) can be expressed as

where C = PSυγ ^-1 and |C| ≠ 0 (since rank (υ) = rank (S) = 2κ+1). Since υ ^H ν = 0, we have

From Eq. (10), we can infer that{s _k ${\}}_{k=0}^{2\kappa }$ ⊥ R(ν), where R(ν) is the span of ν. Thus {ω_k ${\}}_{k=0}^{2\kappa }$ are the unique solutions of S ^H νν ^H S = 0.

The term ${\mathbf{s}}_k^{\mathrm{H}}$ ν _i in Eq. (10) recalls the definition of the inner product of two complex sequences 〈s _k , ν _i 〉 which is given by

Further, from Eq. (10) we know that the signal vectors are orthogonal to the noise vectors which yields the dot product zero (since s _k and ν _i are orthogonal). Hence, combining Eqs. (10) and (11), we obtain

Since s _k = exp(jω_k ) and from the definition of the inner product in Eq. (12), we get

Equation (13) is nothing but the N point Discrete Fourier Transform of the noise vector ν _i evaluated at ω = ω_k . Hence frequencies ω_k are the zeros of 〈s,ν _i 〉. This observation motivates us to form a power spectrum function P_i [exp(jω_f )] which is defined as

where, s _f = [1,exp(jω_f ),…..,exp(jmω_f )]^T and NOF is the number of frequencies. In Eq. (14), we generate ω_f from [0, 2π) in steps of 2π/NOF, such that, at ω_f =ω_k , the inner product in Eq. (14) is zero, which eventually manifests as a peak of infinite amplitude in the power spectrum. Since the eigenvectors corresponding to noise subspace are not of the same power (eigenvalues), usually, a weighted averaging is performed for m-(2κ+1) noise eigenvectors. Hence Eq. (14) for the power spectrum can be written as [20]

Examination of the prominent peaks in P[exp(jω)] gives the estimate of the frequencies present in the signal which in turn gives the information relative to the phase step value applied to the PZT. It is important to emphasize that since the signal subspace is orthogonal to the noise subspace, we generate a set of signal vectors s _f for frequencies between [0, 2π) and perform the inner product with all the noise eigenvectors ν _i such that the resultant is zero. The vector ν _i is obtained from the eigenvalue decomposition of R _I matrix in Eq. (3). Usually, for a large data sample the noise subspace is larger than the signal subspace and hence more information can be extracted from the noise subspace.

A point to note here is that in the Fourier transform method, the frequency resolution is restricted to 2π/N for N data frames. However, the resolution obtained while estimating the power spectrum using Eq. (15) is limited only by the number of frequencies NOF generated between the interval [0, 2π). Hence, in the proposed method the frequency resolution is 2π/NOF for the same number of data frames N. Since, in the proposed method the frequency resolution is no more dependent on the number of data frames acquired but on the size of the NOF, the method is referred to as a super-resolution Fourier transform method. Details on the influence of NOF in estimating the phase step value is presented in Section 4. Next section presents a brief comparison of the proposed concept with the well known discrete Fourier transform (DFT) method.

3. Comparison with DFT methods

The time series harmonics analysis is a crucial problem in signal processing. DFT is one of the methods quite often used to determine the harmonics present in the signal. The DFT for the interferometric data in Eq. (1) for N data frames is given by

The position of the prominent peaks in the amplitude plot of the DFT indicates the frequencies that are present in the signal. These frequencies are nothing but the phase steps applied in the subsequent data frames. If the value of N is small and the separation between the successive phase step values is less than the DFT resolution limit, 2π/N, then these phase steps cannot be retrieved. If the value of N is large enough, then we can achieve good accuracy in the estimation of phase step values. Since the measurement is sensitive to random and systematic error sources, acquiring large number of data frames is not a preferential choice in phase shifting interferometry. One way by which the resolution in DFT can be increased is by padding zeros to the intensity signal in Eq. (16). The DFT of the zero padded signal Î _ZP of length ZP is given by

In Eq. (17), I(t) is zero for N≤t≤ZP-1. The performance of our proposed method over DFT methods is tested by applying it to computer generated fringes. The phase step applied to the PZT is selected as α =0.5760 radian (33°) and two cases are considered: κ=1, N = 8; and κ=2, N = 12.

 

Fig. 1. Plots of power spectrum obtained for κ = 1 using a) Fourier transform b) Zero padded Fourier transform and c) super-resolution Fourier transform for determining the phase step values α.

Download Full Size | PDF

Figures 1(a)–(b) show power spectrum obtained at an arbitrary pixel (x, y) using the DFT and zero padded DFT methods, respectively. For zero padded DFT the intensity signal is padded with 4088 zeros such that the length of the zero padded signal ZP = 4096. It can be observed that determining the phase step values for small data frames is not possible from these two plots. Figure 1(c) shows the power spectrum obtained from the proposed method using Eq. (15). We observe that the peaks are defined clearly and with improved accuracy as compared to the plots obtained in Figs 1(a)–(b). The peaks are observed in this case at α = 0,±0.5768 radian with an error of 0.1389%. The proposed method can thus find the frequency accurately even if the number of frames is limited.

Similarly, Figs. 2(a)–(b) show the power spectrum obtained from DFT and zero padded DFT methods, respectively, for the second case(κ=2). To observe five peaks from small number of data frames is a difficult task. The plot in Fig. 2(c) shows the power spectrum corresponding to frequencies 0,±α, and ±2α obtained using Eq. (15) for NOF = 4096, N = 12, and SNR = 60 dB. In this case the peaks are observed at 0,±0.5768, and ±1.1508 radians with the percentage error of 0.0%, 0.1389%, and 0.0955% respectively

 

Fig. 2. Plots of power spectrum obtained for κ = 2 using a) Fourier transform b) Zero padded Fourier transform and c) super-resolution Fourier transform for determining the phase step values α.

Download Full Size | PDF

Hence, we can state that the super-resolution Fourier transform method is more efficient than both the Fourier transform and zero padded Fourier transform methods.

4. Evaluation of the algorithm

The proposed concept is tested by simulating fringe patterns in Eq. (1) with phase step α = 0.5760 radian (33°), for κ = 1 and κ = 2, and for different number of data frames in the presence of noise. The minimum number of data frames required for retrieving the phase steps are twice the number of frequencies present in the signal. However, the presence of random noise necessitates acquiring large number of data frames. The number frequencies present in the spectrum can be evaluated by well known classical signal processing techniques [22]. The first step in the proposed concept involves the design of a covariance matrix R _I using Eq. (3). The second step involves performing eigenvalue decomposition of matrix R _I which yields signal and noise subspaces as shown in Eq. (6). Further a power spectrum estimation is performed using Eq. (15) for the noise eigenvectors by generating NOF frequencies between 0 and 2π. In the present simulation we consider NOF = 4096 . Finally, the frequencies or the phase step values are estimated pixel wise from the peaks obtained using Eq. (15).

During the simulation, the signal-to-noise ratio (SNR) is varied from 10–80 dB in steps of 0.1 dB. Figures 3(a)–(d) show the influence of an increase in the number of data frames, N = 6,8,10, and 12 , while retrieving the phase step values for κ= 1. As the number of data frames increases, the phase step values can be reliably estimated at lower SNR’s. For instance, Fig. 3(c) shows that phase steps can be estimated, within allowable error in the measurement, for SNR = 40 dB onwards as compared to SNR = 70 dB onwards in Fig. 3(b). Similarly, Fig. 3(d) shows that phase steps can be estimated for SNR = 35 dB onwards.

Figures 3(e)–(h) show plots for retrieving phase steps from a signal with κ=2. The number of data frames is obviously increased in the present case. The number of data frames used during the simulations are N = 10,12,14, and 16. Figures 3(e)–(h) show the influence of an increase in the number of data frames on retrieving the phase step values. Figure 3(h) shows that phase steps can be estimated for SNR = 27 dB onwards as compared to SNR = 45 dB onwards in Fig. 3(f).

 

Fig. 3. Plots of phase step α (in degrees) versus SNR obtained using Eq. (15) at an arbitrary pixel location on a data frame for (a) - (d) κ = 1; N = 6,8,10, and 12, respectively, and (e) - (h) κ= 2; N = 10,12,14, and 16, respectively.

Download Full Size | PDF

Table 1. Phase step estimation for various values of NOF and N = 8 and 12. During simulation the phase step α is taken as 33°, κ=1, and SNR = 40 dB.

View Table

The estimation of phase steps can also be improved with an increase in the number of frequencies NOF generated between [0,2π] in Eq. (15). The analysis in Table 1 shows a reduction in the percentage error as the NOF increases. The simulation also shows a reduction in the percentage error for the same NOF when the number of data frames are increased. For N = 12 and NOF = 4096, typical time for the computation of phase step at a pixel is 0.05 seconds on a desktop PC with Pentium IV 2.66 GHz processor.

 

Fig. 4. Plots of power spectrum for N = 12 data frames and κ= 1 for various orders of nonlinearity in the PZT response to the applied voltage.

Download Full Size | PDF

Interestingly, under the controlled experimental conditions, that is, in the case of high SNR data, the plots in Figs. (1) and (2) are qualitatively indicative of the nonlinearity present in the PZT response to the applied voltage during phase stepping. Analysis of the response of the PZT using Eq. (15) shows that when the PZT exhibits a nonlinear characteristic, the power spectrum obtained is not sparse. Figures 4(a)–(d) demonstrate that as the quadratic nonlinearity of the PZT response increases, the power spectrum shows broadening of the peaks. Hence, for this method to work efficiently, the PZT should be used in the linear region. Usually, if small increments in the voltages are applied to the PZT, then nearly uniform phase steps can be obtained. Ironically, this analysis could also be used as a bench mark for identifying the non-linear characteristics of the PZT.

 

Fig. 5 (a) Fringe pattern and (b) retrieved phase φ.

Download Full Size | PDF

 

Fig. 6. Typical phase error in computation of phase distributions φ (in radians), for the phase steps obtained from Fig 3(d) for SNR=30dB.

Download Full Size | PDF

The phase steps applied for N data frames in Figs 4(a)–(d) is assumed to follow the response α_t =α_t + e ₁(tα)²/2π. During the simulation the phase step α is selected as 0.5760 radian. Figure 4(d) shows that for e ₁ =10%, the power spectrum is not sparse as compared to Fig. 4(a) where sharp peaks can be observed for e ₁ = 0%.

5. Estimation of phase distribution

Once the phase step values have been estimated, phase ϕ can be estimated using the following expression

where, 0,α,2α,…,(N-1)α are the phase steps for frames I ₀,I ₁, I ₂,.., I _N-1 respectively. The phase φ is subsequently computed from the argument of ℓ₁. Recall that ℓ _k has been defined in Eq. (1). Figure 5(b) shows the wrapped phase for the typical simulated fringe map in Fig. 5(a). During the simulation the phase φ in Eq. (1) is assumed to be randomly distributed between 0 and 2π. Figure 6 shows typical error occurring during the estimation of phase φ at SNR = 30 dB, and for the phase step value obtained in Fig. 3(d) for κ = 1.

6. Conclusion

To conclude, we have proposed a novel approach to extract phase information from an interferogram in the presence of harmonics. The proposed method offers the flexibility to use arbitrary phase steps and also does not necessitate the addition of carrier fringes while applying the Fourier transform method. The method thus overcomes, the limitations posed by other well known methods which use tools such as Fourier and wavelet transforms for phase estimation. The resolution in estimating the phase can be enhanced by acquiring large number of data frames and also by increasing the NOF points, albeit at the expense of computational cost. The proposed method is unique in the sense that it not only allows for the selection of arbitrary phase steps and the use of non-collimated waveforms but also allows for the use of multiple PZTs in an optical configuration. Incorporation of two PZTs was first proposed by Rastogi [23] for holographic moiré. Inclusion of two PZTs in an optical set up paves the way for extraction of two displacement components simultaneously. Hence, the proposed method has the potential to cope with configurations containing multiple PZTs.

The proposed method requires a priori assumptions about the frequency contents in the intensity and also the acquired samples should have a relatively sparse spectrum. This, of course, does not matter in this application, since we know the spectrum is rather sparse (in the absence of vibration) and that is the strength of the super-resolution method in this particular application.