1. Introduction

Phase shifting interferometry (PSI) is a widely used technique for topography measurement of specular surfaces with sub-nanometer resolution [1, 2]. Usually, several (at least three) linearly phase shifted interferograms are captured from which the phase of the wavefront is calculated. A large number of different algorithms for evaluation of linearly phase shifted interferograms exist [1, 3–5], but non of them is directly applicable to sinusoidally phase shifted interference (SinPSI) signals.

Most SinPSI evaluation algorithms calculate the phase related to a height value from the ratio of the amplitudes of the harmonics of the interference signal. The first SinPSI algorithm was introduced by Sasaki et al. in 1986 [6, 7]. It uses the first two harmonics of the interference signal captured at four sampling positions within a complete sinusoidal phase shift cycle. Later, a similar technique for the detection of the first two harmonics in sinusoidally modulated interference signals by an analog electronic circuit was used by Minoni et al. [8]. Dubois [9] also uses four sampling points but in combination with the integrating bucket technique at particular phase shift amplitude and initial phase offset. Fallagis et al. [10] increase the phase shift amplitude to achieve a nearly linear movement of the phase shifter at the zero crossings of the sinusoidal phase shift signal and use linear PSI evaluation algorithms with samples acquired at equidistant spatial phase steps. Schulz et. al [11], Knell et al. [12] and Tereschenko et al. [13] use a technique similar to Fallagis at higher phase shift amplitudes but in combination with a lock-in phase detection, which uses all sampling points around the zero crossings at falling as well as rising flanks of the sinusoidal phase shift signal. De Groot [14] extends the original Sasaki method to arbitrary number of harmonics and sampling points with the possibility to increase the robustness of the algorithm even if uncertainties in phase shift amplitude occur, by applying weights to the different harmonics. In this article, we focus on the method of de Groot and show how it can be further generalized. In particular, we show how it can be implemented for practical applications and how by some extensions a significantly higher position data rate can be achieved.

We begin in chapter 2 with a brief review of the original algorithm, then in chapter 3 we explain our extensions and generalizations for the practical implementation of the algorithm. Later, in section 4 we extend the algorithm to a temporal high-speed implementation. In section 5 we show results based on our extensions and finally, in section 6 the paper ends with a conclusion and an outlook to future work.

2. Original algorithm

A monochromatic sinusoidally phase modulated interference signal is given by the following equation:

To take into account the integration time during which each individual intensity value is recorded, we need to integrate over a certain time [14]:

The introduced attenuation factor $B\left(n\right)$ considers only the intensity averaging, also known as integrating bucket. In general, electronic devices used for the signal acquisition are equipped with (anti-aliasing) low-pass filters, which additionally introduce an amplitude reduction as well as a phase shift for higher harmonics in the interference signal. For the sake of simplicity, it is assumed that the cutoff frequency of the measuring devices is high enough not to influence the acquired interference signal and especially those harmonics, which are used for signal processing.

The idea of the algorithm is based on the calculation of the ratio of at least two harmonics, or more precisely their differently weighted relative strengths, so that the phase Θ

can be calculated with the phase quadrature formula:

$Q_{\text{odd},\text{even}}^{\left(n\right)}$ are independent from $\sin (\mathrm{\Theta })$ and $\cos (\mathrm{\Theta })$ strengths of the harmonics of the interference signal for a particular phase shift amplitude a₀.

The algorithm calculates a single value Θ from the interference signal sampled at P temporally equidistant positions corresponding to one period of the sinusoidal excitation signal $\varphi \left(t\right)$. Figure 1 shows both signals for the phase shift amplitude $a=\pi$, the height dependent phase $\mathrm{\Theta }=\pi /2$ and the phase shift offset φ = 0.

 

Fig. 1 Sinusoidally phase modulated interference signal (blue curve) sampled at the temporally equidistant points with the sampling interval Δα during one period of the phase shift function ϕ(t) = π cos(2πt).

Download Full Size | PDF

To detect the height dependent phase Θ by use of Eq. (8) one has to perform a frequency analysis of the interference signal for odd and even harmonics separately. De Groot [14] suggests to adjust the phase offset to φ = 0, so that the sampling positions ${\alpha }_j=j\mathrm{\Delta }\alpha +{\alpha }_0$ with $\mathrm{\Delta }\alpha =2\pi /P$ and ${\alpha }_0=\pi /P$ of the interference signal are symmetrical to the middle of the sampling period (see Fig. 1). In this case one can use a normalized real-valued Fourier (cosine) transform

The calculation of $\sin (\mathrm{\Theta }){\displaystyle \sum }_{n=1,3,\mathrm{...}}^{n_{max}}{\gamma }_n{Q}_{\text{odd}}^{\left(n\right)}\text{and}\cos (\mathrm{\Theta }){\displaystyle \sum }_{n=2,4,\mathrm{...}}^{n_{max}}{\gamma }_n{Q}_{\text{even}}^{\left(n\right)}$ according to Eq. (8) is performed by the following formulas:

The required normalization constants ${\mathrm{\Gamma }}_{\text{odd},\text{even}}$ in Eq. (8) are calculated by

To improve the robustness of the algorithm, de Groot proposes to use so-called filter functions to adjust the weighting factors γ_n:

The design rules for the selection of γ_n are the following: "First, the normalized filter functions $F_{\text{odd}}/{\mathrm{\Gamma }}_{\text{odd}}$ and $F_{\text{even}}/{\mathrm{\Gamma }}_{\text{even}}$ should be near maximum, matched in value, and matched in derivative at the design amplitude a₀. Second, the filter functions should have low values far from the design amplitude a₀, particularly at integer multiples of a₀, so as to reduce sensitivity to nonlinearities in the phase shift and in the detector" [14].

The algorithm achieves a similar measurement performance as the linear PSI algorithms [14], but for a practical implementation there are some undefined or unmentioned parameters that must be known for the correct use of the algorithm. In the next section we will explain these limitations and their overcoming in detail.

3. Practical implementation

If the phase shift amplitude a can be specified exactly, the phase offset φ is known and also the sampling rate can be adjusted in a certain range to an integer multiple of the oscillation frequency f₀, then the examples parameterized by de Groot [14, 16] but also the methods, described in other articles [6, 8, 9], can be used for the measurement of the height dependent phase Θ.

In many cases the sinusoidal phase shift is performed by simple phase shifters such as electromagnetic coils [13] or piezo-electric transducers [11, 12, 17] without any position feedback or even the possibility to adjust the phase shift amplitude a. Depending on the respective settings of the interferometer the amplitude $\widehat{z}$ may vary and is a priori unknown. The same applies to the phase offset φ. But also if these parameters are known, there are still at least two problems: first, how to define the sampling vectors $h_{\text{odd},\text{even}}$ for an arbitrary phase offset $\phi \ne 0$ and second, how to select the proper weighting factors γ_n for a particular amplitude a₀ to fulfill the design rules. In the following we will propose our preprocessing design steps to obtain the desired parameters.

3.1. Quasi-analytical estimation of the phase offset φ and phase shift amplitude a

An interference signal is usually sampled synchronously to the periodic oscillation of the phase shifter by a detector such as a CCD/CMOS camera or a photodiode at equidistant timing intervals. Especially in systems without any feedback the sampling period may begin at any arbitrary phase position with respect to the beginning of the assumed cosine shaped excitation signal. The original algorithm [14] expects an interference signal, symmetrically sampled to the cosine shaped excitation signal, as shown in Fig. 1, so that Eq. (10) can be used. To take into account an arbitrary phase shift offset φ we introduce new sampling vectors, which calculate a shifted cosine transform for a given harmonic n:

These sampling vectors are used in the same manner as Eq. (10) but not restricted to φ = 0. In comparison with Eq. (10) the normalization values are omitted, since they can be taken into account by the weighting factors γ_n. The transform results are a purely real values with correct sign and amplitude.

The phase offset φ is obtained from the measured interference signal by calculating the argument of the complex-valued Fourier transform for a single frequency bin n

of a first portion of the interference signal $\mathrm{\Delta }\overline{I}\left(t\right)$ (see Eq. (7))

Depending on the selected harmonic n the estimated phase offset φ changes with a multiple of n. To achieve an unambiguous result in the range of $\left(-\pi ,\pi \right]$ one should use the first harmonic or some a priori information by using higher harmonics.

The phase shift amplitude a is estimated in three steps: first, coarse estimation of a, second, selection of “optimal” harmonic pair and finally, fine estimation of a.

The coarse estimation of a uses the DFT amplitude spectrum of the interference signal. As described in [11–13] one can use a digital phase lock-in technique to obtain Θ by evaluation of a part of the interference signal between the turning points of the sinusoidal phase shift signal at fringe frequency f_f. The value of f_f corresponds to the frequency with the highest amplitude of the interference signal and depends on the phase shift amplitude a. To use this technique for the estimation of the phase shift amplitude a a set of simulations in the range $3\le a\le 15$ and $0\le \mathrm{\Theta }<\pi$ was performed. A linear dependence of the frequency f_f on the amplitude a was found, so a simple linear regression model $a\left(f_f\right)=m{f}_f+b$, where m is the slope and b the intercept of the line was estimated, which can be used for coarse estimation of a.

The more exact estimation of the phase shift amplitude a (third step) is based on the ratio of two harmonic amplitudes n and $n+2$. A similar technique to estimate the phase shift amplitude a based on the magnitudes of the first and third harmonic for a limited region $a<3.5$ was already introduced by Sasaki et al. in [6]. The same technique as in [6] was also used in [18] in context of spectral domain optical coherence tomography for the removal of complex-conjugate ambiguity resulting in overlapping mirror images.

The quotient $H_n\left(\mathrm{\Theta }\right)/H_{n+2}\left(\mathrm{\Theta }\right)$ of the frequency transform of the interference signal (Eq. (7)) for two harmonics (n and $n+2$) using new sampling vectors (Eq. 17) and conversion to the a-dependent Bessel functions results in:

There is no analytical solution to determine a₀ from Eq. (19), but one can use a numerical optimization algorithm for nonlinear functions limited to a certain search range (constrained optimization) [19] to minimize the difference of the left and the right side of Eq. (19). A necessary condition for the optimization is that the function $f\left(a\right)={\mathrm{J}}_n\left(a\right)/{\mathrm{J}}_{n+2}\left(a\right)$ is injective in the area to be optimized. To achieve a meaningful value for a a proper selection of the harmonic pair is necessary. This is done in the second step of the estimation procedure.

The selection of the “optimal” harmonic pair is based on the coarse value of a and the strengths of the first $P/2-1$ harmonics of the interference signal, which are depending on the Bessel functions for a particular value a₀ and the current value ofΘ. An algorithm selects a harmonic pair with the highest amplitudes but also ensures that their ratio is far away form the poles, where the function values tend towards $\pm \infty$. The detection of the poles is based on the coarse estimate of a from the first step. The correct estimation of the phase shift amplitude a is only possible if the coarse a-estimation from the first step was successful. This assumes that the estimation of the phase offset φ is also correct. A similar dynamic technique, but with three harmonics, was proposed in [20]. THe authors extend a so-called Pernick-Method [21] to a wide dynamic range up to a phase shift amplitude $a=100\pi$. Their quite complex and not completely disclosed search algorithm achieves a small maximum error of 1% in phase shift amplitude estimation, and doesn’t relay on the value of φ, but operates only with low-distorted signals.

3.2. Simultaneous estimation of a and φ

In most cases, the previously described algorithms yielded largely plausible estimates of the parameters. But the estimation of the phase shift offset φ and the amplitude a may fail, especially if the estimation of φ fails due to a current value of $\sin (\mathrm{\Theta })=0$. For certain phase shift amplitudes a, where the value of the first Bessel function is zero (e.g. $J_1\left(a=3.83\right)$ or $J_1\left(a=7.02\right)$), the same problem of inaccurate φ-estimation will also occur. Also the quality of the coarse a-estimation, due to low SNR or poor performance for small amplitudes ($a<5$), may lead to a incorrect determination of the phase shift amplitude. Therefore, it was important to validate the estimated parameters with respect to their plausibility by comparison with the actual interference signal. The review was largely based on the experience of the user and was not always practical. To overcome this problem a brute force simultaneous estimation of these both parameters was developed.

The approach is simple but effective. The first period (or the first few periods) of the interference signal is compared with simulated signals of all possible combinations of the parameters $\phi ,\mathrm{\Theta }$ and a. The optimization parameter is the standard deviation of the difference between the measured and the simulated interference signal. The simulated signal with the smallest deviation from the measured signal provides the desired search parameters. Figure 2 shows the flow chart of the algorithm.

 

Fig. 2 Algorithm for the simultaneous estimation of the phase shift amplitude a, the phase offset φ and the height dependent phase Θ. A discrete search parameter range (up to 36 values for each parameter) limits the total number of iterations, but still covers all possible variations in the parameters that can occur in a given setup.

Download Full Size | PDF

To perform the parameter search in a reasonable time the range of possible amplitudes a is limited to appropriate values between $a_{min}$ and $a_{max}$ which are separated by a value a_step. The ranges for φ and Θ are separated by ${\phi }_{\text{step}}$ and ${\mathrm{\Theta }}_{\text{step}}$ and vary from 0 to π or $-\pi$ to π, respectively. The measured interference signal I, used for the parameter estimation, is free of the DC component and is normalized to $\pm 1$. For each a the standard deviation for all available φ and Θ values is computed. Then the minimal ${\sigma }_{min}$ inside the 2D array ${\sigma }_{\mathrm{\Theta },\phi }$ together with corresponding values of ${\mathrm{\Theta }}_{min}$ and ${\phi }_{min}$ is determined. If ${\sigma }_{min}$ is smaller than a global value ${\sigma }_{min,\text{global}}$ an update of all search parameters occurs and the loop continues with the next value of a. At the end of the a-loop the best matching values $a_{\text{best}},{\phi }_{\text{best}}$ are available, but they are limited to their discrete steps. To increase the accuracy of a, a 1D parabola approximation to the data of the ${\sigma }_{\text{global}}$-array around the location of ${\sigma }_{min,\text{global}}$ is performed. To obtain a more accurate value of φ, a 2D parabola approximation to the data within the ${\sigma }_{\mathrm{\Theta },\phi ,\text{global}}$-array around the position of ${\sigma }_{min,\text{global}}$ is performed. Additionally, the whole search procedure can be repeated with new search ranges located around the values estimated in the first run of the algorithm.

It should be noted, due to the definition range of the cosine function it is not possible to distinguish between $\cos (a\cos (\alpha +\phi )+\mathrm{\Theta })$ and $\cos (a\cos (\alpha +\phi +\pi )-\mathrm{\Theta })$. Thus, as shown in Fig. 2, the phase offset φ ranges only between 0 to π. Negative phase offset values φ can not be identified unambiguously and lead to an opposite sign of the Θ values obtained later by the evaluation algorithm. In some measurement devices, e.q. as described in section 5.3, the sign can be determined by detection of scan direction. Generally, additional knowledge of the sign is necessary.

3.3. Determination of the weighting factors ${\mathrm{\gamma }}_{\mathrm{n}}$

Proper weighting factors γ_n are essential for the reliable function of this algorithm. As mentioned in section 2 and in [14] filter functions $F_{\text{odd}\left(a\right)}$, $F_{\text{even}\left(a\right)}$ are used to adjust the parameters γ_n according to the design rules described therein. Due to high number of different parameters such as the number of harmonics and variable phase shift amplitude a the determination of γ_n is based on a numerical optimization.

The design rules describe a typical multi-objective optimization problem [22]. A classical approach to solve such kind of problems is to assign a weight to each objective function and sum them up, so that the problem is converted to a single-objective problem [23]. This approach together with the unconstrained optimization based on the BFGS algorithm [24] from the ALGLIB library is used to find the weighting factors γ_n. The algorithm minimizes the objective function $f_{min}\left({\gamma }_n,\mathbf{w},\mathbf{a}\right)$ by using the function gradient calculated (by the library) through numerical differentiation. The objective function is defined as follows:

 

Fig. 3 Variation of the normalized filter functions F_odd/Γ_odd and F_even/Γ_even as a function of the phase shift amplitude a for 10 harmonics and the design amplitude a₀ = 5.175. The function values in the ranges a₁ to a₃ are used for the optimization.

Download Full Size | PDF

In the first range $0\le a_1\le 0.5a_0$ and the third range $1.5a_0\le a_3\le 5a_0$ the squared sum of the filter function values is minimized. For both ranges relatively small weights ($w_1+w_3<0.1$) are used. The second range is the most important, here the difference of the function values around the design amplitude a₀ in the range $0.75a_0\le a_2\le 1.25a_0$ is minimized. In some cases, due to multiple reflections, strong interference signal distortions corresponding to integer multiple of a₀ ($2a_0,3a_0,\mathrm{...}$) may appear [25]. To reduce the sensitivity of the algorithm at corresponding phase shift amplitudes, Eq. (20) can be extended in a similar manner as described before by additional weighting parameters w_n.

Due to the complexity of the objective function with a lot of local minima, the resulting weighting factors may vary depending on the starting parameters. Therefore, it is a good strategy to repeat the optimization with different start parameters to improve the final search result.

If the parameter of the measurement device (sampling frequency, phase shift amplitude) can be adjusted in a certain wide range, then an appropriate parameter selection should follow. Due to DFT based mathematical background of the algorithm, the measurement uncertainty decreases with the increasing number of sampling points per period. In the original algorithm the maximum harmonic is selected to fulfill the sampling theorem with at least two sampling points per harmonic signal period. We propose to use an oversampling factor of two for the highest harmonic used for the evaluation: $n_{max}\le P/4$. It decreases the sensitivity of the algorithm to the random noise. A high phase shift amplitude ($a>10$) is not preferable. Higher amplitudes increase the requirements on the phase shift hardware and introduce unpleasant audible frequencies especially in camera based systems were the maximum phase shift frequency is in a lower single-digit kilohertz range. Another advice is to limit the maximum harmonic where its amplitude is reduced up to 10% of the amplitude of the strongest lower harmonic. Another hint is to exclude harmonics from the evaluation those values of the corresponding Bessel functions for a given amplitude are zero or near zero. This will decrease the sensitivity of the algorithm to noise which may be in the same range as the particular amplitude of the harmonic.

4. Temporal high-speed extension of the evaluation algorithm

The evaluation algorithm described in previous sections and in [14] calculates a single height value for each period of the sinusoidal phase shift signal. The extension of the original algorithm to arbitrary phase shift offset φ allows the evaluation to be started at any position within the excitation signal period. By creating a separate set of sampling vectors for each sample position within a reference signal period, taking into account the relative phase offset to the previous value and repeatedly applying them to the interference signal, it is possible to calculate a height value for each sampled intensity value, as shown in Fig. 4.

 

Fig. 4 Simple implementation of the sliding evaluation algorithm with a set of P different sampling vectors, where each sampling vector is shifted by Δα to the previous one.

Download Full Size | PDF

However, the proposed evaluation strategy is computationally ineffective. Therefore, based on the translation / time shifting property of the Fourier transform, a sliding recursive technique is introduced.

The time shifting property of the Fourier transform results in:

In the case of a discrete sampled finite signal shifted by k points it is expressed by:

Consider the sampling window of the length P as a circular buffer, where by the shift of the window by one sample only a single new value x_m is stored in the buffer while the last value $x_{m-P}$ is deleted, the spectral component ${\underset{\_}{X}}_n\left(m\right)$ at the new sampling window position m will differ to the old spectral component ${\underset{\_}{X}}_n\left(m-1\right)$ at previous position $m-1$ only by the contribution of x_m, $x_{m-P}$ and the additional phase shift. Therefore, one can apply the shifting property of the Fourier transform to the intermediate result $\left({\underset{\_}{X}}_n\left(m-1\right)+x_m-x_{m-P}\right)$ to get the following common computation rule for the new spectral component, also known as sliding DFT [26]:

Looking at the previously described evaluation algorithm, it calculates nothing else than cosine transforms or more precisely the strengths of different frequency bins n, which are the harmonics of the sinusoidal phase shift signal, additionally weighted by the factors γ_n to fulfill the algorithm’s requirements. To apply the sliding DFT technique to the evaluation algorithm it is necessary to consider the initial phase shift φ and also compensate for the phase shift introduced by each movement of the sliding window.

Consider the transform according to Eq. (23) for the signals described by the following equation:

To compensate for a shift by p samples we need to extend Eq. (26) to:

To calculate the shifted DFT as in Eq. (27) by using Eq. (24) we need to compensate the result $X_n\left(m=p-1\right)$ from the previous $p-1$ shift by ${\mathrm{e}}^{\mathrm{i}n\phi }{\mathrm{e}}^{\mathrm{i}2\pi n\left(p-1\right)/P}$ and then translate the sum by ${\mathrm{e}}^{-\mathrm{i}n\phi }{\mathrm{e}}^{-\mathrm{i}2\pi np/P}$:

After simplification of Eq. (28) we get:

Due the fact that the signal $x\left(j\right)$ being transformed at the window position $m=p$ has a current phase shift of $n\left(\phi +2\pi p/P\right)$ (compare with Eq. (27)) and the DFT value is purely real, the imaginary part of the Eq. 29 can be omitted and we get the following result:

Now, Eq. (30) can be used to implement the sliding evaluation of the sinusoidally phase shifted interference signals, which we call sliding SinPSI. The evaluation begins just as the previously described algorithm with the estimation of parameters and the generation of the sampling vectors. The first Θ-value is calculated by a formula similar to Eq. (8) by using the new sampling vectors $\tilde{h}\left(n,j\right)$ (see Eq. (17)):

The following Θ-values are calculated by the formula based on the sliding real DFT (Eq. (30)) and on Eq. (31):

The sliding SinPSI evaluation algorithm is computationally very efficient. The sliding sampling vectors $S_{\text{odd},\text{even}}\left(p\right)$ can be precalculated prior to the evaluation during the parameter estimation. Therefore, apart from the computation of the first Θ-value and the arctan-computation, each new height value requires only two real multiplications and four additions.

5. Results

Since the basic structure of the algorithm has not been changed and the height values are calculated in exactly the same way as by the original algorithm by P. de Groot, the formulas for the error analysis still apply [14, 16, 27]. An important requirement for a proper work of this algorithm is a synchronous detector triggering to the sinusoidal phase shift signal as well as integer multiple and constant amount of sampling points per phase shift period. It can be achieved if the sinusoidal phase shift and the detector trigger signals have a common master clock source [12, 13].

Here we will show the influence of the erroneously estimated or assumed parameters on the evaluation result in some worst case scenarios as well as the need of the temporal high-resolution extension of the algorithm.

5.1. Influence of phase offset φ on height determination

First of all, we want to emphasize the need for a reliable determination of the phase offset. As long as the phase offset can be determined with an accuracy of approximately 15 deg, the resulting height error is quite small. An overview of the dependence of height measurement errors on small phase shift errors (${\phi }_{\text{err}<}$ 16 deg) is given in [27]. In case of an algorithmic determination of this parameter and possible problems mentioned in section 3.1, the error in phase offset may cover all possible values in the range of $\pm \pi$. Thus, the measurement error may be also very big. It depends on the number of selected harmonics, shows local maxima at certain error offsets ${\phi }_{\text{err}}$ and increases rapidly till infinity if the phase offset error approaches towards $\pm \pi /2$. Figure 5 summarizes the analytical maximal error in height determination given by the Eq. (162) defined in [16] and shows exemplary four evaluation results at different phase offsets.

 

Fig. 5 (a) Theoretical error in height determination depending on the phase offset error φ_err, given by Eq. (162) in [16]. (b) Height deviation from linear ramp with 2 nm slope per evaluation period for different phase offsets φ_err. Interference signal generated with a = 5, φ = 0, P = 50, n_max = 7 and λ = 850 nm. Maximal deviations: Δh_φ=0 = 0.019 nm; Δh_φ=π/8 = 28.8 nm; Δh_φ=π/4 = 3.85 nm; Δh_φ=π/2 = 104.5 nm.

Download Full Size | PDF

The error in height evaluation is $\lambda /4$-periodic. As expected from Fig. 5(a) the error for ${\phi }_{\text{err}}=\pi /4$ is smaller than for ${\phi }_{\text{err}}=\pi /8$. In the case of ${\phi }_{\text{err}}=\pi /2$ the result has approximately a sawtooth shape, and due to the ambiguity of the phase evaluation the maximum error is limited to values of $\pm \lambda /4$.

It should be noted, that without the use of the new sampling vectors given by Eq. (17), the phase shift offset compensation can only achieved by the adjustment of the sampling interval $\mathrm{\Delta }\alpha$ and relative offset α₀ to the phase shift excitation signal (see Fig. 1) or the change in the phase shift frequency f₀. The introduction of the new sampling vectors relaxes the requirement on the hardware and enables an easier use of this evaluation algorithm.

5.2. Robustness of parameter estimation

As mentioned in section 3.2, the estimation of parameters a₀ and φ may fail, if the current value $\sin (\mathrm{\Theta })\approx 0$ or the values of the Bessel functions for certain amplitudes a are zero. In addition, the signal quality, i.e. the signal-to-noise ratio (SNR) influences the parameter estimation. To compare the techniques described in section 3.1 with the brute force technique described in section 3.2, a set of Monte-Carlo simulations was performed. Interference signals were calculated by Eq. (2), whereby the amplitude a was randomly varied in the range from 3 to 15 and the initial phase offset φ as well as the height dependent phase Θ were varied in the range of $\pm \pi$. The number of sampling points per period was set to $P=50$ and two phase shift periods were used. The noise was given by the SNR value, equally distributed in the range from 0 to 100 dB, which is defined as follows [28]:

To achieve a reliable coverage of all possible parameters, 11 million cycles were performed. Figure 6 shows the results of these simulations.

 

Fig. 6 Comparison of parameter estimation algorithms depending on the signal quality. (a) and (b) standard deviation and maximal error in phase estimation φ. (c) and (d) standard deviation and maximal error in estimation of phase shift amplitude a. "QA": quasi-analytical, "BF 1./(2.) run": brute force 1./(2.) run, "BF 2. run + fit": brute force 2. run with additional LSQ fit.

Download Full Size | PDF

Due to $\sin (\mathrm{\Theta })=0$ values, even nearly ideal interference signals are incorrectly evaluated by the quasi-analytical algorithm, which leads to the maximal error of 90° in the phase offset estimation and consequently in an incorrect estimation of the phase shift amplitude a, as can be seen in Figs. 6(b) and 6(d). The parameter estimation with the brute force algorithm may fail only for very noisy signals with a SNR below 5 dB. In practice, interference signals show SNR values in the range between 10 and 60 dB. For SNR $\ge$ 10 dB the rms error in height determination using this algorithm is $\sigma \le 2.65$ nm. Therefore, in the worst case the error in phase offset determination will be less than 3° and the error in amplitude estimation below 0.4. As described in section 3.2, it is appropriate to repeat the brute force parameter search after the first run with new search ranges located around the values found there. For simulations shown in Fig. 6 $a_{\text{step}}=0.25$ and ${\phi }_{\text{step}}=\pi /30$ were used in the first run and were reduced by one order of magnitude for the phase shift amplitude and by factor of 16.45 for the phase offset estimation in the second run. Additionally, least squares fits (1D and 2D parabola fits) were used in the second run to increase the accuracy of the estimation.

The maximum error in brute force parameter estimation for higher SNR values mainly depends on small phase shift amplitudes $a<5$. The simple normalization technique of the interference signal, as described in section 3.2, removes the DC component from the signal which contains a relatively large amount of the signal energy (see Eqs. (1) and (5)). Thus, the brute force technique, which takes place in the time domain, lacks on information inside of the data due to decreased SNR.

5.3. Application of sliding SinPSI

An application of sinusoidal phase shift interferometry is the distance measurement between interferometer and specimen during the depth-scan in coherence scanning interferometry (CSI) [29, 30]. The distance information is used for reconstruction of CSI signals distorted by mechanical vibrations [13, 17, 31]. With the knowledge of the real distance for each camera image acquisition, it is possible to rearrange captured images according to their measured positions and thereby compensate for the influence of vibrations. The acquisition / integration time of each image takes only a few microseconds. Thus, the corresponding distance values should be acquired at high data rates to perform a proper assignment to the integration interval. In the interferometer described in [17] an ultrasonic piezo transducer is used to perform sinusoidal phase shifts at the frequency of approximately 40 kHz. In some situations, at high vibration frequencies and high amplitudes this distance data rate is not high enough to calculate a continuous distance course without phase jumps. By applying the sliding SinPSI evaluation algorithm it is possible to calculate a continuous distance course and therefore also to compensate for higher frequency / larger amplitude vibrations in CSI. Figure 7 shows evaluation results obtained with conventional and sliding SinPSI algorithms from the same dataset.

 

Fig. 7 Distance measurement during a CSI depth-scan without vibrations and in presence of high aperiodic vibration caused by a hammer stroke. Several phase jumps (shown in lower inset) appear at the time 425 ms in the result obtained by the conventional SinPSI algorithm. Upper inset shows the height data density of the sliding SinPSI result. The SNR of the interference signal is 27 dB. For better comparability, an offset of 1 μm was added to the blue curve.

Download Full Size | PDF

Due to fast distance changes (more than $\pm \lambda /4$ between two consecutive sinusoidal excitation periods), caused by a hammer stroke, several phase jumps appear in the distance course obtained by the conventional SinPSI algorithm. After the vibrationhas receded the remaining distance error is approximately 2 μm (deviation between red and greed curve on the right-hand side). The much higher distance data rate (104 samples per phase shift period) of the sliding SinPSI algorithm allows to obtainthe distance course without any phase jump and therefore can be used for a successful vibration compensation.

The large number of distance data points, generated by the sliding SinPSI algorithm was used in another CSI instrument with lower sinusoidal phase shift frequency (2 kHz) also for vibration compensation. As discussed in [13] fast distance changes compared to the sinusoidal phase shift frequency lead to sparse distance courses with quite big distance changes between consecutive height values. This results in erroneous assignment of the camera images with respect to their acquisition position. By increasing the distance data rate by a factor of 200 with the sliding SinPSI algorithm it was possible to compensate for periodic 50 Hz vibration with 1.12 μm amplitude. The SNR of the interference signal was 44 dB. In this case the distance changes between consecutive height values were reduced from 185 nm, obtained by the conventional SinPSI algorithm, to less than 1 nm, by using the sliding SinPSI technique.

5.4. Error analysis in sliding SinPSI evaluation

A comparison of height profiles obtained by sliding and conventional SinPSI algorithms for measurements with fast distance changes revealed periodic height errors in the sliding evaluation algorithm resulting between points, which are equal in both algorithms. The errors are greater, the faster the height change occurs. To investigate this behavior a set of simulations with different height change rates $\mathrm{\Delta }h$ per sampling point for interference signals of different initial phase offsets φ were performed. The height change was either 1 nm or 0.1 nm per sample point, the phase shifting amplitude was set to a = 5, and the initial phase offset φ was changed from 0 to π in steps of $\pi /4$. To ensure a correct estimation of a and $\check{a}rphi$ the height was not changed during the first phase shift period of the length $P=100$. The obtained results are shown in Fig. 8.

 

Fig. 8 Height errors in sliding SinPSI evaluation in dependance on height change Δh per sampling point and initial phase offset φ. a). Δh = 1 nm; b). Δh = 0.1 nm. Curves with dot markers represent evaluation with conventional SinPSI algorithm, solid curves show sliding SinPSI evaluation results. For better comparability additional height offsets (100, 200 and 300 nm respectively 10, 20 and 30 nm) were added to the profiles with φ ≠ 0.

Download Full Size | PDF

The maximum height error is proportional to the distance change rate $\mathrm{\Delta }h$ and does not depend on the initial phase offset. Height values obtained with the conventional algorithm match exactly with the values obtained with the sliding algorithm at the integer multiple positions of P. If we pick a height value at arbitrary position from the sliding SinPSI result and then investigate values at integer multiple of P, the picked sparse distance course has the same expected distance change. The observed behavior is most likely due to the continuous change in the amplitude of the harmonics (amplitude modulation) of the interference signal and may probably be reduced due to its systematic behavior.

The sliding SinPSI technique is based on the sliding DFT [26]. Due to its recursive implementation, each error in the sampled interference signal (e.g. phase shift error or noise) directly propagates into the result, there is no averaging. In the original implementation the result is an average of all data from the whole phase shift period, but not in the sliding algorithm. Thus, especially deviations (non-linearities) in the sinusoidal phase shift excitation signal cannot be canceled and need to be measured or estimated. It should be noted, that these errors are quite small (less than one nanometer) and only observable because of a much higher position (height) data rate compared with the original algorithm, which provides only a single (averaged) height value per period. Nevertheless, high amount of additional height values (even slightly corrupted) in between the ones measured by conventional implementation, allow to avoid phase jumps, which may appear by fast distance changes and completely distort the measurements.

6. Conclusion and outlook

An extended version of the evaluation algorithm for sinusoidally phase shifted interference signals based on the original de Groot algorithm [14] was developed. It was discussed, how to overcome the limitations, not handled or not considered by the original algorithm. Two different techniques to estimate the phase shift amplitude a and the offset φ from an unknown interference signal were developed. Both parameters need to be known for practical implementation and parametrization of the algorithm. A numerical optimization technique for the calculation of weighting factors γ_n, which are essential for the robustness of the algorithm, based on the design rules proposed by de Groot was introduced. Furthermore, a temporal high-speed technique, which extends the conventional algorithm and allows to obtain a height value for each sampling point of the interference signal was developed. Based on simulations and practical measurements the performance of the proposed parameter estimation and the evaluation of interference signals by the sliding SinPSI algorithm was demonstrated.

The sliding SinPSI technique opens new possibilities for the use of existing measurement devices without any hardware changes. Now it can be used even with slow camera systems (lower 100 fps) for the real-time phase evaluation at the camera rate.

Periodic errors in height evaluation with the sliding SinPSI algorithm, observed by continuous distance changes during signal acquisition, as shown in section 5.4, require further investigation. One possible reason for such deviations could be the amplitude modulation of each harmonic component during the fast height change. Although these deviations can be kept within a reasonable range they can be further examined by Phase-Locked-Loop (PLL) techniques known from power grid control systems [32].