9

我正在尝试插入一些数据以进行绘图。例如，给定 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

（这来自采样理论。本质上，我试图从其样本中重新创建一个信号，并将其上采样到更高的频率。）

矩阵由以下生成：

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”变成一个生成器函数并逐行执行内积，但我认为这不会大大加快我的代码速度，直到我开始将内容分页进出内存。

提前感谢您的建议！

4

8 回答 8

9

这是上采样。有关一些示例解决方案，请参阅重采样/上采样帮助。

一种快速的方法（对于离线数据，如绘图应用程序）是使用 FFT。这就是 SciPy 的原生resample()函数所做的。但是，它假定一个周期性信号，因此并不完全相同。请参阅此参考：

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

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

在此处输入图像描述

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

于 2011-09-21T01:28:08.653 回答
3

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

于 2009-12-05T11:51:51.200 回答
1

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

像这样重写 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)

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

于 2009-12-05T11:10:08.290 回答
1

这是一个使用 scipy 进行一维插值的最小示例——没有重新发明那么有趣，但是。
情节看起来像sinc，这并非巧合：尝试 google spline resample "approximate 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

于 2009-12-05T17:10:34.863 回答
1

我不太确定你想做什么，但你可以做一些加速来创建矩阵。 Braincore 的使用建议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 的数组转换为形状为 ix 1 的数组，并将形状为 j 的数组转换为形状为 1 x j 的数组。在减法步骤中，numpy 将“广播”每个输入以充当 aixj 形状的数组并进行减法。（“广播”是 numpy 的术语，反映了没有制作额外的副本来将 ix 1 延伸到 ix j 的事实。）

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

（如果你在减法之前进行除法，还有一个额外的加速，特别是因为在减法中，除法取消了乘法。）

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

否则，您应该能够使用numpy.convolve. 从我刚刚学到的关于 sinc-interpolation 的一些知识中，我会说你想要类似numpy.convolve(orig,numpy.sinc(numpy.arange(j)),mode="same"). 但我可能对细节有误。

于 2009-12-06T05:57:16.910 回答
1

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

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

于 2009-12-07T06:28:59.987 回答
0

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

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

于 2009-12-05T06:55:06.047 回答
0

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

摘要如下：

提出了一种方法，通过该方法，给定函数的数学描述，生成逼近函数图的圆锥样条曲线。选择圆锥弧作为原始曲线是因为一些设备驱动程序中已经包含简单的圆锥曲线增量绘图算法，并且有简单的圆锥曲线局部逼近算法。介绍了一种基于原函数一阶导数的形状分析自适应选择节点的分裂合并算法。

于 2009-12-07T06:38:17.663 回答