Hankel

Hankel implements the quasi-discrete Hankel transform (QDHT) as introduced in H. F. Johnson, Comput. Phys. Commun., 43, 2 (1987), laid out for 0-order in L. Yu, et al., Optics Letters 23 (1998), and extended to higher orders in M. Guizar-Sicairos and J. C. Gutiérrez-Vega, JOSA A 21, 53 (2004). We generalized the cylindrical QDHT in the above publications to its (hyper)spherical form. The forward $p^{\mathrm{th}}$-order (hyper)spherical Hankel transform of the radially symmetric function $f(r)$ is defined as

\[\tilde{f}(k) = c_n^{-1} \int_0^\infty f(r) j_p^n(kr) r^n\,\mathrm{d}r\,,\]

where $c_n$ is a constant (see sphbesselj_scale), and $j_p^n(x)$ is the $p^{\mathrm{th}}$-order (hyper)spherical Bessel function of the first kind with spherical dimension $n$ (see sphbesselj). Correspondingly, the inverse transform is

\[f(r) = c_n^{-1} \int_0^\infty \tilde{f}(k) j_p^n(kr) k^n\,\mathrm{d}k\,.\]

Note that here, $k$ includes the factor of $2\pi$, i.e. it is spatial angular frequency.

The remainder of the documentation considers without loss of generalization only the cylindrical case ($n=1$, with $j_p^n(x) = J_p(x)$).

The QDHT approximates these transforms under the assumption that $f(r) = 0$ for $r > R$ where $R$ is the aperture size. By expanding $f(r)$ as a Fourier-Bessel series and choosing to sample $f(r)$ at coordinates $r_n = j_nR/j_{N+1}$, where $j_n$ is the $n^{\mathrm{th}}$ zero of $J_p(x)$ and $N$ is the number of samples, the Hankel transform turns into a matrix-vector multiplication.

Hankel follows the AbstractFFTs approach of planning a transform in advance by creating a QDHT object, which can then be used to transform an array using multiplication (* or mul!) and to transform back by left-division (\ or ldiv!). In contrast to AbstractFFTs, however, pre-planning a transform is required since the transform only works properly if the sampling points are chosen for a particular combination of sample number $N$ and aperture size $R$. The workflow to transform a function $f(r)$ is therefore as follows:

j01 = besselj_zero(0, 1)
R = 1
N = 8
q = QDHT(0, R, N)
f(r) = besselj(0, r*j01/R)
fr = f.(q.r)
fk = q * fr # 0th-order QDHT => fk should have only the first entry non-zero

8-element Array{Float64,1}:
  0.13475706165848417
  4.76372956333634e-10
 -6.088187982006206e-10
  7.288009860206225e-10
 -8.36749654060558e-10
  9.16271604004363e-10
 -9.187249659023021e-10
  7.19563708969384e-10

Normalisation

The transform as implemented here is unitary, i.e.

q \ (q * fr) ≈ fr

true

The transform satisfies Parseval's Theorem. To calculate the $L^2$ norm (e.g. energy in a signal) with correct scaling, use integrateR in real ($r$) space and integrateK in reciprocal ($k$) space. For an explanation of how these functions work, see the Derivations page.

On-axis and symmetric samples

The QDHT does not contain a sample on axis, i.e. at $r=0$, but for the $0^{\mathrm{th}}$-order QDHT, it can be obtained using onaxis, which takes the transformed array as its input. This is because the on-axis sample is obtained from

\[f(r=0) = \int_0^\infty \tilde{f}(k) J_0(0) k\,\mathrm{d}k\,.\]

Note here why this does not work for other values of $p$: $J_p(0) = 0 \; \forall \; p>0$.

Given an array A, the convenience method symmetric produces the symmetric array, i.e. given A, sampled at $[r₁, r₂, r₃, ...]$, it generates the array sampled at $[...-r₃, -r₂, -r₁, 0, r₁, r₂, r₃...]$. This works also for higher-dimensional arrays. The corresponding spatial coordinates can be obtained with Rsymmetric.

Function Reference

Hankel.QDHT — Type

QDHT{p, n}(R, N; dim=1)
QDHT([p, ] R, N; dim=1)
QDHT(p, [n, ] R, N; dim=1)

p-th order quasi-discrete Hankel transform over aperture radius R with N samples which transforms along dimension dim. If not given, p defaults to 0, and n the spherical dimension defaults to 1 (cylindrical).

After:

[1] L. Yu, M. Huang, M. Chen, W. Chen, W. Huang, and Z. Zhu, Optics Letters 23 (1998)

[2] M. Guizar-Sicairos and J. C. Gutiérrez-Vega, JOSA A 21, 53 (2004)

[3] H. F. Johnson, Comput. Phys. Commun., 43, 2 (1987)

but with some alterations:

The transform matrix T is not the same as C/T defined in [1, 2] but is more like the form used in equation 14 of [3]. Instead of dividing by $J_1(α_{pn}) J_1(α_{pm})$ we divide by $J_1(α_{pn})^2$. This cancels out the factor between $f$ and $F$ so we do not have to mutltiply (divide) by $J_1(α_{pn})$ ($J_1(α_{pm})$) before and after applying the transform matrix.

Follows AbstractFFT approach of applying fwd and inv transform with mul and ldiv.

To calculate radial integrals of functions sampled using QDHT, use integrateR and integrateK.

The type of the coefficients is inferred from the type of R (but is promoted to be at least Float), so for arbitrary precision use QDHT([p, ] BigFloat(R), ...).