## Frequency responses and resolving power of numerical integration of sampled data

Optics Express, Vol. 13, Issue 8, pp. 2892-2905 (2005)

http://dx.doi.org/10.1364/OPEX.13.002892

Acrobat PDF (381 KB)

### Abstract

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 DFT- based 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

## 1. Introduction

1. W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. **70**, 998–1006 (1980) [CrossRef]

4. R. J. Noll, “Phase Estimates from Slope-Type Wave-Front Slope Measurements,” J. Opt. Soc. Am. **68**, 139–140 (1978) [CrossRef]

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) [CrossRef] [PubMed]

6. F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” Appl. Opt. **30**, 1325–1327 (1991) [CrossRef] [PubMed]

7. W. Zou and Z. Zhang, “Generalized wave-front reconstruction algorithm applied in a Sach-Hartmann test,” Appl. Opt. **39**, 250–268 (2000) [CrossRef]

8. S. Krey, W.D. van Amstel, J. Campos, A. Moreno, and E.J. Lous, “A fast optical scanning deflectometer for measuring the topography of large silicon wafers,” Proc. SPIE, 5523, (*SPIE’s 49*^{th}* Annual Meeting, Denver*), (2004) [CrossRef]

9. J. Campos, L. P. Yaroslavsky, A. Moreno, and 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) [CrossRef]

10. A. Moreno, J. Campos, and L.P. Yaroslavsky, “Frequency response of five integration methods to obtain the profile from its slope,” Optical Engineering, (2005) (to be published) [CrossRef]

## 2. Continuous and numerical integrators and their frequency responses

*a*(

*x*) and

*α*(

*f*) are signal and its Fourier transform spectrum, respectively, and

*a*

_{k}},

*k*=0,1,…,

*N*-1 is a set of

*samples of the input signal to be integrated, DFT{.} and IDFT{.} are operators of direct and inverse Discrete Fourier Transforms, sign*

**N****and**

*N***[11], where**

*N**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.

*” 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.*

**f**## 3. DCT based integrator

*ā*

_{k}}, 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

**should be replaced by 2**

*N**N*. The value of

**-th DFT spectral coefficient of symmetrical signals such as {**

*N**ã*

_{k}} is equal to zero (see Appendix 2, Eq. (A2-4)).

**” 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.**

*f**N*- 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

**samples. We will refer to this modified DFT based integration method as “extended”, or DCT-based method.**

*N*## 4. Experimental comparison

**-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.**

*N*### 4.1 Integration of periodic sinusoidal signals

**is the number of signal samples,**

*N***is the frequency parameter of the sinusoidal signal, the normalized frequency is given by**

*p**corresponds to an initial random phase and*

**2p/N; randphase***represents the domain of the signal (*

**x****0≤**). 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

*x*≤*N**f*

_{i}(

*k*) is the obtained function with the numerical integration and

*f*

_{0}(

*k*) 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

**2**/

*p***, 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.**

*N*### 4.2 Aperiodic signals and boundary effects

*p*in Eq. (22). For different integer frequencies

*p*

_{int}, sinusoidal signals were generated with parameter

*p*=

*p*

_{int}+

*s*/

*n*,

*s*=0,…,

*n*with values equally spaced in the interval

*p*

_{int},

*p*

_{int}+1. For each

*p*

_{int}, the error was evaluated in every signal sample

*k*as

*k*)} are samples of the numerically integrated function for each frequency and {

*k*)} are corresponding samples of the analytical signal.

**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.**

*n**ν*=0.547 (medium frequency) is shown in Fig. 3(b). In this case the boundary effect error for DFT-based 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.

*ν*=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.

*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.

### 4.3 Resolving power of integrators

## 5. Conclusion

## Appendix 1. Discrete representation of convolution integrals and point spread function and frequency response of digital filtering

*a*(

*x*) with shift invariant kernel

*h*(

*x*) (point spread function, PSF):

*a*

_{k}} of signal

*a*(

*x*) as

**weight coefficients {**

*N*_{h}*h*

_{n}} 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

*h*(

*x*) of a continuous filter that corresponds to a given set of digital filter coefficients{

*h*

_{n}}. Results {

*b*

_{k}} of computations according to Eq. (A1-2) are regarded as samples of a continuous signal

*b*(

*x*) that can be generated from the set of samples {

*b*

_{k}} as

*φ*

^{(r}(·)} is a set of interpolation functions that, by convention, reconstruct signal

*b*(

*x*) from samples {

*b*

_{k}} with certain admitted accuracy,

*k̃*=

*k*+

*u*

^{(r)},

*u*

^{(r)}is a position shift (in units of the sampling interval Δ

*x*) of sample

*b*

_{0}on the reconstruction device with respect to the origin of signal

*b*(

*x*) coordinate system,

_{h}**is the number of samples {**

*N**b*

_{k}} used for the reconstruction. By representing in Eq. (A1-3) samples {

*b*

_{k}} as defined by Eq. (A1-2) one obtain

*a*

_{m}} of signal

*a*(

*x*) are obtained by means of a sampling device with point spread function

*φ*

^{(d)}[

*x*-

*m͌*Δ

*x*], where

*m͌*=

*m*+

*u*

^{(d)}and

*u*

^{(d)}is the position shift of sample

*a*

_{0}with respect to the origin of signal

*a*(

*x*) coordinate system:

*x*, ξ:

*DFR*(

*p*)

*h*

_{n}} and is referred to as

**of the digital filter. It is a periodic function with period1/Δ**

*discrete frequency response**x*.

*SV*(

*f*,

*p*)

*N*

_{b}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

*N*

_{b}of samples of the filter output signal {

*b*

_{k}} involved in the reconstruction of continuous output signal

*b*(

*x*). When the number of signal samples

*N*

_{b}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

*N*

_{b}→∞, the filter is shift invariant as

^{(r)}(

*f*) and Φ

^{(d)}(-

*p*) remove all but one periods of the discrete frequency response

*DFR*(

*p*). Because this is not the case for all real devices, one should anticipate aliasing effects similar to those in signal reconstruction.

*SDFT*(

*u*;0) [11].

*u*=-

*u*

^{(d)}=-(

*N*

_{h}-1)2/h, phase shift factors in Eq. (A1-17) are compensated and:

*DFR*(

*f*) of the digital filter at points

*f*=

*r*

*N*

_{h}Δ

*x*,

*r*=0,…,

*N*

_{h}-1 within one period 1 Δ

*x*of its periodicity are proportional to

*SDFT*[-(

*N*

_{h}-1)/2;0] coefficients {

*η*

_{r}) of the interpolation kernel {

*h*

_{n}}. Between these sampling points,

*DFR*(

*f*) is interpolated with the discrete sinc-interpolation function of Eq. (A1-19).

## Appendix 2. Computation of convolution of signals with “mirror reflection” extension using Discrete Cosine and Discrete Cosine/Sine Transforms.

*N*samples and then the inverse DFT should be computed for the first

**samples:**

*N**η*

_{r}} 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:

*η*

_{r}}, we obtain finally:

## Acknowledgments

## References and links

1. | W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. |

2. | D.L. Fried, “Least-Square Fitting a Wave-Front Distortion Estimate to an Array of Phase-Difference Measurements,” J. Opt. Soc. Am. , |

3. | R. H. Hudgin, “Wave-Front Reconstruction for Compensated Imaging,” J. Opt. Soc Am. |

4. | R. J. Noll, “Phase Estimates from Slope-Type Wave-Front Slope Measurements,” J. Opt. Soc. Am. |

5. | A. Dubra, C. Paterson, and C. Dainty, “Wavefront reconstruction from shear phase maps using the discrete Fourier transform,” Appl. Opt. |

6. | F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” Appl. Opt. |

7. | W. Zou and Z. Zhang, “Generalized wave-front reconstruction algorithm applied in a Sach-Hartmann test,” Appl. Opt. |

8. | S. Krey, W.D. van Amstel, J. Campos, A. Moreno, and E.J. Lous, “A fast optical scanning deflectometer for measuring the topography of large silicon wafers,” Proc. SPIE, 5523, ( |

9. | J. Campos, L. P. Yaroslavsky, A. Moreno, and M.J. Yzuel, “Integration in the Fourier Domain for restoration of a function from its slope: Comparison of four methods,” Opt. Lett. |

10. | A. Moreno, J. Campos, and L.P. Yaroslavsky, “Frequency response of five integration methods to obtain the profile from its slope,” Optical Engineering, (2005) (to be published) [CrossRef] |

11. | L. Yaroslavsky, |

12. | J. H. Mathew and K. D. Fink, |

13. | C. Elster and I. Weingärtner, “High-accuracy reconstruction of a function f(x) when only df(x)/dx is known at discrete measurements points,” Proc. SPIE4782, (2002) |

**OCIS Codes**

(070.2590) Fourier optics and signal processing : ABCD transforms

(120.6650) Instrumentation, measurement, and metrology : Surface measurements, figure

(200.3050) Optics in computing : Information processing

(240.6700) Optics at surfaces : Surfaces

**ToC Category:**

Research Papers

**History**

Original Manuscript: January 19, 2005

Revised Manuscript: March 30, 2005

Published: April 18, 2005

**Citation**

L. Yaroslavsky, A. Moreno, and J. Campos, "Frequency responses and resolving power of numerical integration of sampled data," Opt. Express **13**, 2892-2905 (2005)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-8-2892

Sort: Journal | Reset

### References

- W. H. Southwell, �??Wave-front estimation from wave-front slope measurements,�?? J. Opt. Soc. Am. 70, 998-1006 (1980) [CrossRef]
- 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) [CrossRef]
- R. H. Hudgin, �??Wave-Front Reconstruction for Compensated Imaging,�?? J. Opt. Soc. Am. 67, 375-378 (1977) [CrossRef]
- R. J. Noll, �??Phase Estimates from Slope-Type Wave-Front Slope Measurements,�?? J. Opt. Soc. Am. 68, 139-140 (1978) [CrossRef]
- A. Dubra, C. Paterson, and C. Dainty, "Wavefront reconstruction from shear phase maps using the discrete Fourier transform," Appl. Opt. 43, 1108-1113 (2004) [CrossRef] [PubMed]
- F. Roddier and C. Roddier, �??Wavefront reconstruction using iterative Fourier transforms,�?? Appl. Opt. 30, 1325-1327 (1991) [CrossRef] [PubMed]
- W. Zou and Z. Zhang, �??Generalized wave-front reconstruction algorithm applied in a Sach-Hartmann test,�?? Appl. Opt. 39, 250-268 (2000) [CrossRef]
- 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 49th Annual Meeting, Denver), (2004) [CrossRef]
- 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) [CrossRef]
- 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) [CrossRef]
- L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004)
- J. H. Mathew and K. D. Fink, Numerical Methods using MATLAB (Prentice-Hall, Englewood Cliffs, N. J., 1999)
- 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)

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.