2

I'm trying to find any existing implementation for Hankel Transform in Python (actually i'm more into symmetric fourier transform of two 2d radially symmetric functions but it can be easily reduced to hankel transform).

I do know about hankel python module, but it requires lambda function for input whereas I have only 1d-array.

Any thoughts?

Anton Savostyanov
  • 53
  • 1
  • 6
  • 1
    The GNU Scientific Library (GSL) includes [functions for computing the discrete Hankel transform](https://www.gnu.org/software/gsl/manual/html_node/Discrete-Hankel-Transforms.html#Discrete-Hankel-Transforms). [PyGSL](http://pygsl.sourceforge.net/) provides a Python wrapper of GSL. It has been awhile since I tried installing it, so I don't know how challenging that will be. – Warren Weckesser Apr 06 '17 at 16:02

4 Answers4

7

I'm the author of hankel. While I would not recommend to use my code in this case (since as you mention, it requires a callable input function, and its purpose is to accurately compute integrals, not to do a DHT), I will say that it is possible.

All you need to do is interpolate your input 1D array. How to do this is up to you, but typically something like the following works pretty well:

from scipy.interpolate import InterpolatedUnivariateSpline as Spline
import numpy as np

x, y = # Import/create data vectors

# Do this if y is both negative and positive
fnc = Spline(x,y, k=1) #I usually choose k=1 in case anything gets extrapolated.

# Otherwise do this
spl = Spline(np.log(x), np.log(y), k=1)
fnc = lambda x : np.exp(spl(np.log(x)))

# Continue as normal with hankel.transform(fnc, kvec)

The big issue with doing this is in choosing the parameters N and h such that the transform is well approximated for all values of k in kvec. If kvec spans a wide dynamic range, then hankel is very inefficient as it uses the same underlying arrays (of length N) for each k in a transform, meaning the most difficult k sets the performance level.

So again, in short I wouldn't recommend hankel, but if you can't find anything else, it will still work ;-)

StevenMurray
  • 742
  • 2
  • 7
  • 18
5

A friend of mine was looking for a Hankel transform in Python and after deciding there was nothing available - partly on the basis of the answers to this question - approached me to share mine.

I have now tidied, documented and released the code as PyHank (docs).

This package performs Quasi-discrete Hankel transforms, as requested by the original question (if a few years too late!) but hopefully this answer will direct future users who find this question (like my friend) to the right place.

EdR
  • 503
  • 4
  • 13
  • That is a very nice and well-documented package and it seems to do exactly what the question is about. (And I was looking for something exactly like this for my project, so thank you!) – PoorYorick Sep 10 '20 at 15:30
1

If you read the documentation for Hankel all the way through the Transforms section you will see that to perform the transformation you call hankel.transform(function, array, ret_err=bool) so you just need the function for whatever form of the transformation you require. I believe there is a list of transformation functions in the Wikipedia entry for Hankel Transformation.

Grr
  • 15,553
  • 7
  • 65
  • 85
  • I should probably state it more clearly: I've read the docs and the point is it requires __callable__ function for input in `hankel.transform`, but I have 1d-array instead – Anton Savostyanov Apr 06 '17 at 23:01
  • Are you saying your function is represented as a 1d array. I was under the impression that you needed both a callable and an array for the `hankel.tranform` method. – Grr Apr 07 '17 at 17:37
  • 1
    What I need is `hankel.transform(array, array, ret_err)` – Anton Savostyanov Apr 07 '17 at 17:40
0

As of 2023 there is scipy implementation of discrete Hankel transform.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fht.html

I am also interested in computing inverse 2D Fourier of radially symmetric function. I am quite not sure if the worse is sampling error when I try to sample it on the 2D grid and compute inverse 2D Fourier directly or if I compute Hankel transform and interpolate on logarithmically spaced grid.

VojtaK
  • 483
  • 4
  • 13