我有一个实值函数G(t),我想在对数均匀间隔的时间内插值:x1c 0d1x。这是为了获得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
包含ts
中所有时间的函数G(t)
的值,G(t)
是上面在文件中定义的实值函数(我相信它的确切定义并不重要):ValueError: A value in x_new is below the interpolation range.
经过快速的研究,似乎每当x(或在本例中t)中传递的数据不是线性均匀间隔时,就会发生这种错误,所以还没有什么异常。
然而,当我试图通过以下方式定义函数log_interp1d
来应用指定的方法there时:
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
替换为log_interp1d
,会发生完全相同的错误,这让我很困惑。
我已经查找了其他链接,但似乎没有提供可靠的建议,以解决这个问题。我仍然可以使用线性均匀间隔的时间来做这个插值,但这将导致一个非常不均匀的数据间隔一旦在对数尺度,我想避免这种情况。我如何插值我的数据,同时保持对数均匀间隔的时间(因此频率)值?
先谢了!
1条答案
按热度按时间vojdkbi01#
错误与日志或其他数据间隔无关。它只是锡上写的:你正在尝试计算插值器在插值范围(x.min(),x.max())之外。
为了避免错误,将bounds_error=False传递给interp1d。然而,外推的质量可能会有很大的变化,所以YMMV。我建议重新评估这种方法,并弄清楚 * 为什么 * 你要外推,以及是否可以使用一些外部知识来代替。