# Computation of Hankel Transform using FFT (Fourier)

To address circular symmetric cases of 2D Fourier Transformations the so called Hankel Transform can be applied (for a detailed derivation of the relation between the 2D Fourier transform and the 1D 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), which 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? –  Ｊ. Ｍ. 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? –  Ｊ. Ｍ. 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]], 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:

-
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
@Pavel I have converted your "answer" into a comment. Answer posts are intended to be used strictly for that. The designers of StackExchange want you to build some "reputation" points before you get the privilege to comment. I'm not sure I agree with that policy but it is beyond my power to change it. I encourage you to participate by asking and answering questions and soon this will be a non-issue. –  Mr.Wizard Aug 2 '13 at 19:40
@Pavel With only 20 points you will be able to chat which is a better way to discuss things. –  Mr.Wizard Aug 2 '13 at 19:42