Integration in the Fourier domain for restoration of a function from its slope : comparison of four methods

In some measurement techniques the profile, f x , of a function should be obtained from the data on measured slope f 0 x by integration. The slope is measured in a given set of points, and from these data we should obtain the profile with the highest possible accuracy. Most frequently, the integration is carried out by numerical integration methods [Press et al., Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1987)] that assume different kinds of polynomial approximation of data between sampling points. We propose the integration of the function in the Fourier domain, by which the most-accurate interpolation is automatically carried out. Analysis of the integration methods in the Fourier domain permits us to easily study and compare the methods’ behavior. © 2002 Optical Society of America OCIS codes: 120.0120, 120.6650, 200.3050.

In this Letter we address the problem of evaluation of the prof ile, f ͑x͒, of a function from the data on measured slope, f 0 ͑x͒, by integration.The slope is measured in a given set of points x n nDx, n 0, . . ., N 2 1 [to shorten the notation we write f 0 n instead of f 0 ͑x n ͒], and from these data we should obtain the prof ile with the highest possible accuracy. 1Many numerical integration methods have been described in the literature. 2 The best-known numerical integration methods are based on the Newton-Cotes quadrature expressions such as the trapezoidal rule, the Simpson rule, and the 3͞8 Simpson rule.In general, these methods are designed to evaluate an integral in a given interval by some iterative procedure.The integral is evaluated for a given sampling distance Dx of the slope function, f 0 ͑x͒, then this distance is decreased (i.e., Dx͞2), and the integral is evaluated again.The iterations are continued until some criterion of convergence is fulfilled.The higher the order of the rule, the faster the convergence.For the iterative procedure to be applied, the slope function, f 0 ͑x͒, has to be evaluated analytically or numerically at the sampling points.However, for the problem of prof ile measurement that we are dealing with, this cannot be done.The slope function is available only at the measured points, and no more measurements can be performed.Then, the question is which integration method is the most appropriate.In Ref. 3, Elster and Weingärtner tested several integration methods with different functions to determine the best method in terms of reconstruction errors.The best was cubic spline interpolation of the slope data, followed by subsequent integration of the spline function.
In general, all integration methods, such as Newton-Cotes quadrature expressions and spline interpolation, assume some kind of interpolation be-tween the data.It is known that, given sampled data, the best interpolation is sinc interpolation, 4 which has an efficient computational implementation in the Fourier domain.Therefore, here we propose integration of the function in the Fourier domain when the most-accurate interpolation is automatically carried out.The facts that integrations methods such as the Newton-Cotes quadrature expressions are recursive and that they are particular implementations of signal digital filtering permit us to study all these methods in the Fourier domain.We propose this kind of analysis as a tool for comparing the different integration methods and how the noise will affect the final result.
First, we describe the first three Newton -Cotes quadrature rules and the integration in the Fourier domain.Then, we show that all these methods can be studied as some kind of digital filtering in the Fourier domain, and the frequency responses of the integration methods are deduced.These frequency responses or transfer functions give us an adequate tool for analyzing the different integration methods.
The three first Newton-Cotes quadrature rules are the trapezoidal, the Simpson, and the 3͞8 Simpson rules. 5The closed trapezoidal, Simpson, and 3͞8 Simpson rules are given by To obtain the prof ile f i we need to apply these expressions recursively.In all the methods we arbitrarily choose the value of the integral in the first point to be zero.Then, for the trapezoidal rule, we obtain In the Simpson rule the second point needs to be evaluated by the trapezoidal rule; then, In the 3͞8 Simpson rule the second and the third points need to be evaluated by the trapezoidal rule and the Simpson rule, respectively; then, In the cubic spline interpolation, a cubic polynomial is evaluated between every couple of points 3 and then an analytical integration of these polynomials is performed.
As mentioned above, the most-accurate signal interpolation is provided by discrete sinc interpolation. 4 Discrete-signal sinc interpolation can be easily implemented in the domain of a discrete Fourier transform (DFT). 6This suggests the use of DFT-based signal integration.Let ͕F 0 ͑r͖͒ be samples of the DFT spectrum of slope function ͕ f 0 k ͖: Then samples of the function can be obtained as an inverse DFT: where ͕h r ͖ are samples of frequency response of the integration f ilter: ‫ء‬ is a phase conjugate and N is the number of the samples.In this case we assume an even number.
The zero frequency must be real.According to the fourth line of Eq. ( 6), h r h ‫ء‬ N 2r , for r N͞2 the frequency response must be real.There are several possibilities for the value of h N /2 : It can be made equal to zero, but then the information of this frequency is lost; it can be made equal to 21͞p, following the second line of Eq. ( 6); or it can be made equal to 21͞2p.It can be shown that the oscillations of the impulse response of the frequency response [the point-spread function (PSF)] are lower with the last choice, so we adopted this value.
All integration methods can be treated as signal convolution with an integration convolution kernel (PSF).Therefore, the simplest way to compare the integration methods is to compare their PSFs' and (or) their corresponding frequency responses.The PSFs' of the trapezoidal, Simpson, and 3͞8 Simpson integration methods are specified in Eqs. ( 1)-( 3), from which one can see that they are recursive digital filters.Then, their analysis is much easier in the Fourier domain. 4Therefore we compare the methods in terms of their frequency responses for trapezoidal, Simpson, and 3͞8 Simpson, integrators.Their responses can be found from Eqs. ( 1) -( 3) as follows: By DFT of these expressions we obtain where we denote with capital letters the DFTs of the functions in lowercase letters: Therefore the frequency responses of these integrators are, respectively, The knowledge of frequency responses ͕h r ͖ of the integrators also allows us to easily evaluate how they transform additive signal-independent input noise.If the spectral density of the input noise is N inp r , the spectral density of the integrator output noise is If input noise is uncorrelated, its spectral density is a constant proportional to the noise variance, s 2 inp .Variance of the output noise in this case can be computed as In Fig. 1 the absolute value of the frequency responses of the Fourier method of integration [Eq.( 6)] and of the Newton-Cotes rules are represented.Because the absolute value of the frequency response is  symmetric we represent only half this curve.The maximum frequency is normalized to 0.5.From the figure one can see that the integration methods are low-pass filters and that they are almost identical at low frequencies, whereas at higher frequencies the 3͞8 Simpson, Simpson, and trapezoidal integration methods tend to introduce integration errors.The frequency response of the 3͞8 Simpson rule tends to infinity for the 2͞3 of the maximum frequency, and the frequency response of the Simpson rule has almost the same tendency for the maximum frequency.These tendencies mean that the noise in these frequencies will be enhanced, as we show in the next example.In Fig. 2 an example of reconstruction of the prof ile from the slope is shown.Figure 2(a) shows the original slope and the slope with additive Gaussian noise of zero mean and variance 1.In Fig. 2(b) the original prof ile and the prof iles obtained with the four integration methods are shown.As the frequency of the profile is not very high, the reconstructions obtained with the three methods are very similar.Figure 2(c) shows a magnif ied part of Fig. 2(b) in which only the prof iles obtained with the Fourier method, the trapezoidal rule, and the Simpson rule are drawn.The trapezoidal rule and the integration in the Fourier domain give almost the same result for this function, whereas the Simpson rule presents high-frequency oscillating behavior (the period is two pixels) because of the frequency response of the high frequency.Figure 2(d) is similar to Fig. 2(c); in the latter case the prof ile obtained with the 3͞8 Simpson rule is drawn.Here oscillating behavior that is due to the frequency response of the integration method (the period is three pixels) is also visible.As can be seen from Fig. 1, this frequency response tends to infinity for a frequency that is 2͞3 of the maximum frequency.
Then, from this analysis we can conclude that known numerical integration methods are not well adapted for these kinds of problems.For the function that we have used, the trapezoidal rule and the Fourier domain method give similar results, although, as follows from the analysis of their frequency responses, the trapezoidal rule tends to excessively attenuate signal high frequencies.We have also shown that the study of the integration methods in the Fourier domain by means of their frequency response is a powerful tool for studying the behavior of these methods in solving the problem of obtaining the prof ile from a given set of measured slopes.
This work is dedicated to Domingo Gonzales of the University of Zaragoza, on the occasion of his retirement.This work was partially financed by Ministerio de Ciencia y Tecnología (Spain) under project BFM2000-0036-C02-01 and by European Community project CTB556-01-4175.J. Campos's e-mail address is juan.campos@uab.es.

Fig. 1 .
Fig. 1.Comparison of the frequency responses of the Fourier, trapezoidal, Simpson's and 3͞8 Simpson methods of integration.

Fig. 2 .
Fig. 2. (a) Slope of the test function.Bold curve, original function; thin curve, with added noise.(b) Original and reconstructed prof iles obtained with the four methods.In (a) and (b) the curves essentially overlap (c), (d) [enlarged portions of (b): (c) prof iles obtained with the Fourier method and the Simpson rule, (d) prof iles obtained with the Fourier method and the 3͞8 Simpson rule; the oscillating curves correspond to the Simpson and 3͞8 Simpson methods.