我有一个实值函数 G(t),我想在对数均匀间隔时间上进行插值:。这是为了获得 G 的分量作为频率 而不是时间的函数的近似表达式。因此,我这样定义了一个时间列表:
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
exponents = np.linspace(-2, 6, 81, endpoint=True) #81 values
ts = [10**(xp) for xp in exponents] #times (s)
我在那里使用了
numpy.linspace
,因为numpy.logspace
似乎产生了错误。最终结果应该看起来像这样(不同颜色的曲线分别根据另一个与此处无关的参数获得):
为此,我尝试使用
scipy.interpolate.interp1d
,如下面的函数所示:
def G_pr_sec(ts, Goft):
"""
Returns G_prime and G_sec at frequencies 1/ts[i] by interpolating the array Goft
Parameters
----------
ts : 1D array
Array of times
Goft : 1D array
Values of G(t)
Returns
-------
G_pr : 1D array
Array of G_prime(1/t)
G_sec : 1D array
Array of G_second(1/t)
"""
G2 = interp1d(ts, Goft)
G_pr = []
G_sec = []
for t in ts:
res_sec = 0
res_sec += -0.47*(G2(2*t) - G2(4*t)) + 1.674*(G2(t) - G2(2*t)) + 0.198*(G2(t/2) - G2(t))
res_sec += 0.62*(G2(t/4) - G2(t/2)) + 0.012*(G2(t/8) - G2(t/4)) + 0.172*(G2(t/16) - G2(t/8))
res_sec += 0.0433*(G2(t/64) - G2(t/32)) + 0.0108*(G2(t/256) - G2(t/128))
res_sec += (0.0108/4)*(G2(t/1024) - G2(t/512)) + (0.0108/16)*(G2(t/4096) - G2(t/2048))
res_sec += (0.0108/64)*(G2(t/(64*256)) - G2(t/(32*256))) + (0.0108/256)*(G2(t/(256**2)) - G2(t/(128*256)))
G_sec.append(res_sec)
res_pr = 0
res_pr += G2(t) -0.142*(G2(4*t) - G2(8*t)) + 0.717*(G2(2*t) - G2(4*t)) + 0.046*(G2(t) - G2(2*t))
res_pr += 0.099*(G2(t/2) - G2(t)) + 0.103*(G2(t/4) - G2(t/2)) + 0.001*(G2(t/8) - G2(t/4))
res_pr += 0.00716*(G2(t/16) - G2(t/8)) + 0.000451*(G2(t/64) - G2(t/32))
G_pr.append(res_pr)
return [G_pr, G_sec]
但是当使用包含函数
G_pr_sec(ts, Gs)
的值的Gs
调用G(t)
时立即出现以下错误,ts
是文件中上面定义的实值函数定义在那里并不重要,我相信):G(t)
经过快速研究,似乎只要传入 x(或当前情况下的 t)的数据不是线性均匀分布的,就会发生此错误,所以还没有什么异常。
但是,当我尝试通过以下方式定义函数
ValueError: A value in x_new is below the interpolation range.
来应用指定there的方法时:
log_interp1d
在前面的函数体中用
import scipy as sp
import scipy.interpolate
def log_interp1d(xx, yy, kind='linear'):
logx = np.log10(xx)
logy = np.log10(yy)
lin_interp = sp.interpolate.interp1d(logx, logy, kind=kind)
log_interp = lambda zz: np.power(10.0, lin_interp(np.log10(zz)))
return log_interp
替换
interp1d
,发生了完全相同的错误,这让我很困惑。我查看了其他链接,但似乎没有一个提供可靠的建议来解决这个问题。我仍然可以使用线性均匀间隔的时间来进行这种插值,但这会导致数据在对数尺度上的间隔非常不均匀,我想避免这种情况。我如何在保持对数均匀间隔的时间(以及频率)值的同时插值我的数据?
提前致谢!
为避免错误,将 bounds_error=False 传递给 interp1d。然而,外推的质量可能千差万别,所以 YMMV。我建议重新评估该方法并弄清楚
为什么你要推断,以及你是否可以使用一些外部知识。