Frequency responses and resolving power of numerical integration of sampled data

Methods of numerical integration of sampled data are compared in terms of their frequency responses and resolving power. Compared, theoretically and by numerical experiments, are trapezoidal, Simpson, Simpson-3/8 methods, method based on cubic spline data interpolation and Discrete Fourier Transform (DFT) based method. Boundary effects associated with DFTbased and spline-based methods are investigated and an improved Discrete Cosine Transform based method is suggested and shown to be superior to all other methods both in terms of approximation to the ideal continuous integrator and of the level of the boundary effects. ©2005 Optical Society of America OCIS Codes: (070.2590) Fourier transforms, (120.6650) Surface measurements, (200.3050) Information processing, (240.6700) Surfaces References and links 1. W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. 70, 9981006 (1980) 2. D.L. Fried, “Least-Square Fitting a Wave-Front Distortion Estimate to an Array of Phase-Difference Measurements,” J. Opt. Soc. Am., 67, 370-375 (1977) 3. R. H. Hudgin, “Wave-Front Reconstruction for Compensated Imaging,” J. Opt. Soc Am. 67, 375-378 (1977) 4. R. J. Noll, “Phase Estimates from Slope-Type Wave-Front Slope Measurements,” J. Opt. Soc. Am. 68, 139140 (1978) 5. A. Dubra, C. Paterson, and C. Dainty, "Wavefront reconstruction from shear phase maps using the discrete Fourier transform," Appl. Opt. 43, 1108-1113 (2004) 6. F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” Appl. Opt. 30, 1325-1327 (1991) 7. W. Zou and Z. Zhang, “Generalized wave-front reconstruction algorithm applied in a Sach-Hartmann test,” Appl. Opt. 39, 250-268 (2000) 8. S. Krey, W.D. van Amstel, J. Campos, A. Moreno, E.J. Lous, “A fast optical scanning deflectometer for measuring the topography of large silicon wafers,” Proc. SPIE, 5523, (SPIE’s 49 Annual Meeting, Denver), (2004) 9. J. Campos, L. P. Yaroslavsky, A. Moreno, M.J. Yzuel, “Integration in the Fourier Domain for restoration of a function from its slope: Comparison of four methods,” Opt. Lett. 27, 1986-1988, (2002) 10. A. Moreno, J. Campos, L.P. Yaroslavsky, “Frequency response of five integration methods to obtain the profile from its slope,” Optical Engineering, (2005) (to be published) 11. L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004) 12. J. H. Mathew and K. D. Fink, Numerical Methods using MATLAB (Prentice-Hall, Englewood Cliffs, N. J., 1999) 13. C. Elster, I. Weingärtner, “High-accuracy reconstruction of a function f(x) when only df(x)/dx is known at discrete measurements points,” Proc. SPIE 4782, (2002) (C) 2005 OSA 18 April 2005 / Vol. 13, No. 8 / OPTICS EXPRESS 2892 #6362 $15.00 US Received 19 January 2005; revised 30 March 2005; accepted 31 March 2005


Introduction
Numerical integration has numerous applications in several optical fields as the wave-front reconstruction from wave-front slope measurements [1][2][3][4].It has found important optical applications such as, for example, shearing interferometry.Recently, a new wave-front reconstruction method based on the Discrete Fourier transform useful for this interferometric technique is presented in [5].In [6] Roddier et.al presented a novel wave-front reconstruction algorithm using an iterative Fourier transforms.Many of these algorithms are applied in the Sachk-Hartmann test as, for example, the generalized algorithm presented in [7].Yet another application of the numerical integration is the optical surfaces determination by laser deflectometry.This technique is based on the measure of the deviation suffered by the incident light in a test surface [8].This deviation contains the slope data information of the profile of the test surface.For this purpose, we begin investigating a one-dimensional integration operation [9,10].
While carrying out numerical computation with sampled data one should realize that sampled data represent, with certain accuracy, data that are originally continuous functions and those numerical algorithms approximate certain continuous transformation of those functions.Consequently, results of the computation should be evaluated with respect to that continuous transformation.In other words, given computational algorithm applied to sampled data, one should find out to what continuous transformation of continuous functions that are represented by the sampled data this algorithm corresponds.In this paper, we address this problem for the case of numerical integration of functions using sampled representation of their derivatives.
Integration of functions can be regarded as a convolution of the functions with a corresponding integration kernel, or point spread function.Different numerical integration algorithms will correspond to different approximations of the ideal integration point spread function.In particular, they will differ in terms of their resolving power, that is, of their capability to resolve between two close sharp peaks in the functions.
Thanks to the convolution theorem for Fourier Transform, convolution integral can be treated in Fourier transform domain as a product of Fourier spectrum of the function and of that of the convolution kernel called convolution frequency transfer function, or frequency response.On account of that, one can also characterize numerical integration algorithms in terms of the accuracy of approximating the ideal integration frequency response.
The purpose of the paper is to investigate, theoretically and by numerical simulation, frequency responses and resolving power of several numerical integration algorithms.Specifically, we investigate trapezoidal integration formula, two modifications of Simpson formula, integration using cubic splines and two integration methods by '1/f-filtering' in the domain of Discrete Fourier Transform that, in a certain sense, approximate frequency response of the ideal integration most closely.We will also show limitations imposed to the integration accuracy by the finiteness of the number of available function samples and by methods of treating boundary effects in the numerical integration.
In Sect.2 we present and compare analytical formulas for frequency responses of continuous integration and its different numerical approximations, including the method of integration in the domain of Discrete Fourier Transform (DFT) using Fast Fourier Transform algorithm.The latter being the best approximation to the ideal continuous integrator in terms of its frequency response, heavily suffers from boundary effects due to finiteness of the number of signal samples.Therefore in Sect. 3 we introduce an improved modification of this method that works in the domain of Discrete Cosine Transform and is much less vulnerable to the boundary effects.The analytical treatment is supported by experimental results of comparison of the methods presented in Sect. 4. In Sect.4.1 we compare the integration accuracy provided by different numerical algorithms for integration of sinusoidal functions with integer number of periods in sampled data, the case of periodic signals, when no boundary effects are observed for the DFT based method.In Sect.4.2 we investigate boundary effects for aperiodic signals of DFT and DCT based numerical integration algorithms and show that the latter provides results with substantially lower boundary effects.In Sect.4.3 we present results on numerical comparison of resolving power of different integration algorithms.In Conclusion, we summarize the results.The paper contains also two appendices.As we compare numerical integration methods and ideal continuous integrator in terms of their frequency responses, we prove, in Appendix 1, that discrete frequency responses of the numerical integrators are samples of frequency responses of the corresponding continuous integration filters.In Appendix 2 we show how DCT based integration method can be efficiently implemented using fast Discrete Cosine Transform algorithms.

Continuous and numerical integrators and their frequency responses
Digital signal integrating is an operation that assumes interpolation of sampled data.Similarly to signal differentiating, signal integrating operates with infinitesimal signal increments.It can also be treated as signal convolution.In the Fourier transform domain, it is described as where ( ) ( ) is the integrating filter frequency response.
In digital processing, integrating filtering described by Eq. ( 1) can be implemented in the domain of Discrete Fourier Transform (DFT) as: where is a set of N samples of the input signal to be integrated, {} .

{} .
IDFT are operators of direct and inverse Discrete Fourier Transforms, sign • denotes element-wise product of vectors, and for even N and for odd N [11], where h is the signal sampling interval.Note that digital signal integrating according to Eq. (3) automatically implies signal discrete sinc-interpolation [11].We will refer to this integration method as DFT-based method.
A number of numerical integration methods have been described in the literature.Most known numerical integration methods are the Newton-Cotes quadrature rules ( [11], [12] and [13]).The three first rules are the trapezoidal, the Simpson and the 3/8 Simpson ones.In all the methods, the value of the integral in the first point is not defined because it affects to the result constant bias and should be arbitrarily chosen.When it is chosen to be equal to zero, the trapezoidal rule is: (5) 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 these integration methods, a linear, a quadratic, and a cubic interpolation, respectively, are assumed between the sampled slope data.In the cubic spline interpolation, a cubic polynomial is evaluated between every couple of points [13], and then, an analytical integration of these polynomials is performed.
As it is shown in Appendix, discrete frequency response of the digital filter (Discrete Fourier Transform coefficients of its point spread function) are samples of the equivalent continuous filter frequency response.Given signal sampling and reconstruction devices, frequency response of the continuous filter that corresponds to a digital filter is fully determined by the discrete frequency response of the digital filter.Therefore, we will compare the numerical integration methods in terms of their discrete frequency responses.Discrete point spread functions for trapezoidal, Simpson and 3/8 Simpson's integrators are determined from Eqs. ( 5)- (7) as follows: (10) Discrete point spread function for the method based on cubic spline interpolation of the input data information can be found as follows: (11) where the coefficients are determined by the system of linear equations: ( ) ( ) Taking N point Discrete Fourier Transform of these expressions (Eqs.( 8)-( 12)) we obtain Therefore the frequency responses of these integrators are, respectively: In Fig. 1, absolute values of the frequency responses of the DFT based method of integration (Eq.( 4)) and of the Newton-Cotes rules and cubic spline method (Eqs.( 17)-( 20)) are represented with a frequency coordinate normalized to its maximal value.Because the absolute values of discrete frequency responses are symmetric, only half of the curves are shown.DFT method frequency response coefficients are, by the definition, samples of the ideal "1/f" frequency response of the continuous integrator.Therefore, with respect to approximation of the ideal frequency response, it can be regarded as a "gold standard".As it can be appreciated in the figure, frequency responses of all methods are similar in the low frequency region, and the Newton-Cotes rules and cubic spline method began to deviate from the ideal one in the medium and high frequency zones.The Simpson and 3/8 Simpson rules exhibit large deviations and even poles (the frequency response tends to infinity) at the highest frequency and at 2/3 of the maximal frequency, respectively.The cubic spline based method is the closest to the "gold standard" DFT-method.

DCT based integrator
Although DFT-based integration (FI) method is the closest approximation to the continuous integrator, this method suffers from boundary effects since it implements cyclic convolution rather than shift-invariant convolution.Boundary effects exhibit themselves in form of oscillations around signal discontinuities that may occur between samples at the beginning and the end of the available signal realization.One can substantially decrease influence of boundary effects by means of signal extension to double length with its "mirror reflected" copy according to the equation: Such an extended signal { } k a , by the definition, has no discontinuities at its ends and in the middle and can be used in the DFT based integration method instead of the initial signal.Frequency response of the integration filter in this case is defined by Eq. (4a), in which the number of samples N should be replaced by N 2 .The value of ( ) Note, that in this case of the doubled signal length the degree of approximation of the ideal "1/f " frequency response of the continuous integrator is even better than that for the above DFT based method because the doubling of the number of signal samples results in two times more dense grid of samples of the frequency response.
The doubling of the number of signal samples in this implementation of the DFT based method does not necessarily doubles the computational complexity of computation.As it is shown in Appendix 2, 2N-DFT convolution of signals obtained by mirror reflection extension can be carried out using fast algorithms of Discrete Cosine Transform and of associated with it Discrete Cosine/Sine Transform for signals of N samples.We will refer to this modified DFT based integration method as "extended", or DCT-based method.

Experimental comparison
In this section we present results of experimental comparison of the methods.First, we will compare accuracy provided by different numerical integration algorithms for integrating sinusoidal signals with integer number of periods in sampled data that, due to the periodicity do not present any discontinuity at signal borders and therefore are free of boundary effects.Then we investigate boundary effects for two modifications of the DFT-method, regular one that uses N-point DFT, and the DCT based method that assumes even extension of the input data to double length and show that the latter substantially reduces integration boundary effects.Finally, we provide results on numerical comparison of resolving power of different integration algorithms.

Integration of periodic sinusoidal signals
In these experiments, sinusoidal signals and their derivatives are generated as, correspondingly: where N is the number of signal samples, p is the frequency parameter of the sinusoidal signal, the normalized frequency is given by 2p/N; randphase corresponds to an initial random phase and x represents the domain of the signal ( ).The frequencies are selected to have an integer number of periods in the signal length.Different integration methods are applied to the derivatives.Then, the obtained functions were compared with the analytically generated signal by estimating the integration root mean square error as is the obtained function with the numerical integration and ( ) is the analytical function of Eq. ( 22).In the experiments, pseudo-random phase data were generated using pseudo-random number generator and the error was averaged over 1000 realizations.In figures 2(a) and 2(b) we represent the average integration error for each of studied integration methods as a function of

N p 2
, the normalized frequency of the sinusoidal signal.In these experiments, signals were evaluated in N = 256 samples.From figures 2, one can see that all methods give similarly low error for low frequencies.When the frequency of the sinusoidal signal is increased, the different integration methods exhibit different behaviour; for example, for a frequency equals to 2/3 of the maximum frequency, the 3/8 Simpson method gives very high error.A similar behaviour occurs for the Simpson method when the frequency is near to the signal maximal frequency defined by the sampling rate.

Aperiodic signals and boundary effects
In order to study the boundary effects for the best methods (DFT-based method (FI), DCTbased method (Extended method) and method based in the interpolation of the slope data by cubic splines (CSI)) another numerical experiment was carried out with sinusoidal signals of a non -integer number of periods.The number of periods is given by the parameter p in Eq. (  Figure 3(a) shows obtained experimental results for the error for the three methods and the signal normalized frequency ν = 0.273 (corresponding to relatively low frequencies).The number of signal samples N was 256 and the number of divisions of the frequency interval n was 20.From the figure one can see that the boundary effects are more severe for FI method than for the CSI method.DCT-based (Extended) method shows errors that are very close though slightly larger then those for the cubic spline method.In the figure we can appreciate that the boundary effects practically disappear after 10-th signal sample.
The sample-wise integration error obtained for signal frequency increased to ν = 0.547 (medium frequency) is shown in Fig. 3(b).In this case the boundary effect error for DFTbased and spline methods are similar while the error for DCT-based (Extended) method is substantially lower.These boundary effects also last only approximately first 10 samples.By comparing these errors with those for the low frequency region we can see that the boundary errors in the first 10 samples increased for all methods and that for DCT-based (Extended) method they are the lowest.In the stationary region (beyond the first 10 samples), the error for the CSI method turned to be higher than that for both DFT-based (FI) and DCT-based methods and is in agreement with figures 1 and 2.
For higher initial high frequency, ν = 0.820, the error produced by the three studied methods is shown in Fig. 3(c).From the figure, we see that the errors produced by the CSI method are much higher than those obtained with both DFT and DCT-based methods.In the CSI method, the stationary errors predominate over those due to boundary effects.We can also see that, for the DFT-based (FI) and DCT-based methods, boundary effects last same 10 first samples and that their values are higher than in the low and medium frequency regions.Finally, in order to check, how boundary effects depend on the number N of samples where the signal and its derivative are evaluated, we repeat the previous numerical experiment for different signal lengths.Figure 4 shows the error obtained in the first 10 samples for the number of the samples N = 256, 512 and 1024.The initial normalized frequency in this experiment was ν = 0.547.From the figure, one can conclude that the boundary errors are similar in all the cases and they do not last more than about 10 samples independently of the number of pixels, or about 10 sampling intervals.

Resolving power of integrators
Resolving power of integrators characterizes their capability to resolve between close sharp impulses in the integrated data.It is fully defined by the integrator frequency responses.However, it is much more straightforward to compare the resolving power for different integrators directly in signal domain.Figure 5 illustrates results of numerical evaluation of the capabilities of three types of integrators, trapezoidal, cubic spline and DFT-based ones, to reproduce two sharp impulses placed on the distance of one sample one from another for the case when the second impulse is half height of the first impulse.The signals are 8 times subsampled to imitate the corresponding continuous signals at the integrator output.
The figure clearly shows that the tested integrators differ in their resolving power.DFTbased integrator produces the sharpest peaks with the lowest valley between the peaks while cubic spline integrators and trapezoidal integrators exhibit poorer behavior.In particular, the latter seems to be incapable of reliable resolving the half height impulse against oscillations seen in the sub-sampled signal.

Conclusion
We presented results of analytical and experimental comparison of frequency responses and resolving power of five methods for numerical integration of sampled data: trapezoidal method, Simpson method, Simpson-3/8 method, cubic spline method and DFT-based method.
We have shown that DFT based method provides the best numerical approximation to the ideal continuous integrator and outperforms other integrators in terms of the resolving power, although it is vulnerable to boundary effects due to the fact that it implements signal cyclic convolution rather than regular shift invariant convolution.We suggested a modification of this method, a DCT-based one that assumes signal extension by "mirror reflection" on its boundaries and that can be efficiently implemented using Fast DCT and DcST transforms.We have shown that the DCT-based method, being even slightly better than the DFT-based method in terms of approximation of the ideal continuous integrator, is also substantially more robust to boundary effects.We also have shown that boundary effect errors for both methods do not propagate more than to about 10 first and last signal samples, or about 10 sampling intervals.

( ) ( ) ( )
is numerically evaluated in computers using samples { } k a of signal ( ) Then Eq. (A1-4) can be rewritten as (A1-6) Therefore, PSF of the equivalent continuous filter is As one can see from this equation, continuous filter equivalent to the given digital filter of Eq. (A1-2) is not shift invariant.The reason lies, as it will be clear from what follows, in the finiteness of the number of signal samples.
It is more convenient to characterize equivalent continuous filter by its frequency response found as Fourier Transform of its impulse response over both variables x , ξ : .
(A1-8) This expression contains 4 terms:  , the filter is shift invariant as As we will see later, the shift parameter is needed in order to coordinate indices of signal samples with the continuous coordinate system of the discrete frequency response.We refer to this modification of DFT as ( ) 0 ; u SDFT [11].
Replace { } n h in Eq. (A1-10) by inverse ( ) (A2-3) is Discrete Cosine Transform (DCT).Therefore, the DFT spectrum of signal extended by "mirror reflection" can be computed via DCT using Fast DCT algorithm.From properties of DCT it follows that DCT spectra feature the following symmetry property: where asterisk symbolizes complex conjugate.By using Eqs.(A2-4) and (A2-6), obtain that: .Both transforms can be computed using fast algorithms [11].
its Fourier transform spectrum, respectively, and

Fig 1 .
Fig 1.Comparison of frequency responses of trapezoidal, Simpson, 3/8 Simpson, Cubic Splines and Fourier methods of integrations

Fig. 2 .
Fig. 2. Integration error of periodic sinusoidal signals as a function of the normalized frequency: (a) for all methods; (b) only for DFT-based, trapezoidal and cubic spline methodsFigures 2(a), and (b), show that the trapezoidal method and the method based on the cubic spline interpolation give similar errors in all frequencies, lower for the cubic spline based method.The error produced by the Fourier integration method for considered periodic signals is defined only by computation round-off errors.Obviously, all these results are in agreement with the behaviour of the frequency responses of the methods shown in figure1.
of the analytical signal.

Fig. 4 .
Fig. 4. Average error evaluated in the 10 first samples of the domain for the same initial normalized frequency (p = 0.547) but different N: (a) N = 256, (b) N = 512, (c) N = 1024.The black curve corresponds to Fourier integration method, the blue one to the DCT-based method and the red one to the cubic spline based method.

a
operation is commonly called digital filtering and the set of h N weight coefficients { } n h is called point spread function of the digital filter.Let us consider the correspondence between digital filter of Eq. (A1-2) and convolution integral of Eq. (A1-1) by finding the point spread function ( ) x h of a continuous filter that corresponds to a given set of digital filter coefficients { } n h .Results { } k b of computations according to Eq. (A1-2) are regarded as samples of a continuous signal ( ) x b that can be generated from the set of samples { } of interpolation functions that, by convention, reconstruct signal ( ) shift (in units of the sampling interval x ∆ ) of sample 0 b on the reconstruction device with respect to the origin of signal ( ) x b coordinate system, h N is the number of samples { } k b used for the reconstruction.By representing in Eq. (A1-3) samples { } k b as defined by Eq. (A1-with respect to the origin of signal ( ) the number b N of output signal samples.It reflects the fact that digital filter defined by Eqs.(A2) and obtained as a discrete representation of the convolution integral is nevertheless not shift invariant owing to boundary effects associated with the finiteness of the number b N of samples of the filter output signal { } k b involved in the reconstruction of continuous output signal ( ) x b .When the number of signal samples b N increases, contribution of the boundary effects into the reconstructed signal decreases and the quality of the digital filter approximation to the convolution integral improves.In the limit, when ∞ → b N are frequency responses of the signal reconstruction and discretization devices, respectively.If signal reconstruction and discretization devices are, as it is required by the sampling theorem, ideal low pass filters, terms ( ) ( ) is not the case for all real devices, one should anticipate aliasing effects similar to those in signal reconstruction.Show that discrete frequency response of the digital filter ( )p DFRis directly related to coefficients of the Discrete Fourier Transform of the filter discrete PSF { } n h .Compute DFT with a shift parameter u introduced in signal domain as:

Appendix 2 .
19) It follows from Eq. (A1-18) that values of the discrete frequency response ( ) the discrete sinc-interpolation function of Eq. (A1-19).Computation of convolution of signals with "mirror reflection" extension using Discrete Cosine and Discrete Cosine/Sine Transforms.Consider a signal extended to its double length by "mirror reflection":

η
computing convolution, the signal spectrum defined by Eq. (A2-2) should be multiplied by the filter frequency response for signals of N 2 samples and then the inverse DFT should be computed for the first N samples: are filter frequency response coefficients.The coefficients { } r η are the DFT of samples of the filter PSF that are real numbers.They feature the symmetry property: of this expression constitute inverse DCT of the product ( ) is Discrete Cosine/Sine Transform (DcST) of the product( ) ).For different integer frequencies int p , sinusoidal signals were generated with parameter