Abstract

We describe a fast computational algorithm able to evaluate the Rayleigh–Sommerfeld diffraction formula, based on a special formulation of the convolution theorem and the fast Fourier transform. What is new in our approach compared to other algorithms is the use of a more general type of convolution with a scale parameter, which allows for independent sampling intervals in the input and output computation windows. Comparison between the calculations made using our algorithm and direct numeric integration show a very good agreement, while the computation speed is increased by orders of magnitude.

© 2009 Optical Society of America

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  1. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 3.
  2. J. W. Cooley and J. W. Tookey, “An algorithm for the machine computation of the complex Fourier series,” Math. Comput. 19, 297-301 (1965).
    [CrossRef]
  3. L. P. Yaroslavsky, “Discrete transforms, fast algorithms, and point spread functions of numerical reconstruction of digitally recorded holograms,” in Advances in Signal Transforms: Theory and Applications, J. Astola and L. Yaroslavsky, eds. (Hindawi2007), pp. 93-141.
  4. J. A. C. Veerman, J. J. Rusch, and H. P. Urbach, “Calculation of the Rayleigh-Sommerfeld diffraction integral by exact integration of the fast oscillating factor,” J. Opt. Soc. Am. A 22, 636-646 (2005).
    [CrossRef]
  5. F. Shen and A. Wang, “Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula,” Appl. Opt. 45, 1102-1110(2006).
    [CrossRef] [PubMed]
  6. J. E. Harvey, “Fourier treatment of near-field scalar diffraction theory,” Am. J. Phys. 47, 974-980 (1979).
    [CrossRef]
  7. N. Delen and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast Fourier transform approach,” J. Opt. Soc. Am. A 15, 857-867 (1998).
    [CrossRef]
  8. L. B. Klebanov and J. F. Crouzet, “Quasi-convolutions and applications to coded images,” J. Math. Sci. 99, 1120-1126(2000).
    [CrossRef]
  9. D. Mustard, “Fractional convolution,” J. Austral. Math. Soci. Ser. B 40, 257-265 (1998).
    [CrossRef]
  10. P. C. Logofătu and D. Apostol, “The Fourier transform in optics: from continuous to discrete or from analogous experiment to digital calculus,” J. Optoelectron. Adv. Mat. 9, 2838-2846(2007).
  11. M. Mandal and A. Asif, Continuous and Discrete Time Signals and Systems (Cambridge U. Press, 2007).
  12. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 2.
  13. “An introduction to the sampling theorem,” National Semiconductor Application Note 236 (1980), http://www.national.com/an/AN/AN-236.pdf.
  14. D. H. Bailey and P. N. Swarztrauber, “The fractional Fourier transform and applications,” SIAM Rev. 33, 389-404(1991).
    [CrossRef]
  15. A. Bultheel and H. Martinez, “Computation of the fractional Fourier transform,” Appl. Comput. Harmon. Anal. 16, 182-202 (2004).
    [CrossRef]

2007 (1)

P. C. Logofătu and D. Apostol, “The Fourier transform in optics: from continuous to discrete or from analogous experiment to digital calculus,” J. Optoelectron. Adv. Mat. 9, 2838-2846(2007).

2006 (1)

2005 (1)

2004 (1)

A. Bultheel and H. Martinez, “Computation of the fractional Fourier transform,” Appl. Comput. Harmon. Anal. 16, 182-202 (2004).
[CrossRef]

2000 (1)

L. B. Klebanov and J. F. Crouzet, “Quasi-convolutions and applications to coded images,” J. Math. Sci. 99, 1120-1126(2000).
[CrossRef]

1998 (2)

1991 (1)

D. H. Bailey and P. N. Swarztrauber, “The fractional Fourier transform and applications,” SIAM Rev. 33, 389-404(1991).
[CrossRef]

1979 (1)

J. E. Harvey, “Fourier treatment of near-field scalar diffraction theory,” Am. J. Phys. 47, 974-980 (1979).
[CrossRef]

1965 (1)

J. W. Cooley and J. W. Tookey, “An algorithm for the machine computation of the complex Fourier series,” Math. Comput. 19, 297-301 (1965).
[CrossRef]

Apostol, D.

P. C. Logofătu and D. Apostol, “The Fourier transform in optics: from continuous to discrete or from analogous experiment to digital calculus,” J. Optoelectron. Adv. Mat. 9, 2838-2846(2007).

Asif, A.

M. Mandal and A. Asif, Continuous and Discrete Time Signals and Systems (Cambridge U. Press, 2007).

Bailey, D. H.

D. H. Bailey and P. N. Swarztrauber, “The fractional Fourier transform and applications,” SIAM Rev. 33, 389-404(1991).
[CrossRef]

Bultheel, A.

A. Bultheel and H. Martinez, “Computation of the fractional Fourier transform,” Appl. Comput. Harmon. Anal. 16, 182-202 (2004).
[CrossRef]

Cooley, J. W.

J. W. Cooley and J. W. Tookey, “An algorithm for the machine computation of the complex Fourier series,” Math. Comput. 19, 297-301 (1965).
[CrossRef]

Crouzet, J. F.

L. B. Klebanov and J. F. Crouzet, “Quasi-convolutions and applications to coded images,” J. Math. Sci. 99, 1120-1126(2000).
[CrossRef]

Delen, N.

Goodman, J. W.

J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 3.

J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 2.

Harvey, J. E.

J. E. Harvey, “Fourier treatment of near-field scalar diffraction theory,” Am. J. Phys. 47, 974-980 (1979).
[CrossRef]

Hooker, B.

Klebanov, L. B.

L. B. Klebanov and J. F. Crouzet, “Quasi-convolutions and applications to coded images,” J. Math. Sci. 99, 1120-1126(2000).
[CrossRef]

Logofatu, P. C.

P. C. Logofătu and D. Apostol, “The Fourier transform in optics: from continuous to discrete or from analogous experiment to digital calculus,” J. Optoelectron. Adv. Mat. 9, 2838-2846(2007).

Mandal, M.

M. Mandal and A. Asif, Continuous and Discrete Time Signals and Systems (Cambridge U. Press, 2007).

Martinez, H.

A. Bultheel and H. Martinez, “Computation of the fractional Fourier transform,” Appl. Comput. Harmon. Anal. 16, 182-202 (2004).
[CrossRef]

Mustard, D.

D. Mustard, “Fractional convolution,” J. Austral. Math. Soci. Ser. B 40, 257-265 (1998).
[CrossRef]

Rusch, J. J.

Shen, F.

Swarztrauber, P. N.

D. H. Bailey and P. N. Swarztrauber, “The fractional Fourier transform and applications,” SIAM Rev. 33, 389-404(1991).
[CrossRef]

Tookey, J. W.

J. W. Cooley and J. W. Tookey, “An algorithm for the machine computation of the complex Fourier series,” Math. Comput. 19, 297-301 (1965).
[CrossRef]

Urbach, H. P.

Veerman, J. A. C.

Wang, A.

Yaroslavsky, L. P.

L. P. Yaroslavsky, “Discrete transforms, fast algorithms, and point spread functions of numerical reconstruction of digitally recorded holograms,” in Advances in Signal Transforms: Theory and Applications, J. Astola and L. Yaroslavsky, eds. (Hindawi2007), pp. 93-141.

Am. J. Phys. (1)

J. E. Harvey, “Fourier treatment of near-field scalar diffraction theory,” Am. J. Phys. 47, 974-980 (1979).
[CrossRef]

Appl. Comput. Harmon. Anal. (1)

A. Bultheel and H. Martinez, “Computation of the fractional Fourier transform,” Appl. Comput. Harmon. Anal. 16, 182-202 (2004).
[CrossRef]

Appl. Opt. (1)

J. Austral. Math. Soci. Ser. B (1)

D. Mustard, “Fractional convolution,” J. Austral. Math. Soci. Ser. B 40, 257-265 (1998).
[CrossRef]

J. Math. Sci. (1)

L. B. Klebanov and J. F. Crouzet, “Quasi-convolutions and applications to coded images,” J. Math. Sci. 99, 1120-1126(2000).
[CrossRef]

J. Opt. Soc. Am. A (2)

J. Optoelectron. Adv. Mat. (1)

P. C. Logofătu and D. Apostol, “The Fourier transform in optics: from continuous to discrete or from analogous experiment to digital calculus,” J. Optoelectron. Adv. Mat. 9, 2838-2846(2007).

Math. Comput. (1)

J. W. Cooley and J. W. Tookey, “An algorithm for the machine computation of the complex Fourier series,” Math. Comput. 19, 297-301 (1965).
[CrossRef]

SIAM Rev. (1)

D. H. Bailey and P. N. Swarztrauber, “The fractional Fourier transform and applications,” SIAM Rev. 33, 389-404(1991).
[CrossRef]

Other (5)

J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 3.

M. Mandal and A. Asif, Continuous and Discrete Time Signals and Systems (Cambridge U. Press, 2007).

J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 2.

“An introduction to the sampling theorem,” National Semiconductor Application Note 236 (1980), http://www.national.com/an/AN/AN-236.pdf.

L. P. Yaroslavsky, “Discrete transforms, fast algorithms, and point spread functions of numerical reconstruction of digitally recorded holograms,” in Advances in Signal Transforms: Theory and Applications, J. Astola and L. Yaroslavsky, eds. (Hindawi2007), pp. 93-141.

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)

Fig. 1
Fig. 1

Geometric parameters of the Rayleigh–Sommerfeld diffraction formula for the 1D case.

Fig. 2
Fig. 2

Input test function for diffraction computation.

Fig. 3
Fig. 3

Real part of the well-sampled convolution kernel corresponding to a propagation distance d = 140 μm , an output window Δ x = 50 μm (which together require a minimum number of samples N min = 39 ), and an overabundant number of samples N = 1024 . The imaginary part has a similar look.

Fig. 4
Fig. 4

Transformed test function using the well-sampled convolution kernel shown in Fig. 3.

Fig. 5
Fig. 5

Computation errors of the fast algorithm for the case of the well-sampled kernel (the square modulus of the difference with respect to the numeric integration).

Fig. 6
Fig. 6

Undersampled convolution kernel. It corresponds to the same propagation distance d = 140 μm , a much larger output window Δ x = 50 μm (which constitute a set of parameters requiring a minimum number of samples N min = 5006 ), and an insufficient number of samples N = 1024 .

Fig. 7
Fig. 7

Transformed test function using the poorly sampled convolution kernel shown in Fig. 6.

Fig. 8
Fig. 8

2D input test function for diffraction computation.

Fig. 9
Fig. 9

Real part of a well-sampled 2D convolution kernel. The imaginary part has a similar look.

Fig. 10
Fig. 10

Intensity profile of the diffraction figure u . The central high-intensity peak was cut so that the smaller details of the rest of the profile can be observed.

Equations (49)

Equations on this page are rendered with MathJax. Learn more.

u ( r ) = i λ Σ u ( r ) exp ( i k | Δ r | ) | Δ r | cos ( n , Δ r ) d σ ,
Δ r = r r , | Δ r | = [ ( x x ) 2 + ( y y ) 2 + ( z z ) 2 ] 1 / 2 .
u ( x ) = d i λ u ( x ) · exp [ i k r ( x , x ) ] r ( x , x ) 3 / 2 d x ,
r ( x , x ) | Δ r | = [ ( x x ) 2 + d 2 ] 1 / 2 .
u ( x ) = ( u ( x ) * h ( x ) ) ( x ) = u ( x ) h ( x x ) d x ,
h ( x ) = d i λ exp ( i k x 2 + d 2 ) ( x 2 + d 2 ) 3 / 4 .
u n = u ( x n ) = d i λ δ x n = N / 2 N / 2 1 u n exp ( i k r n n ) r n n 3 / 2 ,
r n n = [ ( x n x n ) 2 + d 2 ] 1 / 2 = [ ( n δ x n δ x ) 2 + d 2 ] 1 / 2 .
u n = n = N / 2 N / 2 1 u n h n n ,
h n = d i λ δ x exp ( i k n 2 δ x 2 + d 2 ) ( n 2 δ x 2 + d 2 ) 3 / 4 ,
u n = n = N / 2 N / 2 1 u n h α n α n ,
h α n α n = d i λ δ x exp ( i k ( α n α n ) 2 δ 2 + d 2 ) [ ( α n α n ) 2 δ 2 + d 2 ] 3 / 4 ,
M / 2 min ( α n α n ) ,
M / 2 1 max ( α n α n ) .
δ = max ( δ x , δ x ) ,
α = δ x / δ ,
α = δ x / δ .
α min = min ( α , α ) = min ( δ x , δ x ) / δ .
M = N ( 1 + α min ) .
u ( 0 0 0 ( M N ) / 2 u N / 2 u N / 2 + 1 u N / 2 1 N 0     0 0 ( M N ) / 2 ) ,
h m = d i λ δ x exp ( i k m 2 δ 2 + d 2 ) ( m 2 δ 2 + d 2 ) 3 / 4 ,
U p ( α ) = m = M / 2 M / 2 1 u m w M α m p ,
H p = n = M / 2 M / 2 1 h n w M n p ,
w M = exp ( i 2 π / M ) .
u m = ? 1 M p = M / 2 M / 2 1 U p ( α ) H p w M α m p =
= 1 M p = M / 2 M / 2 1 ( m = M / 2 M / 2 1 u m w M α m p ) ( n = M / 2 M / 2 1 h n w M n p ) w M α m p =
= m = M / 2 M / 2 1 u m [ n = M / 2 M / 2 1 h n ( 1 M p = M / 2 M / 2 1 w M p [ n ( α m α m ) ] ) ] =
= m = M / 2 M / 2 1 u m h α m α m .
g ( x ) = n = g n sinc ( n x / δ x ) ,
1 M p = M / 2 M / 2 1 w M p ( n x / δ x ) = exp [ i π / M ( n x / δ x ) ] sin [ π ( n x / δ x ) ] M sin [ π / M ( n x / δ x ) ] ,
M exp [ i π / M ( n x / δ x ) ] sin [ π ( n x / δ x ) ] M sin [ π / M ( n x / δ x ) ] sinc ( n x / δ x ) ,
Φ ( m ) = k ( m 2 δ 2 + d 2 ) 1 / 2 ,
ν ( m ) = 1 2 π d Φ d m = m δ 2 λ ( m 2 δ 2 + d 2 ) 1 / 2 .
ν max = | ν ( ± M / 2 ) | = max ( Δ x , Δ x ) ( Δ x + Δ x ) / 2 N λ [ ( Δ x + Δ x ) 2 / 4 + d 2 ] 1 / 2 < 1 2 .
N > max ( Δ x , Δ x ) λ 2 ( Δ x + Δ x ) / 2 [ ( Δ x + Δ x ) 2 / 4 + d 2 ] 1 / 2 =
= max ( Δ x , Δ x ) λ 2 sin θ max N min
u n = δ x i λ d exp [ i k ( d + n 2 δ x 2 2 d ) ] n = N / 2 N / 2 1 [ u n exp ( i k n 2 δ x 2 2 d ) ] exp [ i 2 π n n N α F r ] ,
α F r = δ x δ x N λ d .
δ x = λ d N δ x .
N Δ x Δ x λ d =
= Δ x · Δ x λ · 2 Δ x + Δ x · Δ x + Δ x 2 d = 1 λ ( 1 Δ x + 1 Δ x ) 1 2 tan θ max N min F r .
F n = m = M / 2 M / 2 1 f m w M α m n ,
F n = w M n 2 α / 2 m = M / 2 M / 2 1 f m w M m 2 α / 2 w M ( m n ) 2 / 2 .
F n = w M n 2 α / 2 m = M / 2 M / 2 1 g m h n m ,
g m = w M m 2 α / 2 f m ,
h m = w M m 2 α / 2 .
g ( 00 0 M / 2 g M / 2 g M / 2 + 1 g M / 2 1 M 00 0 M / 2 ) ,
h ( h M h M + 1 h M 1 2 M ) .
F n = w M α n 2 / 2 F n ext ,

Metrics