Discrete Hankel transform

Discrete Hankel transform is the numerical analogy of the Bessel transform.

Under construction, this article is not yet finished.

Below, the description by GNU is copypasted with minimal torification.

GNU version

The discrete Hankel transform acts on a vector of sampled data, where the samples are assumed to have been taken at points related to the zeroes of a Bessel function of fixed order; compare this to the case of the discrete Fourier transform, where samples are taken at points related to the zeroes of the sine or cosine function.

Specifically, let \(f(t)\) be a function on the unit interval and \(j_{\nu,m}\) the \(m\)-th zero of the Bessel function \(J_\nu(x)\). Then the finite \(\nu\)-Hankel transform of \(f(t)\) is defined to be the set of numbers \(g_m\) given by,

\( \displaystyle g_m = \int_0^1 t~\mathrm d t ~ J_\nu(j_{\nu,m}t)~ f(t), \)

so that,

\( \displaystyle f(t) = \sum_{m=1}^\infty (2 J_\nu(j_{\nu,m}t) / J_{\nu+1}(j_{\nu,m})^2) g_m. \)

Suppose that \(f\) is band-limited in the sense that \(~g_m\!=\!0~\) for \(~m\! >\! M~\). Then we have the following fundamental sampling theorem.

\( \displaystyle g_m = (2 / j_{\nu,M}^2) \sum_{k=1}^{M-1} f(j_{\nu,k}/j_{\nu,M}) (J_\nu(j_{\nu,m} j_{\nu,k} / j_{\nu,M}) / J_{\nu+1}(j_{\nu,k})^2). \)

It is this discrete expression which defines the discrete Hankel transform. The kernel in the summation above defines the matrix of the \(\nu\)-Hankel transform of size \(M\!-\!1\). The coefficients of this matrix, being dependent on \(\nu\) and \(M\), must be precomputed and stored; the gsl_dht object encapsulates this data. The allocation function gsl_dht_alloc returns a gsl_dht object which must be properly initialized with gsl_dht_init before it can be used to perform transforms on data sample vectors, for fixed \(\nu\) and \(M\), using the gsl_dht_apply function. The implementation allows a scaling of the fundamental interval, for convenience, so that one can assume the function is defined on the interval [0,X], rather than the unit interval.

Notice that by assumption \(f(t)\) vanishes at the endpoints of the interval, consistent with the inversion formula and the sampling formula given above. Therefore, this transform corresponds to an orthogonal expansion in eigenfunctions of the Dirichlet problem for the Bessel differential equation.

Vulgarizarion

The final equation by GNU can be written as follows:

\( \displaystyle g_m = \frac{2}{ j_{\nu,M}^2} \sum_{k=1}^{M-1} ~ f_k ~ \Big(J_\nu(j_{\nu,m} j_{\nu,k} / j_{\nu,M}) / J_{\nu+1}(j_{\nu,k})^2\Big) \)

where

\(\displaystyle f_k=f\Big( \frac{j_{\nu,k}}{j_{\nu,M}} \Big)\)

Introduce the transfer matrix \(T\) with following elements:

\( \displaystyle T_{m,k}= \frac{2}{ j_{\nu,M}^2} \Big(J_\nu(j_{\nu,m} j_{\nu,k} / j_{\nu,M}) / J_{\nu+1}(j_{\nu,k})^2\Big) \) \(\displaystyle = \frac{ 2~J_\nu\!\Big( j_{\nu,m}~ j_{\nu,k} / j_{\nu,M} \Big) }{ J_{\nu+1}\!(j_{\nu,k})^2 ~ j_{\nu,M}^2} \)

This kernel is not symmetric; the careful test is necessary. Since the Bessel transform is its own inverse the sought transformation matrix must be symmetric .

Formulas by Manuel

Different (but perhaps in some sense e`uivalent) form had been suggested at Arizona ; the matrix suggested is explicitly symmetric:

\( \displaystyle T_{m,n} = \frac { 2 ~ J_p( \alpha_{pn} \alpha_{pn} /S ) } { |J_{p+1}(\alpha_{pn})| ~ |J_{p+1}(\alpha_{pm})| ~ S }\)

\(S=2 \pi R V\)

\(F_m=f(\alpha_{pm}/(2 \pi R)) |J_{p+1}(\alpha_{pm})|^{-1}\)

The sense of the scaling constants can be revealed after the implementation.

Modified Hankel and its implementation

The modified transform (with scaled argument) can be written as follows:

\(f(x)=\frac{1}{2\pi} \int_{0}^{\infty} J_{\nu}(2 \pi x y)\, g(y) \, y \, \mathrm d y\)

The transform and \(g\) and its transform \(g\) are represented at the discrete set of \(N\) values of argument

\(\displaystyle X_n= \frac{ \alpha_{0,n} }{ \sqrt{2 \pi \alpha_{\nu,N+1}} }\)

for \(n=1 .. N\) , where \(\alpha_{\nu,n}\) is \(n\)th zero of \(\nu\)th Bessel.

Functions are represented with

\(\displaystyle F_n= f(X_n) \frac{\sqrt{\alpha_{\nu,N+1}}}{|J_{\nu+1}(\alpha_{\nu,n})|}\)

\(\displaystyle G_n= g(X_n) \frac{\sqrt{\alpha_{\nu,N+1}}}{|J_{\nu+1}(\alpha_{\nu,n})|}\)

The transform is approximated with

\(\displaystyle F_m \approx \sum_{n=1}^N T_{m,n} G_n\)

\(\displaystyle G_m \approx \sum_{n=1}^N T_{m,n} F_n\)

where

\(\displaystyle T_{m,n}=\frac {2 J_\nu(\alpha_{nu,n}\alpha_{nu,m} / \alpha_{\nu,N+1} } { |J_{nu+1}(\alpha_{\nu,n})|\, |J_{nu+1}(\alpha_{\nu,m})|\, \alpha_{\nu,N+1} } \)

The implementation with Mathematica for \(\nu=0\) and \(N=M=9\) and its test are copypasted below:

M=10
T= Table[Table[ (2. BesselJ[0, (BesselJZero[0, n] BesselJZero[0, m])/
       BesselJZero[0, M+1]])/(
       Abs[BesselJ[1, BesselJZero[0, n]]]
       Abs[BesselJ[1, BesselJZero[0, m]]]
       BesselJZero[0,M+1]),
{n,1,M}], {m,1,M}];
Id = IdentityMatrix[M];
Print[M]; TT = T.T;
di[M] = Norm[TT-Id]; Print[di[M]] ; Print[TableForm[Chop[TT, 10^(-7)]]

For this implementation, \(T^2 \approx I_M\), where \(I_M\) is identity matrix. The error of this approximation reduces as \(1/M^3\) and can be fitted with function \(\displaystyle d_M=\frac{1}{31000+7900M^2+600 M^3}\)

Application of the discrete Hankel tansform

The discrete Hankel tansform allows to describe "analytically" the propagation of waves with circular (axial) symmetry. Even if the \(N~\log N\) algorithms are not used, the explicit use of the symmetry saves an order of magnitude at the calculus of the distribution of field while it propagates in a free space or in an idealized waveguide. Therefore, the simple (independent on any specific installations) implementation should be considered as the important tool of the scientific research.