Hankel

Hankel implements the $0^{\mathrm{th}}$-order quasi-discrete Hankel transform (QDHT) as first laid out in L. Yu, et al., Optics Letters 23 (1998). The forward $0^{\mathrm{th}}$-order Hankel transform of the radially symmetric function $f(r)$ is defined as

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

where $J_0(x)$ is the $0^{\mathrm{th}}$-order Bessel function of the first kind. Correspondingly, the inverse transform is

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

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

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_0(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(R, N)
f(r) = besselj(0, r*j01/R)
fr = f.(q.r)
fk = q * fr # fk should have only the first entry non-zero

8-element Array{Float64,1}:
  0.13475706165848392
  4.763728629518283e-10
 -6.08818750247802e-10
  7.28800829543002e-10
 -8.367496368984974e-10
  9.162716080424964e-10
 -9.187255203883264e-10
  7.195637162884995e-10

Normalisation

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

q \ (q * fr) ≈ fr

true

To avoid unnecessary multiplications, the transform matrix ($C$ in Yu et al., here it's called T) is altered slightly (see QDHT), and as a consequence the transformation itself does not conserve energy (i.e. satisfy 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.

On-axis and symmetric samples

The QDHT does not contain a sample on axis, i.e. at $r=0$, but it can be obtained.

Function Reference

Hankel.QDHT — Type

QDHT(R, N; dim=1)

Quasi-discrete Hankel transform over aperture radius R with N samples which transforms along dimension dim

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)

but with some alterations:

The transform matrix T is not the same as C/T defined in [1, 2]. Instead of dividing by J₁(αₚₙ)J₁(αₚₘ) we divide by J₁(αₚₙ)^2. This cancels out the factor between f and F so we do not have to mutltiply (divide) by J₁(αₚₙ) (J₁(αₚₘ)) before and after applying the transform matrix. This means T is not symmetric, and does not conserve energy. To conserve energy, use integrateR and integrateK.

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

source

LinearAlgebra.mul! — Method

mul!(Y, Q::QDHT, A)

Calculate the forward quasi-discrete Hankel transform of array A using the QDHT Q and store the result in Y.

Examples

julia> q = QDHT(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R)); Y = similar(A);
julia> mul!(Y, q, A)
8-element Array{Float64,1}:
  4.326937831591551e-6
  2.3341589529175126e-6
  7.689558743828849e-7
  1.546419420523699e-7
  1.8999259906096856e-8
  1.4159642663129888e-9
  7.013670190083954e-11
 -6.07681871673291e-13

source

LinearAlgebra.ldiv! — Method

ldiv!(Y, Q::QDHT, A)

Calculate the inverse quasi-discrete Hankel transform of array A using the QDHT Q and store the result in Y.

Examples

julia> q = QDHT(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R)); Y = similar(A);
julia> mul!(Y, q, A);
julia> YY = similar(Y); ldiv!(YY, q, Y);
julia> YY ≈ A
true

source

Base.:* — Method

*(Q::QDHT, A)

Calculate the forward quasi-discrete Hankel transform of array A using the QDHT Q.

Examples

julia> q = QDHT(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R));
julia> q*A
8-element Array{Float64,1}:
  4.326937831591551e-6
  2.3341589529175126e-6
  7.689558743828849e-7
  1.546419420523699e-7
  1.8999259906096856e-8
  1.4159642663129888e-9
  7.013670190083954e-11
 -6.07681871673291e-13

source

Base.:\ — Method

\(Q::QDHT, A)

Calculate the inverse quasi-discrete Hankel transform of array A using the QDHT Q.

Examples

julia> q = QDHT(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R));
julia> Ak = q*A;
julia> q \ Ak ≈ A
true

source

Hankel.integrateR — Method

integrateR(A, Q::QDHT; dim=1)

Radial integral of A, over the aperture of Q in real space.

Assuming A contains samples of a function f(r) at sample points Q.r, then integrateR(A, Q) approximates ∫f(r)r dr from r=0 to r=∞.

Note

integrateR and integrateK fulfill Parseval's theorem, i.e. for some array A, integrateR(abs2.(A), q) and integrateK(abs2.(q*A), q) are equal, but integrateR(A, q) and integrateK(q*A, q) are not equal.

Examples

julia> q = QDHT(10, 128); A = exp.(-q.r.^2/2);
julia> integrateR(abs2.(A), q) ≈ 0.5 # analytical solution of ∫exp(-r²)r dr from 0 to ∞
true

source

Hankel.integrateK — Method

integrateK(Ak, Q::QDHT; dim=1)

Radial integral of A, over the aperture of Q in reciprocal space.

Assuming A contains samples of a function f(k) at sample points Q.k, then integrateR(A, Q) approximates ∫f(k)k dk from k=0 to k=∞.

Note

Examples

julia> q = QDHT(10, 128); A = exp.(-q.r.^2/2);
julia> integrateR(abs2.(A), q) ≈ 0.5 # analytical solution of ∫exp(-r²)r dr from 0 to ∞
true
julia> Ak = q*A;
julia> integrateK(abs2.(Ak), q) ≈ 0.5 # Same result
true

source

Hankel.onaxis — Method

onaxis(Ak, Q::QDHT; dim=Q.dim)

Calculate on-axis sample in space (i.e. at r=0) from transformed array Ak.

Examples

julia> q = QDHT(10, 128); A = exp.(-q.r.^2/2);
julia> onaxis(q*A, q) ≈ 1 # should be exp(0) = 1
true

source

Hankel.symmetric — Method

symmetric(A, Q::QDHT)

Create symmetric array from samples in A, including on-axis sample.

Given A, sampled at [r₁, r₂, r₃, ...], generates array sampled at [...-r₃, -r₂, -r₁, 0, r₁, r₂, r₃...]

Examples

julia> q = QDHT(10, 128); A = exp.(-q.r.^2);
julia> As = symmetric(A, q);
julia> size(As)
(257,)
julia> As[1:128] == A[128:-1:1]
true
julia> As[129] ≈ 1 # should be exp(0) = 1
true
julia> As[130:end] == A
true

source

Hankel.Rsymmetric — Method

Rsymmetric(Q::QDHT)

Create radial coordinate array to go along with symmetric(A, Q::QDHT).

Examples

julia> q = QDHT(10, 4);
julia> q.r
4-element Array{Float64,1}:
 1.6106347946239767
 3.697078919099734
 5.795844623798052
 7.8973942990196395
julia> Rsymmetric(q)
9-element Array{Float64,1}:
 -7.8973942990196395
 -5.795844623798052
 -3.697078919099734
 -1.6106347946239767
  0.0
  1.6106347946239767
  3.697078919099734
  5.795844623798052
  7.8973942990196395

source