fht#
- scipy.fft.fht(a, dln, mu, offset=0.0, bias=0.0)[源代码][源代码]#
计算快速汉克尔变换。
使用 FFTLog 算法 [1]，[2] 计算对数间隔周期序列的离散 Hankel 变换。
- 参数:
- a类数组 (…, n)
真实的周期性输入数组，均匀对数间隔。对于多维输入，变换是在最后一个轴上进行的。
- dln浮动
输入数组的均匀对数间隔。
- mu浮动
Hankel 变换的阶数，可以是任何正或负的实数。
- 偏移量float, 可选
输出数组均匀对数间隔的偏移量。
- 偏见float, 可选
幂律偏差的指数，任何正或负的实数。
- 返回:
- A类数组 (…, n)
转换后的输出数组，它是实数、周期性的、均匀对数间隔的，并且与输入数组的形状相同。
注释
此函数计算Hankel变换的离散版本
\[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\) 中心。注意，中心索引 \(j_c = (n-1)/2\) 如果 \(n\) 是偶数则为半整数，因此 \(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}\)。对变换进行偏置可以帮助近似 \(a(r)\) 的连续变换，如果存在一个值 \(q\) 使得 \(a_q(r)\) 接近周期序列，在这种情况下，得到的 \(A(k)\) 将接近连续变换。
参考文献
[1]Talman J. D., 1978, J. Comp. Phys., 29, 35
示例
这个例子是
fftlogtest.f的改编版本，该版本在 [2] 中提供。它计算了积分\[\int^\infty_0 r^{\mu+1} \exp(-r^2/2) J_\mu(k, r) 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()
