Hankel Transform

Routine to perform a QDHT with no limit on data size or transform order (other than memory constrain

Description

This file contains 2 functions:
bessel_zeros - finds zeros of a bessel function
hankel_matrix - generates a structure of data to use for Hankel transforms

The algorithm is taken from:
M. Guizar-Sicairos and J. C. Gutierrez-Vega, Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,
J. Opt. Soc. Am. A 21, 53-58 (2004)

An example of beam propagation with circular symmetry has been written in the hankel_matrix help:
Type "help hankel_matrix" or "doc hankel_matrix" in Matlab.

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Integer order Hankel transform

• MATLAB release: MATLAB 7.7 (R2008b)

[Kevin]

This is a nice submission.

[Roy]

The example given inside hankel_matrix.m is incomplete and misguiding. According to it, one can define the Hankel and inverse-Hankel transform as follows:

mat_H = hankel_matrix(0, rmax, samples);
ht = @(f) mat_H.T * (f./mat_H.J);
iht = @(F) (mat_H.T * F) .* mat_H.J;

Recall that the Hankel transform of a Gaussian is a Gaussian:
H{e(-br^2)} = 1/2b * e(-k^2/4b)

However, trying to match this analytic result using the above transform fails. This is because the scaling vector J is defined as Jp1; a possible solution would be to define two vectors, Jp1/rmax and Jp1/vmax, and multiply f before and after the transform, and conversely in the inverse transform.
There is also a factor of 2*pi, but this depends on your formalism of choice.

[Adam Wyatt]

The example given in the code does not scale the amplitudes of the angular distribution because it is not necessary (since it is only necessary to multiply by a phase factor and so the scaling required will be undone again before transforming back). However, strictly speaking this should be done as follows:

H = hankel_matrix(0, rmax, samples);
HT = @(f) (H.T*(f./H.JR)).*H.JV;
IHT = @(F) (H.T*(F./H.JV)).*H.JR;

Note that two scaling vectors are supplied:

H.JR
H.JV

Therefore I agree that the example is misguiding from a physical point of view, but not computationally, but it is not incomplete.

Note that there is a difference between the Hankel transform (HT) and the quasi discrete Hankel transform (qDHT). This code provides a way of performing the qDHT and the scaling vectors to turn it into an accurate estimation of the HT.

Indeed there is a factor of 2pi, but as stated this does depend on the definition of the HT used - so long as f = IHT(HT(f)), its all good. So using the definition consistent with the code:

HT{exp(-br^2)} = pi/b * exp{-k^2/(2*b)];

Further details can be found in the references given in the code.