如何使用scipy.interpolate.interp1d对对数等间距数据进行插值?

o2g1uqev  于 2023-03-30  发布在  其他
关注(0)|答案(1)|浏览(173)

我有一个实值函数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,会发生完全相同的错误,这让我很困惑。
我已经查找了其他链接,但似乎没有提供可靠的建议,以解决这个问题。我仍然可以使用线性均匀间隔的时间来做这个插值,但这将导致一个非常不均匀的数据间隔一旦在对数尺度,我想避免这种情况。我如何插值我的数据,同时保持对数均匀间隔的时间(因此频率)值?
先谢了!

vojdkbi0

vojdkbi01#

错误与日志或其他数据间隔无关。它只是锡上写的:你正在尝试计算插值器在插值范围(x.min(),x.max())之外。
为了避免错误,将bounds_error=False传递给interp1d。然而,外推的质量可能会有很大的变化,所以YMMV。我建议重新评估这种方法,并弄清楚 * 为什么 * 你要外推,以及是否可以使用一些外部知识来代替。

相关问题