OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 13, Iss. 8 — Apr. 18, 2005
  • pp: 2892–2905
« Show journal navigation

Frequency responses and resolving power of numerical integration of sampled data

L. P. Yaroslavsky, A. Moreno, and J. Campos  »View Author Affiliations


Optics Express, Vol. 13, Issue 8, pp. 2892-2905 (2005)
http://dx.doi.org/10.1364/OPEX.13.002892


View Full Text Article

Acrobat PDF (381 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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

Numerical integration has numerous applications in several optical fields as the wave-front reconstruction from wave-front slope measurements [1

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

4

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

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

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]

]. In [6

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

] 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

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

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

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 49th Annual Meeting, Denver), (2004) [CrossRef]

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

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

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]

].

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.

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

a¯(x)=a(x)dx=i2πfα(f)exp(i2πfx)dx=Hint(f)α(f)exp(i2πfx)df.
(1)

where a(x) and α(f) are signal and its Fourier transform spectrum, respectively, and

Hint(f)=i2πf
(2)

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:

{a¯k}=IDET{ηr(int)DET{ak}},
(3)

where {ak }, k=0,1,…,N-1 is a set of N samples of the input signal to be integrated, DFT{.} and IDFT{.} are operators of direct and inverse Discrete Fourier Transforms, sign

• denotes element-wise product of vectors, and

ηr(int)={0,r=0Nhi2πr,r=1,2,N21h2π,r=N2ηNr*,r=N2+1,,N1,
(4a)

for even N and

ηr(int)={0,r=0Nhi2πr,r=1,2,,(N1)2,ηNr*,r=N2+1,,N1
(4b)

for odd N [11

11. L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004)

], where h is the signal sampling interval. Note that digital signal integrating according to Eq. (3) automatically implies signal discrete sinc-interpolation [11

11. L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004)

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

11. L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004)

], [12

12. J. H. Mathew and K. D. Fink, Numerical Methods using MATLAB (Prentice-Hall, Englewood Cliffs, N. J., 1999)

] and [13

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)

]). 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:

a¯1(T)=0,a¯k(T)=a¯k1(T)+h2(ak1+ak)
(5)

In the Simpson rule, the second point needs to be evaluated by the trapezoidal rule, then

a¯1(S)=0,a¯k(S)=a¯k2(S)+h3(ak2+4ak1+ak)
(6)

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

a¯0(38S)=0,a¯k(38S)=ak3(38S)+3h8(ak3+3ak2+3ak1+ak)
(7)

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

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)

], 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:

a¯k(Tr)a¯k1(Tr)=h2(ak1+ak)
(8)
a¯k(S)a¯k2(S)=h3(ak2+4ak1+ak)
(9)
a¯k(38S)a¯k3(38S)=3h8(ak3+3kk2+3ak1+ak)
(10)

Discrete point spread function for the method based on cubic spline interpolation of the input data information can be found as follows:

a¯k+1(CS)a¯k(CS)=h[12(ak+ak+1)h224(mk+mk+1)]
(11)

where the coefficients are determined by the system of linear equations:

h(mk1+4mk+mk+1)=6h(ak+12ak+ak1)
(12)

Taking N point Discrete Fourier Transform of these expressions (Eqs. (8)(12)) we obtain

α¯r(T)(1exp(i2πrN))=h2αr[1+exp(i2πrN)];
(13)
α¯r(S)(1exp(i4πrN))=h3αr[1+4exp(i2πrN)+exp(i4πrN)];
(14)
α¯r(38S)[1exp(i6πrN)]=3h8αr[1+3exp(i2πrN)+3exp(i4πrN)+exp(i6πrN)],
(15)
α¯r(CS)=h4cos(πrN)isin(πrN)[1+3cos(2πrN)+2]αr
(16)

Therefore the frequency responses of these integrators are, respectively:

ηr(T)=α¯r(Tr)αr={0,r=0,hcos(πrN)2isin(πrN),r=1,,N1,
(17)
ηr(S)=α¯r(S)αr={0,r=0,hcos(2πrN)+23isin(2πrN),r=1,,N1,
(18)
ηr(3S)=α¯r(3S)αr={0,r=0,hcos(3πrN)+3cos(πrN)isin(3πrN),r=1,,N1.
(19)
ηr(CS)=α¯r(CS)αr={0,r=0h4icos(πrN)sin(πrN)[1+3cos(2πrN)+2],r=1,,N1
(20)

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.

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

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.

3. 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:

a˜k={ak,k=0,1,,N1a2N1k,k=N,,2N1
(21)

Such an extended signal {ā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 N should be replaced by 2N. The value of ηN(int) is inessential as N-th DFT spectral coefficient of symmetrical signals such as {ã k } is equal to zero (see Appendix 2, Eq. (A2-4)).

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.

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

4.1 Integration of periodic sinusoidal signals

In these experiments, sinusoidal signals and their derivatives are generated as, correspondingly:

f0(x)=cos(2πpNx+randomphase)
f0(x)=2πpNcos(2πpNx+randomphase)
(22)

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 (0≤xN). 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

error=1Nk=0N[fi(k)f0(k)]2
(23)

where fi (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 2p/N, 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.

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 methods

Figures 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 figure 1.

4.2 Aperiodic signals and boundary effects

In order to study the boundary effects for the best methods (DFT-based method (FI), DCT-based 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. (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

errorpint(k)=1ns=0n[fipint+sn(k)f0pint+sn(k)]2
(24)

where {fip (k)} are samples of the numerically integrated function for each frequency and {f0p(k)} are corresponding samples of the analytical signal.

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

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.

Fig. 3. Experimentally obtained integration error versus sample k for DFT-based (FI) method (black) CSI method (red) and DCT-based (Extended) method (blue). Normalized initial frequency: (a) ν=0.273. (b) ν=0.547 and (c) ν=0.820

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.

4.3 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 sub-sampled to imitate the corresponding continuous signals at the integrator output.

The figure clearly shows that the tested integrators differ in their resolving power. DFT-based 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.

Fig. 5. Theoretical Profile(a) and Integrated Profiles with Trapezoidal Rule (b), Cubic Spline method (c), and DFT- Method (d). All are shown subsampled.

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

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

The convolution integral of a signal a(x) with shift invariant kernel h(x) (point spread function, PSF):

b(x)=a(ξ)h(xξ)dξ.
(A1-1)

is numerically evaluated in computers using samples {ak } of signal a(x) as

bk=n=0Nh1hnakn
(A1-2)

This operation is commonly called digital filtering and the set of Nh weight coefficients {hn } 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{hn }. Results {bk } 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 {bk } as

b(x)=k=0Nb1bkφ(r)(xk˜Δx),
(A1-3)

where {φ(r (·)} is a set of interpolation functions that, by convention, reconstruct signal b(x) from samples {bk } with certain admitted accuracy, =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 N is the number of samples {bk } used for the reconstruction. By representing in Eq. (A1-3) samples {bk } as defined by Eq. (A1-2) one obtain

b(x)=k=0Nb1(n=0Nh1hnakn)φ(r)(xk˜Δx),
(A1-4)

Let us suppose that samples {am } of signal a(x) are obtained by means of a sampling device with point spread function φ (d)[x-Δx], where =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:

am=a(ξ)φ(r)[ξmΔx]dξ
(A1-5)

Then Eq. (A1-4) can be rewritten as

b(x)=k=0Nb1[n=0Nh1hna(ξ)φ(d)[ξ(k˜n)Δx]dξ]φ(r)(xk˜Δx)=
a(ξ)[k=0Nbqn=0Nh1hnφ(r)(xk˜Δx)φ(d)[ξ(k˜n)Δx]]dξ.
(A1-6)

Therefore, PSF of the equivalent continuous filter is

heq(x,ξ)=k=0Nb1n=0Nh1hnφ(r)(xk˜Δx)φ(d)[ξ(k˜n)Δx]
(A1-7)

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, ξ:

Heq(f,p)=heq(x,ξ)exp[i2π(fxpξ)]dxdξ=
[n=0Nh1hnexp(i2πpnΔx)][k=0Nb1exp[i2π(fp)k˜Δx]]×
φ(r)(x)exp(i2πfx)dxφ(d)(ξ)exp(i2πpξ)dξ.
(A1-8)

This expression contains 4 terms:

Heq(f,p)=DFR(p)·SV(f,p)·Φ(r)(f)·Φ(r)(p).
(A1-9)

The term DFR(p)

DFR(p)=n=0Nh1hnexp(i2πpn˜Δx)=exp(i2πpu(d)Δx)n=0Nh1hnexp(i2πpnΔx)
(A1-10)

is determined solely by digital filter discrete point spread function {hn } and is referred to as discrete frequency response of the digital filter. It is a periodic function with period1/Δx.

Term SV(f, p)

SV(f,p)=k=0Nb1exp[i2π(fp)k˜Δx]=sin[π(fp)NbΔx]sin[π(fp)Δx]exp[iπ(fp)(Nb+u(r)1)Δx]
(A1-11)

depends only on the number Nb 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 Nb of samples of the filter output signal {bk } involved in the reconstruction of continuous output signal b(x). When the number of signal samples Nb 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 Nb →∞, the filter is shift invariant as

limNbNbsin[π(fp)NbΔx]Nbsin[π(fp)Δx]=δ(fp)
(A1-12)

Last two terms

Φ(r)(f)=φ(r)(x)exp(i2πfx)dx
(A1-13)

and

Φ(d)(p)=φ(d)(x)exp(i2πpx)dx
(A1-14)

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Φ(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.

ηr=1Nhn=0Nh1hnexp(i2πn+uNhr)
(A1-15)

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 SDFT(u;0) [11

11. L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004)

].

Replace {hn } in Eq. (A1-10) by inverse SDFT(u;0) of {ηr }

hn=r=0Nh1ηnexp(i2πn+uNhr)=1Nhr=0Nh1[ηrexp(i2πurNh)]exp[i2πnrNh]
(A1-16)

and obtain:

DFR(f)n=0Nh1r=0Nh1ηrexp{i2π[(fΔxrNh)n+(fΔxu(d)urNh)]}=
r=0Nh1ηnexp[i2π(fΔxu(d)urNh)]n=0Nh1exp[i2π(fΔxrNh)n]=
r=0Nh1ηrsin[π(fNhΔxr)]sin[π(fNhΔxr)Nh]exp[iπ(Nh12u(d))fΔx]exp[i2π(u+Nh12)rNh]
(A1-17)

Now one can see that, with settings u=-u (d)=-(Nh -1)2/h, phase shift factors in Eq. (A1-17) are compensated and:

DFR(f)r=0LN1ηrsincd[Nh;π(fNhΔxr)],
(A1-18)

where

sincd(N;x)=sinxNsin(xN).
(A1-19)

It follows from Eq. (A1-18) that values of the discrete frequency response DFR(f) of the digital filter at points f=r Nh Δx, r=0,…, Nh -1 within one period 1 Δx of its periodicity are proportional to SDFT[-(Nh -1)/2;0] coefficients {ηr ) of the interpolation kernel {hn }. 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.

Consider a signal extended to its double length by “mirror reflection”:

a˜k={ak,k=0,1,,N1a2N1k,k=N,,2N1.
(A2-1)

Its DFT spectrum is as follows:

α˜r=12Nk=02N1a˜kexp(i2πkr2N)=
{22Nk=0N1akcos[π(k+12)rN]}exp(iπr2N)=αr(DCT)exp(iπr2N),
(A2-2)

where

αr(DCT)=DCT{ak}=22Nk=0N1akcos(πk+12Nr)
(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:

αN(DCT)=0;αu(DCT)=α2Nu(DCT).
(A2-4)

For computing convolution, the signal spectrum defined by Eq. (A2-2) should be multiplied by the filter frequency response for signals of 2N samples and then the inverse DFT should be computed for the first N samples:

bk=12Nr=02N1αr(DCT)exp(iπr2N)ηrexp(i2πkr2N),
(A2-5)

where {η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=η2Nr*},
(A2-6)

where asterisk symbolizes complex conjugate. By using Eqs. (A2-4) and (A2-6), obtain that:

bk=12N{α0(DCT)η0+r=1N1αr(DCT)ηrexp(i2πk+122Nr)+
+r=1N1α2Nr(DCT)η2Nrexp[i2πk+122N(2Nr)]=
12N{α0(DCT)η0+r=1N1αr(DCT)ηr[exp(i2πk+122Nr)+ηr*exp(i2πk+122Nr)]
(A2-7)

As

ηrexp(i2πk+122Nr)+ηr*exp(i2πk+122Nr)=
2Re[ηrexp(i2πk+122Nr)]=ηrrecos(πk+12Nr)ηrimsin(πk+12Nr)
(A2-8)

where {ηrre } and {ηrim } are real and imaginary parts of {ηr }, we obtain finally:

bk=12N{α0(DCT)η0+r=1N1αr(DCT)ηrrecos(πk+12Nr)r=1N1αr(DCT)ηrimsin(πk+12Nr).
(A2-9)

First two terms of this expression constitute inverse DCT of the product {αr(DCT) ηrre } while the third term is Discrete Cosine/Sine Transform (DcST) of the product {αr(DCT) ηrim }. Both transforms can be computed using fast algorithms [11

11. L. Yaroslavsky, Digital Holography and Digital Image Processing (Kluwer Academic Publishers, 2004)

].

Acknowledgments

We would like to thank M. J. Yzuel for useful discussions. This work has been partially financed by Ministerio de Ciencia y Tecnología, under project BFM2003-006273-C02-01, and by the European Community project CTB556-01-4175. Alfonso Moreno acknowledges the grant receive by this European project. L. Yaroslavsky acknowledges the grant SAB2003-0053 by Spain Ministry of Education, Culture and Sport.

References and links

1.

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

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

3.

R. H. Hudgin, “Wave-Front Reconstruction for Compensated Imaging,” J. Opt. Soc Am. 67, 375–378 (1977) [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 49th 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]

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

  1. W. H. Southwell, �??Wave-front estimation from wave-front slope measurements,�?? J. Opt. Soc. Am. 70, 998-1006 (1980) [CrossRef]
  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) [CrossRef]
  3. R. H. Hudgin, �??Wave-Front Reconstruction for Compensated Imaging,�?? J. Opt. Soc. Am. 67, 375-378 (1977) [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, 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]
  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) [CrossRef]
  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) [CrossRef]
  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)

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.

CrossCheck Deposited