scipy.fft.

fht#

scipy.fft.fht(a, dln, mu, offset=0.0, bias=0.0)[源代码]#

计算快速汉克尔变换。

使用 FFTLog 算法计算对数间隔周期序列的离散汉克尔变换 [1], [2]。

参数:

aarray_like (…, n): 实周期输入数组，均匀对数间隔。对于多维输入，变换在最后一个轴上执行。
dlnfloat: 输入数组的均匀对数间距。
mufloat: 汉克尔变换的阶数，任何正或负实数。
offsetfloat, 可选: 输出数组的均匀对数间距的偏移量。
biasfloat, 可选: 幂律偏差的指数，任何正或负实数。

返回:

Aarray_like (…, n): 变换后的输出数组，它是实数的、周期的、均匀对数间隔的，并且与输入数组的形状相同。

另请参阅

ifht: fht 的逆。
fhtoffset: 为 fht 返回最佳偏移量。

说明

此函数计算汉克尔变换的离散版本

\[A(k) = \int_{0}^{\infty} \! a(r) \, J_\mu(kr) \, k \, dr \;,\]

其中 \(J_\mu\) 是阶数为 \(\mu\) 的贝塞尔函数。索引 \(\mu\) 可以是任何实数，正数或负数。请注意，数值汉克尔变换使用 \(k \, dr\) 的被积函数，而数学汉克尔变换通常使用 \(r \, dr\) 定义。

输入数组 a 是长度为 \(n\) 的周期序列，均匀对数间隔，间距为 dln，

\[a_j = a(r_j) \;, \quad r_j = r_c \exp[(j-j_c) \, \mathtt{dln}]\]

以 \(r_c\) 点为中心。请注意，如果 \(n\) 为偶数，则中心索引 \(j_c = (n-1)/2\) 是半整数，因此 \(r_c\) 落在两个输入元素之间。类似地，输出数组 A 是长度为 \(n\) 的周期序列，也均匀对数间隔，间距为 dln

\[A_j = A(k_j) \;, \quad k_j = k_c \exp[(j-j_c) \, \mathtt{dln}]\]

以 \(k_c\) 点为中心。

周期间隔的中心点 \(r_c\) 和 \(k_c\) 可以任意选择，但通常选择乘积 \(k_c r_c = k_j r_{n-1-j} = k_{n-1-j} r_j\) 为 1。可以使用 offset 参数更改此设置，该参数控制输出数组的对数偏移量 \(\log(k_c) = \mathtt{offset} - \log(r_c)\)。为 offset 选择最佳值可以减少离散汉克尔变换的振铃。

如果 bias 参数非零，则此函数计算有偏汉克尔变换的离散版本

\[A(k) = \int_{0}^{\infty} \! a_q(r) \, (kr)^q \, J_\mu(kr) \, k \, dr\]

其中 \(q\) 是 bias 的值，并且将幂律偏差 \(a_q(r) = a(r) \, (kr)^{-q}\) 应用于输入序列。如果存在值 \(q\) 使得 \(a_q(r)\) 接近周期序列，则偏差变换有助于近似 \(a(r)\) 的连续变换，在这种情况下，结果 \(A(k)\) 将接近连续变换。

参考文献

[1]

Talman J. D., 1978, J. Comp. Phys., 29, 35

[2] (1,2)

Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191)

示例

此示例是 fftlogtest.f 的改编版本，该版本在 [2] 中提供。它计算积分

\[\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

变换的参数。

>>> 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

定义解析函数。

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

在 r 处评估函数，并使用 FFTLog 计算 k 处的相应值。

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

对于此示例，我们实际上可以计算分析响应（在这种情况下，它与输入函数相同）以进行比较并计算相对误差。

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

绘制结果。

>>> 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()