{ } Raw JSON

bundles / scipy 1.17.1 / scipy / fft / _fftlog / fht

_Function

`scipy.fft._fftlog:fht`

source: /scipy/fft/_fftlog.py :18

Signature

def fht ( a , dln , mu , offset = 0.0 , bias = 0.0 )

Summary

Compute the fast Hankel transform.

Extended Summary

Computes the discrete Hankel transform of a logarithmically spaced periodic sequence using the FFTLog algorithm ^[1], ^[2].

Parameters

a : array_like (..., n): Real periodic input array, uniformly logarithmically spaced. For multidimensional input, the transform is performed over the last axis.
dln : float: Uniform logarithmic spacing of the input array.
mu : float: Order of the Hankel transform, any positive or negative real number.
offset : float, optional: Offset of the uniform logarithmic spacing of the output array.
bias : float, optional: Exponent of power law bias, any positive or negative real number.

Returns

A : array_like (..., n): The transformed output array, which is real, periodic, uniformly logarithmically spaced, and of the same shape as the input array.

Notes

This function computes a discrete version of the Hankel transform

A (k) = \int_{0}^{\infty} a (r) J_{μ} (k r) k d r,

where $J_{μ}$ is the Bessel function of order $μ$ . The index $μ$ may be any real number, positive or negative. Note that the numerical Hankel transform uses an integrand of $k d r$ , while the mathematical Hankel transform is commonly defined using $r d r$ .

The input array a is a periodic sequence of length $n$ , uniformly logarithmically spaced with spacing dln,

a_{j} = a (r_{j}), r_{j} = r_{c} exp [(j - j_{c}) dln]

centred about the point $r_{c}$ . Note that the central index $j_{c} = (n - 1) /2$ is half-integral if $n$ is even, so that $r_{c}$ falls between two input elements. Similarly, the output array A is a periodic sequence of length $n$ , also uniformly logarithmically spaced with spacing dln

A_{j} = A (k_{j}), k_{j} = k_{c} exp [(j - j_{c}) dln]

centred about the point $k_{c}$ .

The centre points $r_{c}$ and $k_{c}$ of the periodic intervals may be chosen arbitrarily, but it would be usual to choose the product $k_{c} r_{c} = k_{j} r_{n - 1 - j} = k_{n - 1 - j} r_{j}$ to be unity. This can be changed using the offset parameter, which controls the logarithmic offset $lo g (k_{c}) = offset - lo g (r_{c})$ of the output array. Choosing an optimal value for offset may reduce ringing of the discrete Hankel transform.

If the bias parameter is nonzero, this function computes a discrete version of the biased Hankel transform

A (k) = \int_{0}^{\infty} a_{q} (r) (k r)^{q} J_{μ} (k r) k d r

where $q$ is the value of bias, and a power law bias $a_{q} (r) = a (r) (k r)^{- q}$ is applied to the input sequence. Biasing the transform can help approximate the continuous transform of $a (r)$ if there is a value $q$ such that $a_{q} (r)$ is close to a periodic sequence, in which case the resulting $A (k)$ will be close to the continuous transform.

Array API Standard Support

fht has experimental support for Python Array API Standard compatible backends in addition to NumPy. Please consider testing these features by setting an environment variable SCIPY_ARRAY_API=1 and providing CuPy, PyTorch, JAX, or Dask arrays as array arguments. The following combinations of backend and device (or other capability) are supported.

====================  ====================  ====================
Library               CPU                   GPU
====================  ====================  ====================
NumPy                 ✅                     n/a                 
CuPy                  n/a                   ✅                   
PyTorch               ✅                     ✅                   
JAX                   ✅                     ✅                   
Dask                  ⚠️ computes graph     n/a                 
====================  ====================  ====================

See dev-arrayapi for more information.

Examples

This example is the adapted version of ``fftlogtest.f`` which is provided in [2]_. It evaluates the integral .. math:: \int^\infty_0 r^{\mu+1} \exp(-r^2/2) J_\mu(kr) k dr = k^{\mu+1} \exp(-k^2/2) .

import numpy as np
from scipy import fft
import matplotlib.pyplot as plt

✓

Parameters for the transform.

mu = 0.0                     # Order mu of Bessel function
r = np.logspace(-7, 1, 128)  # Input evaluation points
dln = np.log(r[1]/r[0])      # Step size
offset = fft.fhtoffset(dln, initial=-6*np.log(10), mu=mu)
k = np.exp(offset)/r[::-1]   # Output evaluation points

✓

Define the analytical function.

def f(x, mu):
    """Analytical function: x^(mu+1) exp(-x^2/2)."""
    return x**(mu + 1)*np.exp(-x**2/2)

✓

Evaluate the function at ``r`` and compute the corresponding values at ``k`` using FFTLog.

a_r = f(r, mu)
fht = fft.fht(a_r, dln, mu=mu, offset=offset)

✓

For this example we can actually compute the analytical response (which in this case is the same as the input function) for comparison and compute the relative error.

a_k = f(k, mu)
rel_err = abs((fht-a_k)/a_k)

✓

Plot the result.

figargs = {'sharex': True, 'sharey': True, 'constrained_layout': True}
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), **figargs)

✓

ax1.set_title(r'$r^{\mu+1}\ \exp(-r^2/2)$')
ax1.loglog(r, a_r, 'k', lw=2)
ax1.set_xlabel('r')
ax2.set_title(r'$k^{\mu+1} \exp(-k^2/2)$')
ax2.loglog(k, a_k, 'k', lw=2, label='Analytical')
ax2.loglog(k, fht, 'C3--', lw=2, label='FFTLog')
ax2.set_xlabel('k')
ax2.legend(loc=3, framealpha=1)
ax2.set_ylim([1e-10, 1e1])

✗

ax2b = ax2.twinx()

✓

ax2b.loglog(k, rel_err, 'C0', label='Rel. Error (-)')
ax2b.set_ylabel('Rel. Error (-)', color='C0')

✗

ax2b.tick_params(axis='y', labelcolor='C0')

✓

ax2b.legend(loc=4, framealpha=1)
ax2b.set_ylim([1e-9, 1e-3])

✗

plt.show()

✓

`scipy.fft._fftlog:fht`

Signature

Summary

Extended Summary

Parameters

Returns

Notes

Examples

See also

Aliases

Referenced by

This package