如何使用 scipy.interpolate.interp1d 插值对数均匀间隔的数据?

问题描述 投票:0回答:1

我有一个实值函数 G(t),我想在对数均匀间隔时间上进行插值:formula1。这是为了获得 G 的分量作为频率 formula2 而不是时间的函数的近似表达式。因此,我这样定义了一个时间列表:

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

提前致谢!

python scipy interpolation linear-interpolation
1个回答
0
投票

为避免错误,将 bounds_error=False 传递给 interp1d。然而,外推的质量可能千差万别,所以 YMMV。我建议重新评估该方法并弄清楚

为什么

你要推断,以及你是否可以使用一些外部知识。

© www.soinside.com 2019 - 2024. All rights reserved.