重采样，插值矩阵

Question

我正试图插入一些数据以便绘图。例如，给定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矩阵。我可能会把‘重采样矩阵’变成一个生成器函数，然后一行一行地执行内部产品，但我认为在开始对内存中的内容进行分页之前，这不会大大加快代码的速度。

提前谢谢你的建议！

Answer 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

（在Python3.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)

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

Answer 2

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

一种快速的方法是使用fft（对于脱机数据，如绘图应用程序）。这就是SciPy的本地^{} function所做的。不过，它假设有一个周期性信号，so it's not exactly the same。见this reference：

Here’s the second issue regarding time-domain real signal interpolation, and it’s a big deal indeed. This exact interpolation algorithm provides correct results only if the original x(n) sequence is periodic within its full time inter­val.

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

三次插值对于重采样来说是不正确的。这个例子是一个极端的例子（接近采样频率），但是正如你所看到的，三次插值甚至都不接近。对于较低的频率，它应该相当准确。

Answer 3

如果您想以一种非常通用和快速的方式插值数据，样条曲线或多项式是非常有用的。Scipy有Scipy.interpolate模块，非常有用。你可以在官方网页上找到many examples。