Local frequency and envelope estimation by Teager-Kaiser energy operators in white-light scanning interferometry

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 crossenergy 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 OCIS codes: (100.0100) Image processing; (100.3175) Interferometric imaging; (100.4992) Pattern, nonlinear correlators; (100.5070) Phase retrieval. 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). 2. K. Larkin, “Efficient nonlinear algorithm for envelope detection in white light interferometry,” JOSA A 13, 832– 843 (1996). 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). 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). 5. P. Sandoz, “Wavelet transform as a processing tool in white-light interferometry,” Opt. Lett. 22, 1065–1067 (1997). 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). 7. Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” App. Opt. 43, 2695–2702 (2004). 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). 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). 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). #211263 $15.00 USD Received 12 May 2014; revised 30 Jun 2014; accepted 3 Jul 2014; published 22 Jul 2014 (C) 2014 OSA 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018325 | OPTICS EXPRESS 18325 12. P. Maragos and A. Bovik, “Image demodulation using multidimensional energy separation,” J. Opt. Soc. Am. A 12, 1867–1876 (1995). 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). 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). 16. J. Cexus and A. Boudraa, “Link between cross-Wigner distribution and cross-Teager energy operator,” Elect. Lett. 40, 778–780 (2004). 17. A. Boudraa, “Relationships between ΨB-energy operator and some time-frequency representations,” IEEE Sig. Proc. Lett.17, 527–530 (2010). 18. A. Boudraa, T. Chonavel, and J. Cexus, “ ΨB-energy operator and cross-power spectral density,” Sig. Proc. 94, 236–240 (2014). 19. A. Boudraa, J. Cexus, and K. Abed-Meraim, “CrossΨB energy operator-based signal detection,” JASA 123, 4283–4289 (2008). 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). 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). 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). 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). 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). 28. J. Estrada, M. Servin, and J. Quiroga, “A self-tuning phase-shifting algorithm for interferometry,” Opt. Express 18, 2632–2638 (2010). 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).


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][11][12][13][14][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.

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]: where 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.

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 whereẋ(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 Thus the local envelope a(t) and the instantaneous frequencyφ (t) can be estimated using the Energy separation algorithm (ESA) [11]: This operator has been extended to signals of higher dimensions, rendering it applicable to real-valued gray level images I(x 1 , x 2 ) as a 2D TK operator In a similar manner to the 1D TK operator, for a narrowband image 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]: Note that the FSA method [2] corresponds to the discrete TK operator, applied to a differentiated signal. The 2D TK operator has shown its efficiency for image demodulation and particularly for local envelope extraction [12][13][14]. This result has been recently extended to multidimensional signals by introducing a nD compact form of the TK operator using directional derivatives [15]. In addition, it has been shown that this new operator is able to track the envelope a(u) and frequency vector w = (w 1 , w 2 , . . . , w n ) T for any locally narrow band multidimensional signal such as:

Ψ 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]: where This operator is a symmetric bilinear form and Ψ 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.

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: 2. Estimate the envelopeĈ(z) =Ĉ(x 0 , z) and the instantaneous frequencyν(z) = ν z of the signal using the TK energy operator along the optical z-axis.
3. Smooth the obtained envelope using an interpolating spline function.
4. Identify the local maxima(s) z max of the envelopeĈ(z) on the optical z-axis.

Around each of these local maxima on the optical z-axis:
• Fit the model, s θ (z), to the extracted parameters: • Estimate the phase shift angle, θ , with Ψ B as follows: • Update z max by z where ε is a threshold value. Equation (9) gives the maximum of interaction between signals s(z) and s θ (z) by the Ψ B operator, located at θ =θ .

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 e value, given by [28]: where ω 0 = π/2 corresponds to 4πT e λ 0 = π/2 and ω 0 = π/4 corresponds to 4πT e λ 0 = π/4 (for example an average wavelength of λ 0 = 640 nm corresponds to T e = 80 nm for ω 0 = π/2 and T e = 40 nm for ω 0 = π/4). The intensities I 0 , I −2 , I −1 , I 1 , I 2 represent respectively the intensity value at z = z max and their values around I 0 spaced by T e . The intensity I 0 has been initialized by identifying the local maximum along the optical z-axis using the 1D TK method. Once θ is estimated, z max is updated by z ′ max = z max −θ/(2πν z ) using an estimated value of the carrier frequencyν 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 e = 80 and λ 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  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 e 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 T e = 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-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 Ψ 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 256x256 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 e = 80 nm. These measurements were made on a system developed in our laboratory consisting of a Leitz-Linnik interference microscope with x50 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µ 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.

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 e = 80 nm and T e = 40 nm along the optical 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.