EURASIP Journal on Applied Signal Processing 2005:17, 2804–2815 c ○ 2005 Hindawi Publishing Corporation Teager-Kaiser Energy and Higher-Order Operators in White-Light Interference Microscopy for Surface Shape Measurement

In white-light interference microscopy, measurement of surface shape generally requires peak extraction of the fringe function envelope. In this paper the Teager-Kaiser energy and higher-order energy operators are proposed for efficient extraction of the fringe envelope. These energy operators are compared in terms of precision, robustness to noise, and subsampling. Flexible energy operators, depending on order and lag parameters, can be obtained. Results show that smoothing and interpolation of envelope approximation using spline model performs better than Gaussian-based approach.


INTRODUCTION
Different signal processing methods have been proposed in coherence probe microscopy (CPM), also known as whitelight scanning interference microscopy (WLSI), for roughness surface measurement [1,2]. A basic problem in CPM consists in developing an efficient and precise peak detection process of the fringe envelope that corresponds to the axial position of the surface. Moreover, fast and robust methods against noise are required. Most of the methods are based on an AM-modulated signal model, which represents the variation in light intensity measured along the optical axis of an interference microscope. Envelope detection can be performed using a centroid calculation [3], a demodulation procedure [1,2], or measurement of the fringe visibility [4] at a given pixel along the optical axis, z. The Demodulation method [2] requires carrier frequency information, while the five sample adaptive (FSA) method, proposed by Larkin [4], detects the peak by using only five adjacent samples along the optical axis. Methods based on Fourier transform [1,5] or wavelets [6] have also been proposed but require more involved calculation means. The FSA method corresponds to the Teager-Kaiser energy [7], applied to the differentiated signal (in order to eliminate the added lowfrequency components) [8]. The TKEO (Teager-Kaiser energy operator) has found applications in speech processing [9], image processing [10], and pattern recognition [11]. This operator is an energy-tracking operator which does not measure the signal energy but the energy of the source or the system that produces this signal. For demodulating an AM-FM signal into its varying amplitude and instantaneous frequency, Maragos et al. [12] have used the TKEO to develop the energy separation algorithm (ESA) to isolate the AM and FM components. Maragos and Potamianos have proposed an extension and generalization of the TKEO called higher-order energy operators [13]. Some highorder methods were applied in the signal processing field. For example, higher-order instantaneous moments have been introduced in [14] and further used by Barbarossa et al. [15]. Also the quadratic Wigner distribution has been extended to higher-order and found applications in signal processing [16]. However, to our knowledge, there is few work in the literature based on the generalized higher-order operators introduced in [13]. In this paper the discrete higherorder energy operators are applied in CPM for surface profiling. Attention is focused on the envelope detection technique instead of the phase extraction method for measuring surface height. It should be noted that a general technique has also been proposed where we first fit the sampled signal with a spline in order to apply the continuous TKEO energy operator and improve the noise rejection [17]. Here, we propose the use of spline smoothing and interpolation methods after the envelope detection, instead of using a Gaussian shape as many authors do. In Section 2, we present the general model of an interferometric signal. We briefly expose the higher-order energy operators technique in Section 3. We also propose an improvement in the choice of the sampling step, exhibiting a relationship between the sample period, the lag parameter, and the order of the operator (Section 3.3). In Section 4 we present results of tests of the operators on experimental data according to various criteria such as the quality of detection and robustness to the spline or Gaussian techniques (Section 4.3), noise (Section 4.4), and subsampling (Section 4.5).

General model of a signal along the optical z-axis
In using white-light illumination of a sample in an interference microscope, the visibility of the fringes drops off rapidly from a maximum value at the minimum OPD (optical path difference) [2]. Figure 1 shows a typical intensity signal obtained from a CCD (couple-charged device) sensor as the OPD is varied through focus, at a given point (x, y) on the material surface. The signal can be approximated by a modulated sinusoid [4]: s(x, y, z) = a(x, y, z) The factor b(x, y) is proportional to the reflected beam intensity and c(z − 2 · h(x, y)) corresponds to the envelope along the z-axis. The position along the optical axis at the peak of the fringe envelope corresponds to the position of the surface measured at that point. In certain conditions, the offset a(x, y, z) varies slowly. The common technique consists in differentiating the signal in order to eliminate this lowfrequency component.

Spline smoothing and interpolation of the envelope
In general once the envelope is detected, the peak is located using the least squares fitting (LSF) technique along z-axis. In the peak region, the calculated envelope can be approximated by a Gaussian shape exp(−βz 2 ). The envelope is often not symmetrical along the optical axis. To avoid this problem, we propose a spline-based approach which decomposes the envelope into piecewise polynomials with pieces that are smoothly connected together. This approach combines two parameters: the quadratic distance and the smoothness η [18]. The spline function minimizes λ + (1 − λ)η. The more the parameter λ tends to 1, the greater is the tendency towards the least squares solution. Practically, a value of λ = 0.1 is chosen (see Section 4.3). The last step consists in performing an envelope interpolation using the cubic spline method. It should be noted that enough samples should be chosen in order to give a high enough sensitivity at the peak position.

Continuous and discrete higher-order operators
For a continuous real-valued signal s(t) the TKEO, Ψ, is defined via The TKEO is useful to provide the instantaneous frequency and envelope information from a monochromatic AM-FM signal [12]. An approximation of the derivatives by onesample differences provides the discrete-time counterpart of the TKEO [7]: For an AM monochromatic signal s(n) = a(n) cos(Ωn), with the assumption that the envelope varies more slowly than the carrier signal, the discrete TKEO of s(n) yields Ψ[s(n)] a 2 (n) sin 2 (Ω) [12]. An efficient demodulation is thus provided with only three samples. A generalization of the continuous TKEO has been proposed in [13]. This operator can be seen as a particular case of the continuous kth-order differential energy operator (CEO) Ψ k defined by Hence, the TKEO corresponds to k = 2. There is a recurrence step algorithm between one kth-order operator and the previous orders [13]: It is also possible to demodulate monochromatic signals combining different kth-order operators. Finally, one can unify the discrete versions of Teager-Kaiser and higher-order operators under the same class of discrete energy operators (DEO) [13] Ψ km s(n) = s(n)s(n + k) − s(n − m)s(n + k + m).
It can be noticed that the particular case of Ψ 01 defines a generalized TKEO when m > 1 and the classical TKEO when m = 1. More generally, k represents the order of the operator, while m adjusts the distance between the chosen samples.
We have also recently proposed the 2D discrete extensions of the higher-order operators [19]. Moreover, a cross-discrete higher-order energy operator (CDEO) Φ km can be defined between two signals s and w [13]: In Section 4.4, we study the DEO and CDEO in the presence of noise. It can be noticed that the CDEO provides interesting properties in the presence of two uncorrelated signals s and w, while the higher-order DEO better filters the uncorrelated noise. Some properties of the continuous cross-Teager energy operators and their complex extensions have been studied in [20]. We now explore an application of the DEO to the demodulation of interference signals.

Demodulation of real interference signals
A set of discrete ESA (DESA) has been proposed by Maragos et al. [12] for tracking the instantaneous frequency and envelope of monocomponent AM-FM signal. We consider a local AM signal s(n) = a(n) cos(nΩ + θ), with Ω = ωT e , T e being the sampling step. A local constant amplitude The backward derivative y(n) = (s(n)−s(n−1)) of the signal is Hence, Applying (8) and (10), we derive the formula Moreover, where P 2m+k (x) is a 2m + k degree polynomial in cos Ω. Finally, from the estimated Ω value, we deduce the amplitude for the higher-order DESA: Various higher-order DESA can also be derived from other continuous derivative approximations (backwards, forwards, or symmetric differences). This leads to the following observations.
(1) The higher-order DESA is effective for local AMsignals, unless the denominators in (11) and (13) are null. For example, an estimation of an infinite amplitude is possible for the zeros of P 2m+k (x), even though the real signal has a finite energy. It is shown that the TKEO is not effective for certain signals and can violate the conditions of physical continuity required for the amplitude and detected phase [21]. Although the work mentioned deals with the continuous TKEO and argues that the results are acceptable for narrowband signals, such a situation can occur for the DEO. (2) The higher-order DESA is complex and not competitive in terms of computing time.
To overcome the previous problems, a more efficient and faster algorithm based on a signal interference model and adapted sampling step is proposed. In white-light interferometry, the signal can be modeled by an AM signal where the envelope is Gaussian. We compute the DEO of the discrete AM signal s(n) = a(n) cos(ωnT e ) and a(n) = e −α(nTe−β) 2 . The maximum value of this function occurs at β/T e and corresponds to the surface position. We have In surface detection by CPM, the step T e is a multiple of 10 nm, and α, m, k are small enough to write the following approximations: Using relation (15), Ψ km [a(n) cos(ωnT e )] can be approximated as follows: Ψ km a(n) cos ωnT e + θ a 2 n + k 2 sin mωT e sin (m + k)ωT e .
Finally, |Ψ km [a(n) cos(ωnT e + θ)]| is proportional to the envelope a(n) translated by k/2 to the left-hand side. This fact does not affect the results provided that relative heights of the surface are measured and not absolute positions. An interesting point concerns the generalized Teager-Kaiser approach (6) with k = 0, for which the resulting envelope is not translated. Hence, under certain conditions relying on the parameters T e , k, and m (see Section 3.3), such an operator is able to perform a fast demodulation procedure with only a few samples.

Orders and optimal sampling time
Unless the condition (17) is satisfied, the detected envelope is attenuated: sin mωT e sin (m + k)ωT e = 1.
In particular, (18) and (19) give an optimal detection: where p, q ∈ N. Dividing (19) by (18) yields This factor should contain primary odd numerator and denominator terms. A particular solution consists in choosing an even integer k and odd integer for m. We notice k = 0 is always a solution. For example, the Teager-Kaiser algorithm with k = 0 and m = 1 is a solution that corresponds to the FSA algorithm [4,8]. Moreover, (18) shows the link between the lag parameter m and the sampling time T e for a given carrier frequency ω = 2πν: Thus, an optimal sampling time for a given lag parameter can be computed, and inversely the lag parameter can be adjusted to take into account the sampling period. For instance, an average wavelength of 640 nm in WLSI corresponds to a value of ν = 1/320 nm −1 . Then, q = 0 gives one solution: T e = 80 nm, k = 0, m = 1. Inserting ν = 1/320 nm −1 , (21) becomes (1) If a sampling period T e is given, (a) compute m according to (22), choosing a value for q, (b) compute k according to (20), (c) while m is too big, go to (a) and compute m according to the nearest value T e of T e (see below). (2) If a lag parameter m is given, (a) choose q and compute T e according to the formula (22), (b) compute the order k according to (20), (c) while k is too big, go to (a) and change the q value.
The aim of our work is to choose optimal values of T e , m, and k in order to preserve the information. This also demonstrates the possibility of undersampling the signal and preserving the envelope. The general procedure is shown in Algorithm 1. In general, for steps 1(c) and 2(c), any values of k and m which are too big are avoided. Actually, the detected envelope is sensitive to the artifacts outside the frequency peak, due to a low SNR in that area. To obtain adequate parameters, the choice of m is based on the nearest integer T e of T e . For example, a chosen sampling period of T e = 85 nm gives a value of m = 16, while periods T e = 90 nm and T e = 80 nm give, respectively, m = 8 and m = 1. Finally, when T e = 85 nm, the choice of the parameter m is based on the nearest value T e = 80 nm.

Experimental context
The proposed algorithms have been first applied to the profiling of a step etched in silicon (Figure 2a) using the interference microscopy system described in [22]. A series of xy images were scanned along z and a single xz image retrieved for use in the present work. The image corresponds to a modulated signal (1D or 2D) typical of such an optical system. The peak of the fringe envelope gives the position of the surface.
Converting the OPD to time along the z-axis gives a sampling  Figure 2a is 128 × 130 pixels. Secondly, in order to compare these procedures, an objective measure of detection quality is introduced. The assumption is based on the regular shape of the studied material: this leads us to compare the final surface with two regular lines corresponding, respectively, to the top level and the bottom level of the step. Thus, the profile in Figure 2b represents one possible reference profile. The reference is computed using a linear regression algorithm. Once the reference shape is obtained, the quality detection factor (error rate) consists of the absolute mean error between the processed and the reference surface. The measured surface appears to be sloping because of the introduction of tilt fringes by slightly tilting the reference mirror. The reasons for using a sloping sample instead of one that is parallel to the reference mirror is to observe how the envelope detection algorithm behaves towards variations in the positions of the sampling points on the envelope, which is particularly important in the subsampling case (see Section 4.5).

Surface profiling and optimal orders
A sampling step of T e = 80 nm leads to the following rule: show the spline interpolation effect on the same operators.

Surface detection using spline and Gaussian techniques
The best previous operators, Ψ 01 , Ψ 03 , Ψ 21 , and Ψ 23 , are selected to study their behavior regarding the spline and Gaussian approaches. The robustness of Gaussian and spline methods are measured relatively to window size and smoothing parameter, respectively. For Gaussian approximation, in practice a window W of several points around the top of the envelope is used. The envelope is obtained using LSF method. Figure 6a (resp., Figure 6b) compares the error rate (nm) for Ψ 01 and Ψ 21 (resp., Ψ 03 and Ψ 23 ) according to the size of W. According to Figure 6a, the Gaussian approach performs better for Ψ 01 (minimum error rate of 12.2 nm), but gives the worst results for Ψ 23 (minimum error rate of 16.8 nm). Figure 7a (resp., Figure 7b) compares the error rate (nm) using splines for the operators Ψ 01 and Ψ 21 (resp., Ψ 03 and Ψ 23 ) according to the smoothing parameter λ. In any case, the best smoothing parameter λ equals 0.1. One may notice that the spline-based approach is more adapted to the  generalized operators Ψ 21 and Ψ 23 . Moreover these energy operators are robust in regards to variations of λ: the position error drops faster when using Ψ 01 or Ψ 03 . Remark that the operator Ψ 23 is the most competitive when using splines. Finally, the previous results can be explained as follows: the more the profile tends to be asymmetric-which is the case in the presence of higher-order functions-the less adapted is the Gaussian model. In the next section, we show the mutual influence between the symmetric assumption and the noise level.

Robustness to noise data
The robustness to noise is an important aspect of any envelope detection algorithm used in practical CPM systems for surface roughness measurement. Indeed, the quality of signals can vary greatly depending on the nature of the sample. For a given signal s(n) with an added noise w(n), the output of Ψ km of the noisy signal s(n) + w(n) is given by The algorithms are studied in the presence of noise according to different envelope estimations. A white noise w(n) with a relative standard deviation of 2% (i.e., a value of 5.0 for 255 grey levels) is added to the original measurement data. Figure 8 shows a square of a typical interpolated envelope |Ψ km [s(n)]| superimposed with the crossterms |Φ km [s(n); w(n)] + Φ km [w(n); s(n)]|. The cross-energy terms of (26) are lower than the function Ψ km [s(n)]. So, the CDEO seems to follow the lack of correction between the samples s(n) and w(n). As shown in Figure 8, the CDEO are comparable for k = 0 and k = 2. On the other hand, the noise being uncorrelated, we expect Ψ 0m [w(n)] w(n) 2   Figure 9a (resp., Figure 9b) shows a comparison between the DEO Ψ 01 and Ψ 21 (resp., Ψ 03 and Ψ 23 ) in the presence of a single white noise signal w(n) along the optical axis. The energy operators Ψ 21 and Ψ 23 clearly provide a lower noise level with a zero mean, whereas Ψ 01 and Ψ 03 are proportional to w 2 (n). Thus, the higher-order functions seem to better filter the noise components than the lower-order operators. Finally, we compared both Gaussian and spline methods in  the presence of noise. In order to process the randomized effect due to the noise, each result is based on an average of 50 different simulations. Figure 10a (resp., Figure 10b) compares the error rate (nm) obtained by Gaussian approach, for Ψ 01 and Ψ 21 (resp., Ψ 03 and Ψ 23 ) according to the size of W. In the same manner, Figure 11a (resp., Figure 11b) compares the error rate (nm) obtained by spline approach for Ψ 01 and Ψ 21 (resp., Ψ 03 and Ψ 23 ) according to λ. For Gaussian method, the performances are sensitive to the window size. This is particularly true for the DEO Ψ 03 and Ψ 23 . Again, the worst position errors are obtained by the third-order DEO Ψ 23 (5 nm higher), while the DEO Ψ 21 is more stable.
The operator Ψ 01 gives the best results for the Gaussian approximation (minimum error rate of 17 nm versus 18.5 nm for splines). For spline approach, the performances related to the smoothness parameter confirm the previous results: the optimal value of λ equals 0.1, while Ψ 21 and Ψ 23 are more robust to λ. The operator Ψ 23 is the most competitive when splines are used and gives here the best results (minimum error rate of 15.3 nm versus 21 nm for Gaussians). Finally, the spline technique seems to be well adapted or more competitive for the higher-order operators Ψ 21 and particularly Ψ 23 , while the Gaussian approach suits better both other operators. However, the more the detected profile along the optical z-axis is regular, the better the surface detection will be. Basically the splines based only on the distance criteria are not adapted to this problem because they tend to follow the noise. The global advantage of the Gaussian approach consists in finding a given regular symmetric shape and this is helpful to overcome the presence of noise, particularly when such operators like Ψ 01 or Ψ 03 are more sensitive to the noise (Figure 9). This is also its limit because a strong hypothesis concerning the shape of the profile might miss the correct surface position in the presence of an asymmetric operator such as Ψ 23 .

Sample step
Using the method described in Section 3.3, possible values of k and m in regards to different sampling steps T e , which are multiples of 20 nm (i.e., 20·p nm where p is an integer), are computed. The performances of the associated functions Ψ km were then tested. Some results are summarized in Figures 12, 13, and 14. For sampling steps of T e greater than 80 nm, the DEO Ψ 0m gives the best results for different lag parameters m when T e = 100 nm (m = 4, Figure 12a), T e = 120 nm (m = 2, Figure 12b Figure 14a). Moreover, the energy operators Ψ 0m are computationally more interesting, because only three samples are required. Other higher orders (k = 0) give errors in position between 15 and 25 nm for sampling step T e = 120 nm (e.g., (k, m) = (4, 2), Figure 13a) or T e = 240 nm (Ψ 21 and Ψ 25 give respective position errors of 22.8 and 24.6 nm). However, improvements in the precision of the acquisition system, and the use of adapted prefiltering could improve these results for the undersampling case.

CONCLUSION
In this work, higher-order energy operators based on the Teager-Kaiser function for surface relief analysis using interference microscopy are used. These functions are well adapted to this problem since the data consists of signals modulated along the optical z-axis. The discrete versions of these functions are described by a set of two parameters, one parameter corresponding to the order and the other one to a lag, from which the choice of the adjacent samples around the peak of the envelope can be adjusted. In particular these   operators generalize the FSA technique. Secondly, instead of using only an approach based on a Gaussian assumption to extract the envelope, the cubic splines based on distance and smoothness criteria are proposed. Both the Gaussian and spline methods have been tested on a step etched in silicon and sloping with respect to the reference mirror. For the higher-order functions, the spline technique is competitive in the presence of noise and should offer more flexibility for a nonsymmetrical envelope. A general rule relating the order, the lag parameter, and the sampling step is introduced. For a given sampling period, there exists an optimal choice of these energetic operators which leads to flexible operators in terms of sampling step, precision, and noise level. This flexibility is important for choosing the optimal algorithm for a given sample and measurement conditions.