-
13-09-2019 - |
题
我正在尝试插入一些数据以进行绘图。例如,给定 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”制作成生成器函数并逐行进行内积,但我认为这不会大大加快我的代码速度,直到我开始将内容分页进出内存。
预先感谢您的建议!
解决方案
这是上采样。请参阅与一些示例性解决方案的重采样/上采样帮助。
有一个快速的方法来做到这一点(用于离线数据,如绘图应用程序)是使用的FFT。这是SciPy的原生 resample()
功能做什么。它假设一个周期信号,但是,,所以它不是完全一样。请参阅此参考一个>:
下面是关于时域实信号的插值第二个问题,这是一个大问题,确实如此。这个精确的内插算法提供正确的结果只有当原始x(n)的序列是其完整的时间间隔内周期性的。
您函数假定所述信号的样本都是0所定义的范围之外,所以这两种方法将发散从中心点路程。如果你垫零分的信号首先,它会产生一个非常接近的结果。有几种多个零过去情节的边缘未在此示出:
三次插值不会对重新采样的目的是正确的。这个例子是一个极端的情况下(靠近采样频率),但你可以看到,立方插值甚至还没有接近。对于较低频率应该是相当准确的。
其他提示
如果你想在一个相当普遍,快捷的方式来插值数据,样条或多项式是非常有用的。 SciPy的有scipy.interpolate模块,这是非常有用的。你可以找到在官方页面中的范例。
您的问题是不完全清楚;你试图优化您发布的代码吧?
像这样重新写入正弦应极大地加快速度。这种实现避免了检查数学模块在每次调用进口,不会执行属性访问三次,并取代异常与条件表达式处理:
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及以上)
最后,可以创建具有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)
这应该比你原来的执行显著更快。试试这些想法不同的组合,并测试其性能,看看哪出最好的!
我不太清楚你想要做什么,但还是有一些的加速可以做,以创建矩阵。 Braincore的使用numpy.sinc
建议是第一步,但第二个是实现numpy的功能希望在numpy的阵列,在那里他们可以位于C speen do循环工作,并且可以比单独的元件做得更快。
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的aranges,numpy的具有形状i中的阵列转换成一个具有形状I X 1,和具有形状j中的阵列,塑造1×Ĵ。在减法步骤,将numpy的“广播”的每个输入到充当I XĴ成形阵列和做减法。 (“广播”是numpy的的术语,反映没有额外的拷贝被制成拉伸I X 1到i X j中的事实。)
现在的numpy.sinc可以通过在编译代码的所有元素,比任何for循环您可以编写更快的迭代。
(如果减法之前做除法有一个额外的加速可用,特别是因为在矿井后者分割取消乘法。)
唯一的缺点是,你现在支付额外NX10 * N阵列持有的差异。这可能是一个因素在于如果N较大,存储器是一个问题。
否则,你应该能够写这个使用numpy.convolve
。从一点点我刚刚了解正弦插值,我说你想要的东西,像numpy.convolve(orig,numpy.sinc(numpy.arange(j)),mode="same")
。但我可能是错的有关细节。
如果你唯一的兴趣是“生成‘平滑’阴谋”我只想去一个简单的多项式样条曲线拟合:
有关的任何两条相邻的数据点的三次多项式函数的系数可以从这些数据点的坐标和所述两个额外的点到其左侧和右侧被计算(不考虑边界点)。这将在一个很好的产生点与连续的第一dirivitive平滑曲线。有转换4个坐标,以4个多项式系数直线前进公式,但我不想剥夺你的寻找它的乐趣; O)。
这是使用 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 几行中的基本样条插值
小的改善。使用运行在编译的C代码内置numpy.sinc(x)函数。
可能的改进较大:你能做到在飞行插值(如绘图时)?或者是你绑只接受一个矩阵绘图库?
我建议你检查你的算法,因为它是一个不平凡的问题。具体而言,我建议你访问的文章“功能绘制使用圆锥样条”胡锦涛和Pavlidis(1991)(IEEE计算机图形和应用程序)。他们的算法实现允许功能的自适应采样,使得绘制时间大于与规则间隔的方法要小。
抽象如下:
一个方法被提出,其中,给定一个 的数学描述 功能,圆锥样条逼近 函数的曲线图被产生。 圆锥弧被选为 原始曲线,因为有 简单的增量绘制算法 为二次曲线已经包括在一些 设备驱动程序,并有简单的 算法局部逼近由 圆锥曲线。一个分裂 - 合并算法 用于自适应选择疙瘩, 根据形状的分析 基于原有功能的 一阶导数,是 引入的。