logo标签列表

Resampling，插值矩阵

pythonmatrixinterpolationsignal-processinggenerator
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'变成生成器函数，并逐行执行内积，但我认为在开始对内存进行分页之前，这不会显著加快我的代码速度。
提前感谢您的建议！
2
完全不管你在代码中想做什么，仅仅通过插值额外的点而没有数据生成模型的想法是错误的。如果你想以任何统计原则的方式进行这样的操作，你需要执行某种回归分析。请参见http://en.wikipedia.org/wiki/Generative_model。 - twolfe18
看起来Phil只想使用插值来绘图。只要这些插值点没有被用于其他目的，我就不明白为什么还需要一个生成模型。 - Jitse Niesen
1
@Phil：你想使用sinc插值法有特别的原因吗？毕竟它是O(N^2)算法，而其他方法如三次样条只有O(N)。 - Jitse Niesen
2
@twole18：数据的模型是根据http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem进行采样的。您可以使用sinc函数*完全*恢复原始数据。 - endolith
顺便说一下，numpy已经有了sinc()函数。http://docs.scipy.org/doc/numpy/reference/generated/numpy.sinc.html - endolith
8个回答
9
这是上采样。有一些示例解决方案，请参见帮助重新采样/上采样。
一个快速的方法（用于离线数据，如您的绘图应用程序）是使用FFT。这就是SciPy的本机resample()函数所做的。它假定一个周期性信号，所以它不完全相同。请参见此参考文献：

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

您的函数假定信号的样本在定义范围之外全部为0，因此这两种方法将远离中心点。如果您首先用许多零填充信号，则会产生非常接近的结果。在这里未显示的绘图边缘有几个更多的零：
立方插值不适用于重新采样。这个例子是一个极端情况（接近采样频率），但是正如您所看到的，立方插值甚至都不接近。对于较低的频率，它应该相当准确。
1
谢谢回答！@endolith 我注意到你下面的评论。你说得对，我应该从一开始就让我的问题更清晰明了。 - Phil
3
如果您想以一种通用且快速的方式插值数据，样条或多项式非常有用。Scipy拥有scipy.interpolate模块，非常实用。您可以在官方页面中找到许多示例。
1

你的问题不是很清楚，你是想优化你发布的代码，对吗？

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

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
@endolith，“在Python中，异常比条件语句更快”这种说法并不正确，它实际上取决于异常情况发生的频率。无论如何，与避免在每个函数调用时进行导入和属性查找相比，这应该是相当微不足道的。在这里不使用try/except只是一种风格和代码清晰度的问题。 - taleinat
@endolith 关于 pi * x，我不确定为了避免单个浮点数乘法而创建一个新的局部变量是否有益。这是你必须测试的事情之一。然而，与我建议的其他更改相比，这真的微不足道，那些更改将产生巨大的影响。 - taleinat
是的，异常比条件语句更快，因此如果它们很少发生，使用它们的代码也会更快。在这种情况下，条件语句只会在输入恰好为0时发生，这是非常罕见的，因此使用异常更快。在一个快速测试中，对于随机输入，异常版本大约快30％，而使用pix = pi*x也可以加快40％左右的速度。 - endolith
1
这是一个使用Scipy进行一维插值的最小示例 - 不如重新发明有趣，但是。 图表看起来像sinc函数，这不是巧合： 尝试谷歌样条重采样“近似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月新增：另请参阅使用少量NumPy代码进行基本样条插值

1

我不太确定你想做什么，但是有一些加速方法可以用来创建矩阵。Braincore的建议使用numpy.sinc是第一步，但第二步是意识到numpy函数希望在numpy数组上工作，在那里它们可以以C速度循环，并且比在单个元素上更快。

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对数组进行索引时，numpy会将形状为i的数组转换为形状为i x 1的数组，将形状为j的数组转换为形状为1 x j的数组。在减法步骤中，numpy将“broadcast”每个输入以作为一个i x j形状的数组并执行减法。（“Broadcast”是numpy的术语，反映了没有额外的复制来将i x 1拉伸到i x j。）
现在，numpy.sinc可以在编译代码中迭代所有元素，比您写的任何for循环都要快得多。
（如果您在减法之前进行除法，尤其是在后者中，除法可以取消乘法，那么还可以获得额外的加速。）
唯一的缺点是现在需要支付额外的Nx10*N数组来保存差异。如果N很大且内存是问题，则可能会造成交易破裂。
否则，您应该能够使用numpy.convolve编写此内容。从我刚学到的有关sinc插值的少量知识中，我认为您需要像这样的东西：numpy.convolve（orig，numpy.sinc（numpy.arange（j）），mode =“same”）。但是我可能对具体情况错了。
我正在尝试进行卷积计算，因此我认为使用numpy.convolve可能是正确的方向。 - Phil
1

如果您只是想要“生成一个平滑的图形”，我建议使用简单的多项式样条曲线拟合：

对于任意相邻的两个数据点，可以从这些数据点的坐标和它们左右两侧的另外两个点（忽略边界点）计算出三次多项式函数的系数。这将生成一条漂亮平滑的曲线，并具有连续的一阶导数。有一个简单的公式可以将4个坐标转换为4个多项式系数，但我不想剥夺您查找它的乐趣; o)。

0

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

可能的更大改进：您能否在绘图发生时动态进行插值？或者您是否被绘图库约束只能接受矩阵？

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

我建议您检查算法，因为这是一个非平凡的问题。具体而言，我建议您获取胡和帕夫利迪斯（1991年）在IEEE计算机图形和应用中发表的文章“使用圆锥样条绘制函数”的访问权限。他们的算法实现允许自适应采样函数，使得渲染时间比常规间隔方法更短。

摘要如下：

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

1
我的算法来自采样理论。本质上，我试图从样本中重新创建信号，并以更高的频率重新采样它。出于绘图的目的，我确信我的解决方案不是最佳方法... - Phil
@Phil：你应该在问题中说清楚。 - endolith
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接