File Exchange

## Hankel Transform

version 1.1.0.0 (2.9 KB) by Adam Wyatt

### Adam Wyatt (view profile)

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

Updated 04 Mar 2009

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.

### Cite As

Adam Wyatt (2020). Hankel Transform (https://www.mathworks.com/matlabcentral/fileexchange/15623-hankel-transform), MATLAB Central File Exchange. Retrieved .

Garima Nagar

### Garima Nagar (view profile)

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

anping Yu

RAVI KUMAR

### RAVI KUMAR (view profile)

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

Ahmed Fasih

### Ahmed Fasih (view profile)

Answered my own question: https://gist.github.com/fasiha/16f796be6a55cc7dae98 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

### Ahmed Fasih (view profile)

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

Ismail

### Ismail (view profile)

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

### Adam Wyatt (view profile)

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.

Roy

### Roy (view profile)

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.

Kevin

### Kevin (view profile)

This is a nice submission.