Computation of Hankel Transform using FFT (Fourier)

Question

To address circularly symmetric cases of 2-D Fourier Transformations, the so-called Hankel Transform can be applied (for a detailed derivation of the relation between the 2-D Fourier transform and the 1-D Hankel transform see Link).

For more complex or numerically tabulated radial functions, a numeric implementation of the Hankel transform (similiar to FFT) would come in handy. There has been an interesting proposal by Oppenheim from 1980 (full text is available on the Web) that explains the use of the highly optimized FFT algorithm to do the Hankel transform.

Is anybody aware of such an implementation of the Hankel transform in Mathematica?

Answer 1

Finally I found the most promising algorithm proposed in this really good reference Manuel Guizar-Sicairos and Julio C. Gutiérrez-Vega, "Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields," J. Opt. Soc. Am. A 21, 53-58 (2004). The authors call the algorithm pth-order quasi-discrete Hankel Transform (pQDHT) and it is perfectly suited for vector implementation within Mathematica.

The references uses the following definitions for the Hankel Transform:

$$f_2(\nu)=2 \pi \int_0^{\infty } f_1(r) J_p(2 \pi \nu r) r \, dr$$

and its inverse

$$f_1(r)=2 \pi \int_0^{\infty } f_2(\nu) J_p(2 \pi r \nu) \nu \, d\nu$$

My code for the Transform and its inverse looks as follows:

DiscreteHankelTransform[f1_, p_, Np_, rmax_] := 
  Module[{a, aNp1, rv, uv, res, umax, T , J, F1, F2},
  a = Table[N[BesselJZero[p, n]], {n, 1, Np}];
  aNp1 = BesselJZero[p, Np + 1];
  umax = aNp1/(2 Pi rmax);
  rv = a/(2 Pi umax);
  uv = a/(2 Pi rmax);
  J = Abs[BesselJ[p + 1, a]];
  T = BesselJ[p, TensorProduct[a, a]/(2 Pi  rmax umax)]/(
  TensorProduct[J, J] Pi rmax umax) ;
  F1 = (If[MatchQ[Head[f1], List], f1, f1[rv]] rmax)/J;
  F2 = T.F1;
  Return[Transpose[{uv, J/umax F2}]];
];

InverseDiscreteHankelTransform[f2_, p_, Np_, umax_] := 
 Module[{a, aNp1, rv, uv, res, rmax, T , J, F1, F2},
 a = Table[N[BesselJZero[p, n]], {n, 1, Np}];
 aNp1 = BesselJZero[p, Np + 1];
 rmax = aNp1/(2 Pi umax);
 rv = a/(2 Pi umax);
 uv = a/(2 Pi rmax);
 J = Abs[BesselJ[p + 1, a]];
 T = BesselJ[p, TensorProduct[a, a]/(2 Pi rmax umax)]/(
  TensorProduct[J, J] Pi rmax umax) ;
 F2 = (If[MatchQ[Head[f2], List], f2, f2[uv]] umax)/J;
 F1 = T.F2;
 Return[Transpose[{rv, J/rmax F1}]];
];

with f1 and f2 defined as pure functions of r and u respectively. Both functions can be specified as lists of x and y values too. p is the order of the transform, Np is the number of points used for the non zero part of the functions and their transform and rmax and umax are the maximum radial values for which f1 is nonzero and umax the maximum spatial frequency for which f2 is nonzero respectively. A test case is the example function $\frac{\text{sinc}(2 \pi \gamma r)}{2 \pi \gamma r}$ mentioned in the above reference. This test can be carried out as

f = Sin[2 Pi g #]/(2 Pi g #) &;
F0[u_, p_, g_] = 
 Piecewise[{
   {HankelTransform[f, u, p, {p >= 0, 0 < u < g}], 0 <= u < g},
   {Sin[p ArcSin[g/u]]/(2 Pi g Sqrt[u^2 - g^2]), u > g},
   {Infinity, u == g}
  }, Null]
SetAttributes[F0, Listable];

dF01 = DiscreteHankelTransform[f /. g -> 5, 1, 256, 3.0];
dF04 = DiscreteHankelTransform[f /. g -> 5, 4, 256, 3.0];
dF01max = Max[Abs[dF01]];
dF04max = Max[Abs[dF04]];
dynerr1 = Transpose[{dF01[[;; , 1]], 20 Log10[Abs[F0[#[[1]], 1, 5] - #[[2]]]/dF01max] & /@dF01}];
dynerr4 = Transpose[{dF04[[;; , 1]], 20 Log10[Abs[F0[#[[1]], 4, 5] - #[[2]]]/dF04max] & /@dF04}];
GraphicsGrid[{
 {Show[
  Plot[F0[u, 1, 5], {u, 0, 20}, 
    PlotRange -> {-0.005, 0.05}, 
    Exclusions -> None, 
    Frame -> True, 
    Axes -> False, 
    FrameLabel -> {"u", "\!\(\*SubscriptBox[\(f\), \(2\)]\)[\[Nu]]"}, 
    PlotLabel -> "HT, p\[Equal]1"],
  ListPlot[dF01]], 
  ListLinePlot[dynerr1, 
   Frame -> True, 
   Axes -> False, 
   FrameLabel -> {"\[Nu]", "dB"}, 
   PlotLabel -> "e(\[Nu]), p\[Equal]1", 
   PlotRange -> {{0, 20}, {-180, 0}}]}, 
{Show[
  Plot[F0[u, 4, 5], {u, 0, 20}, 
   PlotRange -> {-0.03, 0.03}, 
   Exclusions -> None, 
   Frame -> True, 
   Axes -> False, 
   FrameLabel -> {"\[Nu]", "\!\(\*SubscriptBox[\(f\), \(2\)]\)[\[Nu]]"}, 
   PlotLabel -> "HT, p\[Equal]4"], ListPlot[dF04]], 
  ListLinePlot[dynerr4, 
   Frame -> True, 
   Axes -> False, 
   FrameLabel -> {"\[Nu]", "dB"}, 
   PlotLabel -> "e(\[Nu]), p\[Equal]4", 
   PlotRange -> {{0, 20}, {-180, 0}}]}}, 
   ImageSize -> Large]

Which yields a comparison between the analytical Hankel transform of the sinc function and the discrete version including the dynamic error of the discrete transform as outlined in Fig. 1 of the reference. .

This discrete Transform is perfectly suited for the simulation of the wave propagation through radial symmetric apertures. One can get very nice results of the propagation of Bessel beams as depicted in Fig. 3 of Sicairos: