重采样，插值矩阵答案

【问题标题】：resampling, interpolating matrix重采样，插值矩阵
【发布时间】：2023-03-06 14:25:02
【问题描述】：

我正在尝试插入一些数据以进行绘图。例如，给定 N 个数据点，我希望能够生成一个“平滑”图，由 10*N 左右的插值数据点组成。

我的方法是生成一个 N×10*N 矩阵并计算原始向量和我生成的矩阵的内积，得到一个 1×10*N 向量。我已经计算出我想用于插值的数学，但我的代码很慢。我对 Python 很陌生，所以我希望这里的一些专家可以给我一些想法，让我可以尝试加快我的代码速度。

我认为部分问题在于生成矩阵需要对以下函数进行 10*N^2 次调用：

def sinc(x):
    import math
    try:
        return math.sin(math.pi * x) / (math.pi * x)
    except ZeroDivisionError:
        return 1.0

（这个comes from sampling theory。本质上，我正在尝试从其样本中重新创建一个信号，并将其上采样到更高的频率。）

矩阵由以下生成：

def resampleMatrix(Tso, Tsf, o, f):
    from numpy import array as npar
    retval = []

    for i in range(f):
        retval.append([sinc((Tsf*i - Tso*j)/Tso) for j in range(o)])

    return npar(retval)

我正在考虑将任务分解成更小的部分，因为我不喜欢将 N^2 矩阵放在内存中的想法。我可能可以将“resampleMatrix”变成一个生成器函数并逐行执行内积，但我认为这不会大大加快我的代码速度，直到我开始将内容分页进出内存。

提前感谢您的建议！

【问题讨论】：

除了您尝试对代码执行的操作之外，您可以在没有数据生成模型的情况下仅插入额外点的想法是错误的。如果您想以任何统计原则的方式执行此操作，则需要执行某种回归。见en.wikipedia.org/wiki/Generative_model
看起来 Phil 只想使用插值进行绘图。只要插值点不用于其他目的，我不明白为什么需要生成模型
@Phil：考虑到它是 O(N^2) 算法，而其他方法（如三次样条）只有 O(N)，您想使用 sinc 插值的任何特殊原因？
@twole18：数据的模型是按照en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem采样的。您可以使用 sinc 函数完全恢复原始数据。
numpy 已经有一个 sinc() 函数，顺便说一下。 docs.scipy.org/doc/numpy/reference/generated/numpy.sinc.html

标签： python matrix interpolation signal-processing generator

【解决方案1】：

小改进。使用在编译后的 C 代码中运行的内置 numpy.sinc(x) 函数。

可能的更大改进：您可以动态进行插值（当绘图发生时）？还是您绑定到只接受矩阵的绘图库？

【讨论】：

感谢您的评论。奇怪的是，当我使用 numpy.sinc(x) 时，代码运行速度慢了大约 10 倍。我很惊讶！
描述的绘图部分仅用于说明目的。我并不真的担心绘制情节，只是让实际计算更快。最终这将更像是一个“即时”类型的任务，因为我将处理大型数据集的切片。然而，就目前而言，运行我认为最小的有用数据片需要比下一个数据集到达所需的时间更多......
Tso = 初始采样时间，Tsf = 最终采样时间。因此，如果我从一个以 1kHz 采样的信号开始，并且我想为每个采样生成 10 个插值点（新的采样率为 10kHz），则 Tso = 0.001，Tsf = 0.0001。

【解决方案2】：

您的问题并不完全清楚；您正在尝试优化您发布的代码，对吗？

像这样重写 sinc 应该会大大加快速度。这种实现避免了在每次调用时检查数学模块是否被导入，不会进行三次属性访问，并且将异常处理替换为条件表达式：

from math import sin, pi
def sinc(x):
    return (sin(pi * x) / (pi * x)) if x != 0 else 1.0

您还可以尝试通过直接创建 numpy.array（不是从列表列表）来避免创建两次矩阵（并在内存中并行保存两次）：

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.zeros((f, o))
    for i in xrange(f):
        for j in xrange(o):
            retval[i][j] = sinc((Tsf*i - Tso*j)/Tso)
    return retval

（在 Python 3.0 及更高版本上将 xrange 替换为 range）

最后，您可以使用 numpy.arange 创建行，也可以在每一行甚至整个矩阵上调用 numpy.sinc：

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.zeros((f, o))
    for i in xrange(f):
        retval[i] = numpy.arange(Tsf*i / Tso, Tsf*i / Tso - o, -1.0)
    return numpy.sinc(retval)

这应该比您的原始实现快得多。尝试这些想法的不同组合并测试它们的性能，看看哪个效果最好！

【讨论】：

“用条件表达式替换异常处理”但异常比python中的条件更快。同样，使用一次pi*x 并使用两次会更快，对吧？
@endolith “异常比 Python 中的条件更快”是不正确的，这实际上取决于异常情况发生的频率。无论如何，与避免在每个函数调用上进行导入和属性查找相比，这应该是微不足道的。此处不使用 try/except 是样式和代码清晰度的问题。
@endolith 至于pi * x，我不确定创建一个新的局部变量以避免单个浮点乘法是否有益。这是您只需要测试的事情之一。不过，与我建议的其他更改相比，这确实微不足道，但会产生很大的影响。
是的，异常比条件更快，所以如果它们很少发生，使用它们的代码也会更快。在这种情况下，条件只会在输入恰好为 0 时发生，这种情况非常罕见，因此使用异常会更快。在快速测试中，异常版本的随机输入速度提高了约 30%，而使用pix = pi*x 的速度也提高了约 40%。

【解决方案3】：

如果您想以一种非常通用且快速的方式对数据进行插值，则样条曲线或多项式非常有用。 Scipy 有 scipy.interpolate 模块，非常有用。您可以在官方页面找到many examples。

【讨论】：

【解决方案4】：

这是一个使用 scipy 进行一维插值的最小示例 - 没有重新发明那么有趣，但是。
情节看起来像sinc，这并非巧合：尝试谷歌样条重新采样“近似正弦”。
（大概更少的局部/更多的抽头⇒更好的近似，但我不知道本地 UnivariateSplines 是怎样的。）

""" interpolate with scipy.interpolate.UnivariateSpline """
from __future__ import division
import numpy as np
from scipy.interpolate import UnivariateSpline
import pylab as pl

N = 10 
H = 8
x = np.arange(N+1)
xup = np.arange( 0, N, 1/H )
y = np.zeros(N+1);  y[N//2] = 100

interpolator = UnivariateSpline( x, y, k=3, s=0 )  # s=0 interpolates
yup = interpolator( xup )
np.set_printoptions( 1, threshold=100, suppress=True )  # .1f
print "yup:", yup

pl.plot( x, y, "green",  xup, yup, "blue" )
pl.show()

2010 年 2 月添加：另请参阅 basic-spline-interpolation-in-a-few-lines-of-numpy

【讨论】：

【解决方案5】：

我不太确定您要做什么，但是您可以通过一些加速来创建矩阵。 Braincore's suggestion 使用 numpy.sinc 是第一步，但第二步是要意识到 numpy 函数希望在 numpy 数组上工作，它们可以在 C speen 中执行循环，并且可以比单个元素更快地完成。

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.sinc((Tsi*numpy.arange(i)[:,numpy.newaxis]
                         -Tso*numpy.arange(j)[numpy.newaxis,:])/Tso)
    return retval

诀窍在于，通过使用 numpy.newaxis 对 arange 进行索引，numpy 将形状为 i 的数组转换为形状为 i x 1 的数组，并将形状为 j 的数组转换为形状为 1 x j 的数组。在减法步骤中，numpy 将“广播”每个输入以充当 i x j 形状的数组并进行减法。（“广播”是 numpy 的术语，反映了没有额外复制将 i x 1 拉伸到 i x j 的事实。）

现在 numpy.sinc 可以遍历编译代码中的所有元素，比您编写的任何 for 循环都快得多。

（如果您在减法之前进行除法，则可以提供额外的加速，特别是因为在减法中，除法取消了乘法。）

唯一的缺点是您现在需要支付额外的 Nx10*N 数组来保存差异。如果 N 很大并且内存是一个问题，这可能会破坏交易。

否则，您应该可以使用numpy.convolve 编写此内容。从我刚刚学到的关于 sinc-interpolation 的知识来看，我想说你想要像 numpy.convolve(orig,numpy.sinc(numpy.arange(j)),mode="same") 这样的东西。但我可能对细节有误。

【讨论】：

我正在尝试卷积，所以我认为 numpy.convolve 可能是正确的方向。

【解决方案6】：

如果您唯一的兴趣是“生成“平滑”图”，我会选择一个简单的多项式样条曲线拟合：

对于任何两个相邻的数据点，三次多项式函数的系数可以根据这些数据点的坐标以及它们左右两侧的两个附加点（不考虑边界点）来计算。这将生成一个很好的点具有连续一阶导数的平滑曲线。有一个直接的公式可以将 4 个坐标转换为 4 个多项式系数，但我不想剥夺您查找它的乐趣；o)。

【讨论】：

【解决方案7】：

我建议您检查您的算法，因为它是一个重要的问题。具体来说，我建议您访问由 Hu 和 Pavlidis（1991 年）撰写的文章“使用圆锥样条进行函数绘图”（IEEE 计算机图形学和应用程序）。他们的算法实现允许对函数进行自适应采样，因此渲染时间比规则间隔的方法更短。

摘要如下：

提出了一种方法，给定一个 a的数学描述函数，一个圆锥样条逼近生成函数图。圆锥弧被选为原始曲线，因为有简单的增量绘图算法对于一些已经包含的圆锥曲线设备驱动程序，并且有简单的局部近似算法圆锥曲线。拆分合并算法用于自适应地选择结，根据形状分析基于其原始功能一阶导数，是介绍了。

【讨论】：

我的算法来源于抽样理论。本质上，我正在尝试从其样本中重新创建一个信号，并以更高的频率对其进行重新采样。出于绘图的目的，我确信我的解决方案不是最好的方法......
@Phil：你应该在问题中这么说

【解决方案8】：

这是上采样。有关示例解决方案，请参阅 Help with resampling/upsampling。

执行此操作的一种快速方法（对于离线数据，例如您的绘图应用程序）是使用 FFT。这就是 SciPy 的原生 resample() function 所做的。不过，它假设一个周期性信号so it's not exactly the same。见this reference：

这是关于时域实信号插值的第二个问题，这确实是一件大事。只有当原始 x(n) 序列在其整个时间间隔内是周期性的时，这种精确的插值算法才能提供正确的结果。

您的函数假定信号的样本在定义范围之外都是 0，因此这两种方法将偏离中心点。如果你先用很多零填充信号，它将产生非常接近的结果。图中未显示的边缘还有几个零：

三次插值不适用于重采样。此示例是一个极端情况（接近采样频率），但正如您所见，三次插值甚至不接近。对于较低的频率，它应该非常准确。

【讨论】：

感谢您的回答！ @endolith 我注意到您在下面的评论。你说得对，我应该从一开始就让我的问题更清楚。