## Abstract

In this work, a new method for surface extraction in white light scanning interferometry (WLSI) is introduced. The proposed extraction scheme is based on the Teager-Kaiser energy operator and its extended versions. This non-linear class of operators is helpful to extract the local instantaneous envelope and frequency of any narrow band AM-FM signal. Namely, the combination of the envelope and frequency information, allows effective surface extraction by an iterative re-estimation of the phase in association with a new correlation technique, based on a recent TK cross-energy operator. Through the experiments, it is shown that the proposed method produces substantially effective results in term of surface extraction compared to the peak fringe scanning technique, the five step phase shifting algorithm and the continuous wavelet transform based method. In addition, the results obtained show the robustness of the proposed method to noise and to the fluctuations of the carrier frequency.

© 2014 Optical Society of America

## 1. Introduction

White light scanning interferometry (WLSI) is an important technique for carrying out 3D roughness and surface profile measurements of microscopic surfaces. A high precision is required in order to control the fabrication techniques of new materials, microelectronic devices and MicroElectroMechanical Systems (MEMS) [1]. Moreover, the strength of the WLSI technique lies in its ability to provide fast, robust and high resolution measurements. Another main advantage of WLSI is its ability to perform unambiguous phase extraction at step heights greater than *λ*/2 compared with monochromatic i.e, single-wavelength interferometry, which suffers from the problem of phase ambiguity due to the periodicity of the fringes. Most of the WLSI methods are based on the AM-FM signal model which represents the variations in light intensity measured along the optical axis of an interference microscope. The processed signal carries information related to the axial position of the surface simultaneously within the envelope and the phase. An example of such methods is the five sample adaptive technique which detects the peak of the envelope by using only five adjacent samples along the optical axis [2]. Signal processing transformations such as the Fourier transform [3], the Windowed Fourier transform (WFT) [4] or the Continuous wavelet transform (CWT) [5, 6] have also been used as processing tools in WLSI. Based on these transformations, phase retrieving strategies have been developed for Zero optical path difference (ZOPD) compensation by extracting the local maxima (ridge detection) of the CWT (resp. WFT). For example, the CWT and the WFT are used to compute the phase information at a local maximum, and have been applied to fringe pattern analysis [7,8]. A class of non-linear and instantaneous methods based on the Teager-Kaiser (TK) energy operator [9] and its extended versions has shown its efficiency for tracking local features such as the instantaneous amplitude and frequency of AM-FM signals [10–15]. Thus, different techniques based on the TK operators have been developed for envelope extraction of locally 1D and 2D orientated patterns [13, 14]. These methods have recently been extended to multidimensional signals and to higher orders where it has been shown that instantaneous frequency estimation can be performed in an efficient way for any local AM-FM signal [15]. All these TK operators are limited to the tracking of the energy of a single signal.

Recently, a cross-energy operator based on the TK operator, called Ψ_{B}, has been introduced to track the interaction energy between two signals [16]. This non-linear energy tracking operator deals with complex signals and its usefulness for non-stationary signals analysis has been demonstrated [17,18]. More precisely, Ψ_{B} can be used as a correlator in both the time [19] and frequency [18] domains. Up until now applications of TK operators in WLSI have been limited to envelope estimation [20]. To our knowledge, there has been no attempt in WLSI to retrieve phase information from fringe patterns exploiting the potentialities of the TK energy operators i.e, instantaneous frequency and envelope estimation, and the Ψ_{B} correlator as local and non-linear methods. In this work, we show how this class of operators can be used in WLSI, computing the local phase by correlation, according to the estimated envelope and frequency to improve the surface measurement.

## 2. Experimental system and interferometric signal

WLSI makes use of a series of white light fringes superimposed on an image of the sample on the camera target, the central fringe corresponding to the position of the surface along the optical axis. The fringes are scanned over the whole depth of the sample surface by modifying the distance between the objective and the sample. A series of images is acquired with a camera at regular intervals to give a stack of *xyz* images. The fringe signals along the optical *z*-axis at each pixel in an image are processed to determine the fringe envelope, with the peaks indicating the positions of the surface. For polychromatic interferometry, the total intensity is the sum of the interferences at each wavelength. A typical intensity signal obtained from a digital camera as the OPD is varied in the interferometer at a given point (*x*, *y*) on the material surface, can be approximated along the optical *z*-axis by a modulated sinusoid as follows [2, 21]:

*z*is a vertical scanning position along the optical axis,

*h*(

*x*,

*y*) represents the height of the surface,

*a*(

*x*,

*y*,

*z*) is an offset intensity containing low frequency components,

*b*(

*x*,

*y*) is a factor proportional to the reflected beam intensity, and

*α*(

*x*,

*y*) is an additional phase offset and

*C*(

*x*,

*y*,

*z*) is the envelope. The parameter

*lc*represents the coherence length of the light source and

*λ*

_{0}is the average wavelength of the light source. Generally the phase offset varies slowly from one point (

*x*,

*y*) to the next, and can be neglected, since only relative heights of the surface matter.

The main challenge consists in determining the height at each point of the surface by exploiting the information provided by both the envelope and the phase simultaneously. Since *λ*_{0} is an averaged value, its possible fluctuations have been taken account using a flexible approach. Most often, the analysis of such signals is limited to the optical *z*-axis. To overcome this problem and for more robust analysis, the signal could be processed according to slice (*x*, *z*) to take into account the neighborhood information.

## 3. Principles of the Teager-Kaiser energy operators

The TK energy operators [9, 10] are efficient non-linear and local algorithms for envelope detection and phase retrieval from AM-FM signals such as those given by Eq. (1). The output of the continuous 1D TK energy operator, noted Ψ, for a real-valued signal *x*(*t*) is given by

*ẋ*(

*t*) and

*ẍ*(

*t*) denote the first and the second time derivative of

*x*(

*t*) respectively. Under realistic conditions [11], when applied to AM-FM signal

*s*(

*t*) =

*a*(

*t*) cos(

*ϕ*(

*t*)), the 1D TK energy operator yields as output Ψ[

*s*(

*t*)] ≈ [

*a*(

*t*)

*ϕ̇*(

*t*)]

^{2}. Thus the local envelope

*a*(

*t*) and the instantaneous frequency

*ϕ̇*(

*t*) can be estimated using the Energy separation algorithm (ESA) [11]:

*I*(

*x*

_{1},

*x*

_{2}) as a 2D TK operator

*I*(

*x*

_{1},

*x*

_{2}) =

*a*(

*x*

_{1},

*x*

_{2}) cos(

*ϕ*(

*x*

_{1},

*x*

_{2})) the output of the 2D TK is given by Ψ[

*x*

_{1},

*x*

_{2}] ≈ [

*a*(

*x*

_{1},

*x*

_{2})

*ϕ̇*(

*x*

_{1},

*x*

_{2})]

^{2}. The spatially-varying amplitude

*a*(

*x*

_{1},

*x*

_{2}) can be interpreted as modeling local image contrast and the instantaneous frequency

*ω*(

*x*

_{1},

*x*

_{2}) = ∇

*ϕ*(

*x*

_{1},

*x*

_{2}) describes locally energy spatial frequencies. Applying Φ to

*∂I/∂x*

_{1},

*∂I/∂x*

_{2}and combining all energies, yields the 2D ESA [12]:

*a*(

**u**) and frequency vector

**w**= (

*w*

_{1},

*w*

_{2},...,

*w*)

_{n}*for any locally narrow band multidimensional signal such as:*

^{T}*s*(

**u**) ≃

*a*(

**u**) cos(

*ϕ*(

**u**)), where

**u**= (

*x*

_{1},

*x*

_{2},...,

*x*) and ∇

_{n}*ϕ*=

**w**(

**u**) ≃

**w**.

## 4. Ψ_{B} energy operator

Being related to the instantaneous cross correlation, the Ψ_{B} operator is well suited to measure the interaction between two signals [22]. This operator has found applications in time series analysis [23], data clustering [24], transient detection [22] and time delay estimation [19]. For two complex-valued signals *s*_{1}(*t*) and *s*_{2}(*t*), Ψ_{B} operator is defined as follows [16]:

_{C}is its associated quadratic form [16]. Compared to the cross-correlation function, Ψ

_{B}accounts for the relative changes of the signal by using the first and the second derivatives relative to time. We see in the following how the ΨB operator can be used to estimate the phase shift angle of the interferogram. This operator is helpful to construct a new correlation operator, which overcomes the limitations of the classic correlation function regarding the non stationary signals.

## 5. The proposed method

The TK operator and the Ψ_{B} correlator are combined for envelope extraction and phase retrieval of a WLSI signal. This algorithm, referred to as PETKB (Phase and Envelope estimated by the TK and the Ψ* _{B}* operators) is detailed as follows:

- Select a slice
*s*(*z*) =*s*(*x*_{0},*y*_{0},*z*) corresponding to the site*x*=*x*_{0},*y*=*y*_{0}. - Estimate the envelope
*Ĉ*(*z*) =*Ĉ*(*x*_{0},*z*) and the instantaneous frequency*ν̂*(*z*) =*ν*of the signal using the TK energy operator along the optical_{z}*z*-axis. - Smooth the obtained envelope using an interpolating spline function.
- Identify the local maxima(s)
*z*of the envelope_{max}*Ĉ*(*z*) on the optical*z*-axis. - Around each of these local maxima on the optical
*z*-axis:where- Update
*z*by ${z}_{\mathit{max}}^{\prime}={z}_{\mathit{max}}-\widehat{\theta}/(2\pi {\widehat{\nu}}_{z})$_{max} - If $\left|{z}_{\mathit{max}}^{\prime}-{z}_{\mathit{max}}\right|>\epsilon $ go back to step 5.

*ε*is a threshold value. Equation (9) gives the maximum of interaction between signals*s*(*z*) and*s*(_{θ}*z*) by the Ψ_{B}operator, located at*θ*=*θ̂*.

## 6. Results

The proposed demodulation scheme, PETKB, is illustrated on both synthetic and real data, and the results compared to those of the CWT, 1D TK, 2D TK, Five step phase shifting (FSPS) [25] and the modified Peak fringe scanning microscopy (PFSM) technique [26]. The PFSM algorithm takes into account the maximum value of the raw signal in the neighborhood of the local maxima and is initialized by the TK method. For CWT the complex Morlet mother wavelet is used with the scale depth set 4 and the carrier frequency *ν̂* (*z*) is estimated at the local maxima *z _{max}*. The 1D TK algorithm extracts the height of the surface by computing the local maxima(s) of the 1D envelope

*C*

_{1}(

*z*) along the optical

*z*-axis. The 2D TK algorithm first computes the 2D local envelope

*Ĉ*(

*x*,

*z*) using the 2D TK energy operator, and then extracts the 1D envelope

*Ĉ*(

*x*

_{0},

*z*) by projection along the optical z-axis. The derivatives of the TK operators are computed using the Savitzky-Golay filter [27]. This filter is well suited for smoothing noisy data and for efficient derivative estimation. It provides added noise robustness due to its low-pass filtering nature and produces smooth values of the derivatives. The FSPS algorithm consists in calculating the phase shift term

*θ*(Eq. (8)) using five consecutive samples along the optical

*z*-axis. We used a generalization of this method, to any

*T*value, given by [28]:

_{e}*ω*

_{0}=

*π*/2 corresponds to $\frac{4\pi {T}_{e}}{{\lambda}_{0}}=\pi /2$ and

*ω*

_{0}=

*π*/4 corresponds to $\frac{4\pi {T}_{e}}{{\lambda}_{0}}=\pi /4$ (for example an average wavelength of

*λ*

_{0}= 640 nm corresponds to

*T*= 80 nm for

_{e}*ω*

_{0}=

*π*/2 and

*T*= 40 nm for

_{e}*ω*

_{0}=

*π*/4). The intensities

*I*

_{0},

*I*

_{−2},

*I*

_{−1},

*I*

_{1},

*I*

_{2}represent respectively the intensity value at

*z*=

*z*and their values around

_{max}*I*

_{0}spaced by

*T*. The intensity

_{e}*I*

_{0}has been initialized by identifying the local maximum along the optical

*z*-axis using the 1D TK method. Once

*θ*is estimated,

*z*is updated by ${z}_{\mathit{max}}^{\prime}={z}_{\mathit{max}}-\widehat{\theta}/(2\pi {\widehat{\nu}}_{z})$ using an estimated value of the carrier frequency

_{max}*ν̂*.

_{z}The PETKB method is first tested on a synthetic interferometric image shown in Fig. 1(a) which is typical of that found on a transparent surface on a reflecting substrate resulting in two reflected signals (two layers). The threshold *ε* is set to 1 nm. The results in Figs. 1(b) and 1(c) show respectively the reference surface shape and the profile of a noisy signal along the optical *z*-axis with an additive offset low frequency signal sampled at *T _{e}* = 40 nm. Figures 1(d) and 1(e) show respectively a signal

*s*(

_{θ}*z*) superimposed on the original one (the offset is removed) and the estimated envelope

*Ĉ*(

*z*) along the optical axis, before and after the step 5 of the PETKB. This highlights the ability of the local phase retrieval to correct the position of the maximum first detected by the TK1D. Table 1 summarizes the results obtained for different Gaussian noise levels with

*T*= 80 and

_{e}*λ*

_{0}= 640 nm. We show the rate of relative errors compared with the reference surfaces for both surfaces: ”surface1” corresponds to the highest surface i.e, the higher level signal, and ”surface2” corresponds to the lower surface i.e, the weaker level optical signal. Clearly, in the noise free case the best results are provided by the PETKB, CWT and FSPS methods. In a noisy environment (10% or 20%) the PETKB performs slightly better than the CWT method and outperforms the other methods. Figure 2 shows that the profiles of surface1 and surface2 are well extracted by PETKB and FSPS methods [Figs. 2(a) and 2(b)], while 1D TK and 2D TK yield the same profiles but with noticeable fluctuations [Figs. 2(c) and 2(d)]. For surface2 the PETKB [Fig. 2(a)] performs slightly better than the FSPS [Fig. 2(b)]. The robustness of the methods has also been studied with respect to the fluctuations of the central carrier frequency value, which is also equivalent to assume that the

*T*measurement is not constant (for example a piezoelectric stepper leading to an inaccuracy in sampling). In this case randomized carrier frequencies around the mean value with a 5% margin interval vs additive Gaussian noise has been chosen. Results reported in Table 2 show that, on average, our method outperforms the other methods. Also the PFSM algorithm appears to be less robust than the PETKB and the CWT methods, specially for surface2 where it yields the highest error rate values. Furthermore, we have examined the impact of smaller sampling rates on the surface extraction. Table 3 shows the results for

_{e}*T*= 40 nm and with a fixed carrier frequency. The PETKB method consistently outperforms the other methods and particularly the CWT method for both the noiseless and noisy conditions. Table 4 shows the results for a non-constant carrier frequency (with a 5% margin). Again, as in Table 3 the same conclusions can be drawn. Even for a smaller sampling rate, we still show the effectiveness of the PETKB method in term of the error rate of surface extraction. We observe from the results of Tables 1 and 4 that, on average, the PETKB method performs better than the other methods. This performance can be attributed to the combined actions of the TK and the Ψ

_{e}_{B}operators. Indeed, the only action of the TK (1D TK or 2D TK) yields moderate performance in term of surface extraction and particularly in noisy conditions. The main advantages of the Ψ

_{B}operator are, first its efficient local ability as an instantaneous differential operator and second its ability as a correlator to accurately estimate the phase shift angle of the interferogram. In practice no more than 3 iterations are needed for the convergence of the PETKB algorithm, which is quite reasonable and helpful for the computing time performances. Simulations were carried out on a personal computer with a 2.13GHz Intel Core i3 and 4GB RAM. The computational time, for 256×256 interferometric image, of the CWT and the PETKB methods are compared. The PET KB takes 97 seconds while the CWT takes 95 seconds. We expect that if the scale depth of the CWT increases, then the associated CPU time will too. Even if in terms of computational time the two methods perform similarly, the PETKB provides better surface extraction results than the CWT. Finally, we show in Fig. 3 the demodulation results of real data with

*T*= 80 nm. These measurements were made on a system developed in our laboratory consisting of a Leitz-Linnik interference microscope with ×50 objectives (numerical aperture = 0.85) under white light illumination as described in [29]. The piezo for scanning the sample along the optical axis was the 100

_{e}*μ*m PIFOC model (from PI) with integrated LVDT position sensor to give a linear response and a sensitivity of 5 nm. Figures 3(a) and 3(b) display respectively the original image and its 2D gradient. This derivative method allows us to eliminate the ”batwing effect” or envelope skewing, caused by interference from the top and bottom of the steps for heights less than the coherence length of the light used. Figures 3(c)–3(f) compare the robustness of the profile extraction using the PETKB, FSPS, 1D TK and 2D TK methods. Except for the few artefacts, the extracted surface profile by the PETKB is more accurate than the ones obtained by the other methods, and especially the FSPS method.

## 7. Conclusion

In this paper a new fringe envelope retrieval technique in WLSI is introduced. This demodulation approach exploits the potential of the TK and the Ψ_{B} energy operators as instantaneous energy tracking tools. The TK envelope detection is combined with the TK frequency estimation to retrieve the envelope and instantaneous frequency simultaneously. Once the reference signal is extracted, a correlation technique based on the non-linear operator Ψ* _{B}* is used to identify the surface height information contained simultaneously in the phase and envelope of the signal along the optical

*z*-axis. Computations performed on synthetic and real data with

*T*= 80 nm and

_{e}*T*= 40 nm along the optical

_{e}*z*-axis, show the interest of this TK class of operators, and highlight their robustness to the changes in carrier frequencies, compared to the CWT, FSPS and PFSM methods. In terms of execution time and error rate of surface extraction results, the PETKB method is more competitive than the CWT approach. The promising results obtained are a motivation to take into account more complex information and to apply the 2D and/or 3D detection (correlation) associated with the 2D and/or 3D envelope/frequency estimation. In future work we plan to extend this approach to multidimensional signals.

## Acknowledgments

The authors would like to thank the reviewers for their fruitful suggestions and their comments that helped to improve the paper. Thanks also to Denis Montaner for useful discussions.

## References and links

**1. **C. O’Mahony, M. Hill, M. Brunet, R. Duane, and A. Mathewson, “Characterization of micromechanical structures using white-light interferometry,” Measurement Sci. Technol. **14**, 1807 (2003). [CrossRef]

**2. **K. Larkin, “Efficient nonlinear algorithm for envelope detection in white light interferometry,” JOSA A **13**, 832–843 (1996). [CrossRef]

**3. **P. de Groot, X. C. de Lega, J. Kramer, and M. Turzhitsky, “Determination of fringe order in white-light interference microscopy,” App. Opt. **41**, 4571–4578 (2002). [CrossRef]

**4. **S. Ma, C. Quan, R. Zhu, C. Tay, L. Chen, and Z. Gao, “Micro-profile measurement based on windowed Fourier transform in white-light scanning interferometry,” Opt. Comm. **284**, 2488–2493 (2011). [CrossRef]

**5. **P. Sandoz, “Wavelet transform as a processing tool in white-light interferometry,” Opt. Lett. **22**, 1065–1067 (1997). [CrossRef]

**6. **M. Li, C. Quan, and C. Tay, “Continuous wavelet transform for micro-component profile measurement using vertical scanning interferometry,” Opt. Laser Technol. **40**, 920–929 (2008). [CrossRef]

**7. **Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” App. Opt. **43**, 2695–2702 (2004). [CrossRef]

**8. **H. Niu, C. Quan, and C. Tay, “Phase retrieval of speckle fringe pattern with carriers using 2D wavelet transform,” Opt. Lasers Eng. **47**, 1334–1339 (2009). [CrossRef]

**9. **J. Kaiser, “On a simple algorithm to calculate the energy of a signal,” in Proc. ICASSP, (1990), pp. 381–384.

**10. **P. Maragos and A. Potamianos, “Higher order differential energy operators,” IEEE Sig. Proc. Lett. **2**, 152–154 (1995). [CrossRef]

**11. **P. Maragos, T. Quatieri, and J. Kaiser, “Energy separation in signal modulations with applications to speech analysis,” IEEE Trans. Sig. Proc. **41**, 3024–3051 (1993). [CrossRef]

**12. **P. Maragos and A. Bovik, “Image demodulation using multidimensional energy separation,” J. Opt. Soc. Am. A **12**, 1867–1876 (1995). [CrossRef]

**13. **A. Boudraa, F. Salzenstein, and J. Cexus, “Two-dimensional continuous higher-order energy operators,” Opt. Eng. **44**, 7001–7010 (2005).

**14. **K. Larkin, “Uniform estimation of orientation using local and nonlocal 2-D energy operators,” Opt. Express **13**, 8097–8121 (2005). [CrossRef]

**15. **F. Salzenstein, A. Boudraa, and T. Chonavel, “A new class of multi-dimensional Teager-Kaiser and higher order operators based on directional derivatives,” Multidimensional Sys. Sig. Proc. **24**, 543–572 (2013). [CrossRef]

**16. **J. Cexus and A. Boudraa, “Link between cross-Wigner distribution and cross-Teager energy operator,” Elect. Lett. **40**, 778–780 (2004). [CrossRef]

**17. **A. Boudraa, “Relationships between Ψ_{B}-energy operator and some time-frequency representations,” IEEE Sig. Proc. Lett. **17**, 527–530 (2010). [CrossRef]

**18. **A. Boudraa, T. Chonavel, and J. Cexus, “Ψ_{B}-energy operator and cross-power spectral density,” Sig. Proc. **94**, 236–240 (2014). [CrossRef]

**19. **A. Boudraa, J. Cexus, and K. Abed-Meraim, “Cross-Ψ_{B}-energy operator-based signal detection,” JASA **123**, 4283–4289 (2008). [CrossRef]

**20. **F. Salzenstein, P. Montgomery, D. Montaner, and A. Boudraa, “Teager-Kaiser energy and higher order operators in white light interference microscopy for surface shape measurement,” J. Appl. Sig. Proc. **17**, 2804–2815 (2005). [CrossRef]

**21. **S. Petitgrand, “Méthodes de microscopie interférométrique 3D statiques et dynamiques pour la caractérisation de la technologie et du comportement des microsystèmes,” Ph.D. thesis (2005).

**22. **A. Boudraa, S. Benramdane, J. Cexus, and T. Chonavel, “Some useful properties of cross-Ψ_{B} energy operator,” Int. J. Electron. Comm. **63**, 728–735 (2009). [CrossRef]

**23. **A. Boudraa, J. Cexus, M. Groussat, and P. Brunagel, “An energy-based similarity measure for time series,” Adv. Sig. Proc. **135892**, 1–8 (2008).

**24. **W. Zhang, C. Liu, and H. Yan, “Clustering of temporal gene expression data by regularized spline regression an energy based similarity measure,” Patt. Recong. **43**, 3969–3976 (2010). [CrossRef]

**25. **P. Hariharan, B. Oreb, and T. Eiju, “Digital phase-shifting interferometry: a simple error-compensating phase calculation algorithm,” Appl. Opt. **26**, 2504–2506 (1987). [CrossRef]

**26. **P. Montgomery and J. Fillard, “Peak fringe scanning microscopy (pfsm): submicron 3d measurement of semiconductor components,” Interferometry: Techniques and Analysis pp. 12–23 (1755).

**27. **R. Schafer, “What is a savitzky-golay filter?[lecture notes],” IEEE Sig. Proc. Mag. **28**, 111–117 (2011). [CrossRef]

**28. **J. Estrada, M. Servin, and J. Quiroga, “A self-tuning phase-shifting algorithm for interferometry,” Opt. Express **18**, 2632–2638 (2010). [CrossRef]

**29. **E. Halter, P. Montgomery, D. Montaner, R. Barillon, M. D. Nero, C. Galindo, and S. Georg, “Characterization of inhomogeneous colloidal layers using adapted coherence probe microscopy,” Appl. Surf. Sci. **256**, 6144–6152 (2010). [CrossRef]