File Exchange

image thumbnail

Hankel Transform

version (2.9 KB) by Adam Wyatt
Routine to perform a QDHT with no limit on data size or transform order (other than memory constrain


Updated 04 Mar 2009

View License

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.

Comments and Ratings (9)

I want to do discrete Hankel transformation of the 2D image/ 2D matrix. How can I use this function for the 2D matrix?

anping Yu


how can we use it for the Hankel transformation of 2D image?

Ahmed Fasih

Answered my own question: but essentially, to produce Fig 1(c):

h = Hankel.hankel_matrix(4, 3, 256, 1e-13);
gamma = 5;
f1fun = @(r) sinc(2*gamma*r); % our sinc includes pi
r = h.J_roots / (2*pi*h.vmax);
F1 = f1fun(r) ./ h.J * h.vmax; F2 = h.T * F1;
f2 = F2 .* h.J / h.vmax;
v = h.J_roots / (2*pi*h.rmax);
figure; plot(v, f2)

(See the link for properly-commented version! Don't write code like this!)

Ahmed Fasih

Can you show how to use this function to reproduce the sinc plot, Figure 1 in Guizar-Sicairos & Gutierrez-Vega?


Thanks for this submission.
I have a question: can one uses this code for non normalized parameters?

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:


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.


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.


This is a nice submission.


Slight changes to the hankel structure to simplify usage.

Minor modifications.

MATLAB Release Compatibility
Created with R2008b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Inspired by: Integer order Hankel transform