Computation of Hankel Transform using FFT (Fourier)

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?

• Have you seen this, by any chance? – J. M. will be back soon Jun 1 '13 at 15:05
• Interesting reference, mentions multiple algorithms including the algorithm by Oppenheim, which I referenced above, thanks for pointing me on it. Unfortunately all of the algorithms would need some coding. I'll need to check which of the algorithms is most suited for my problem (radial symmetric diffraction at complex apertures). – Rainer Jun 1 '13 at 16:02
• "Unfortunately all of the algorithms would need some coding." - I know, I'm just saying you've plenty of options. Is using NIntegrate[] too slow for your needs? – J. M. will be back soon Jun 1 '13 at 16:08
• Sure, there are plenty of options. I'm just looking for something optimized readily available in Mathematica since for diffraction calculations I would need hundreds of iterations over the Hankel transform to calculate the optical near field of cylindrical apertures. NIntegrate[] is for sure one option, I'll just would like to do a review on available options before I reinvent the wheel. – Rainer Jun 1 '13 at 16:56

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, 5] - #[]]/dF01max] & /@dF01}];
dynerr4 = Transpose[{dF04[[;; , 1]], 20 Log10[Abs[F0[#[], 4, 5] - #[]]/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: • Might be a good idea to name the authors and the article title in the text, and use the DOI URL as the link target. The bare word "reference" is not very indicative of where the link will go, and you never know when opticsinfobase.org might change their URL scheme in the future. – Rahul Jun 21 '13 at 7:03
• Done, I was a bit sloppy in the reference you're right.. – Rainer Jun 21 '13 at 7:09
• Rainer, can you share anything regarding performance of your implementation of the Henkel transform? I.e., how much time does it take to calculate a field profile at one z coordinate. I am glad I found this post, as I have run in to the same problem myself. Thanks. P. – user8844 Aug 2 '13 at 14:30
• @xzczd damn you are right! on a related note: Rainer mentioned that Both functions can be specified as lists of x and y values too but then the points in which the functions are evaluated have to correspond to the zeros of the Bessel function, right? – yuriyi Nov 21 '18 at 8:14
• Responding to your comment on arbitrary discretized data: i am pretty sure the function to be transformed has to be evaluated at specific points (defined by the zeros of Bessel function). That is the idea behind the proposed algorithm: define Hankel tranform through Bessel series. Please consult the original paper. – yuriyi Nov 26 '18 at 16:59

protected by J. M. will be back soon♦Aug 18 at 3:26

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?