为什么 Matlab interp1 产生的结果与 numpy interp 不同？答案

【问题标题】：Why does Matlab interp1 produce different results than numpy interp?为什么 Matlab interp1 产生的结果与 numpy interp 不同？
【发布时间】：2019-09-25 17:18:53
【问题描述】：

编辑： 编辑代码以产生与 Matlab 一致的结果。见下文。

我正在将 Matlab 脚本转换为 Python，在某些情况下线性插值结果会有所不同。我想知道为什么以及是否有任何方法可以解决这个问题？

这是 Matlab 和 Python 中的代码示例以及结果输出（请注意，在这种情况下 t 恰好等于 tin）：

MATLAB：

t= [ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667, 736696.0125];
tin =[ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667,  736696.0125];
xin = [   nan , 1392.,  1406. , 1418.  ,  nan , 1442. ,   nan];

interp1(tin,xin,t)

ans =

 NaN        1392        1406        1418         NaN        1442         NaN

Python (numpy)：

（scipy interpolate.interp1d 产生与 numpy 相同的结果）

t= [ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667, 736696.0125];
tin =[ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667,  736696.0125];
xin = [   nan , 1392.,  1406. , 1418.  ,  nan , 1442. ,   nan];

x = np.interp(t,tin,xin)

array([   nan,  1392.,  1406.,    nan,    nan,    nan,    nan]) 

# Edit
# Find indices where t == tin and if the np.interp output 
# does not match the xin array, overwrite the np.interp output at those 
# indices 
same = np.where(t == tin)[0]
not_same = np.where(xin[same] != x[same])[0]
x[not_same] = xin[not_same]

【问题讨论】：

我不知道 Python 但 Matlab 有 9 种不同的方法可供interp1 使用... Python 也必须有几种方法。由于您没有在上面的任何代码中明确指定方法，您确定两种语言的 default 设置相同吗？
顺便说一句，您的问题与that one 非常相似，在那里您已经有 cmets 要求您显示一些代码。如果您是同一个人，最好通过添加此处显示的信息和代码来编辑另一个问题（当另一个问题完成后，您可以删除这个以避免多余的问题）。
这是因为你的 y 值有 nans。据推测，numpy 实现做了类似y1 + (y2-y1)*(x-x1)/(x2-x1) 的事情，nan 被传播到结果，matlab 代码有一个特殊的检查if (x == x1) return y1。
@Hoki 我知道这个问题，但不是问这个问题的人。我通过添加我的代码来回复它（因为其中一位评论者要求它）但版主删除了它，所以我决定制作一个新的。另外，numpy 和 matlab 默认都是线性插值。
@user545424 你是绝对正确的。我现在正在检查 x == x1 并返回 y1 的情况，我得到的结果与 matlab 一致。

标签： python matlab numpy interpolation

【解决方案1】：

看起来，Matlab 在它的插值中包含了一个额外的相等性检查。

线性一维插值通常通过找到两个跨越输入值x的x值来完成，然后将结果计算为：

y = y1 + (y2-y1)*(x-x1)/(x2-x1)

如果您传入的x 值恰好等于输入x 坐标之一，则例程通常会计算正确的值，因为x-x1 将为零。但是，如果您的输入数组有一个 nan 作为 y1 或 y2 这些将传播到结果。

根据您发布的代码，我最好的猜测是 Matlab 的插值函数有一个额外的检查，类似于：

if x == x1:
    return y1

并且numpy函数没有这个检查。

要在 numpy 中实现相同的效果，您可以这样做：

np.where(t == tin,xin,np.interp(t,tin,xin))

【讨论】：