9

プロットする目的でいくつかのデータを補間しようとしています。たとえば、N個のデータポイントが与えられた場合、10*N程度の補間されたデータポイントで構成される「スムーズな」プロットを生成できるようにしたいと思います。

私のアプローチは、N x 10 * Nの行列を生成し、元のベクトルと生成した行列の内積を計算して、1 x 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()関数が行うことです。ただし、周期的な信号を想定しているため、まったく同じではありません。このリファレンスを参照してください:

これは、時間領域の実信号補間に関する 2 番目の問題であり、これは非常に重要です。この正確な補間アルゴリズムは、元の x(n) シーケンスが完全な時間間隔内で周期的である場合にのみ、正しい結果を提供します。

関数は、信号のサンプルが定義された範囲外ではすべて 0 であると想定しているため、2 つの方法は中心点から分岐します。最初に多数のゼロで信号をパディングすると、非常に近い結果が得られます。ここには示されていませんが、プロットの端を越えてさらにいくつかのゼロがあります。

ここに画像の説明を入力

リサンプリングの目的では、3 次補間は正しくありません。この例は極端なケース (サンプリング周波数に近い) ですが、ご覧のとおり、キュービック補間はそれに近いものではありません。より低い周波数では、かなり正確なはずです。

于 2011-09-21T01:28:08.653 に答える
3

非常に一般的で高速な方法でデータを補間したい場合は、スプラインまたは多項式が非常に便利です。Scipy には、非常に便利な scipy.interpolate モジュールがあります。公式ページで多くの例を見つけることができます。

于 2009-12-05T11:51:51.200 に答える
1

'「滑らかな」プロットを生成することだけが目的の場合は、単純な多項式スプラインカーブフィットを使用します。

任意の2つの隣接するデータポイントについて、3次多項式関数の係数は、それらのデータポイントの座標と、左右の2つの追加ポイント（境界ポイントを無視）から計算できます。これにより、次のような滑らかな曲線上にポイントが生成されます。継続的な一次微分。4つの座標を4つの多項式係数に変換するための簡単な式がありますが、それを調べる楽しみを奪いたくありません; o）。

于 2009-12-07T06:28:59.987 に答える
1

これは、scipy を使用した 1 次元補間の最小限の例です。再発明ほど楽しいものではありませんが。
プロットは のようになりますがsinc、これは偶然ではありません。Google スプラインのリサンプル「近似 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は最初のステップですが、2 つ目は、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 の配列を形状 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

あなたの質問は完全に明確ではありません。投稿したコードを最適化しようとしていますよね？

このように sinc を書き直すと、かなり高速になるはずです。この実装では、すべての呼び出しで math モジュールがインポートされているかどうかのチェックが回避され、属性へのアクセスが 3 回行われず、例外処理が条件式に置き換えられます。

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

また、(リストのリストからではなく) numpy.array を直接作成することで、行列を 2 回作成すること (およびメモリ内で並列に 2 回保持すること) を回避することもできます。

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 に答える
0

重大な問題であるため、アルゴリズムを確認することをお勧めします。具体的には、Hu and Pavlidis（1991）の記事「円錐スプラインを使用した関数プロット」（IEEEコンピューターグラフィックスとアプリケーション）にアクセスすることをお勧めします。それらのアルゴリズムの実装により、関数の適応サンプリングが可能になり、レンダリング時間は規則的な間隔のアプローチよりも短くなります。

要約は次のとおりです。

関数の数学的記述が与えられると、関数のプロットを近似する円錐スプラインが生成される方法が提示されます。一部のデバイスドライバーにすでに含まれている円錐曲線の単純な増分プロットアルゴリズムがあり、円錐曲線による局所近似の単純なアルゴリズムがあるため、円錐曲線がプリミティブ曲線として選択されました。一次導関数に基づく元の関数の形状分析に従って、ノットを適応的に選択するためのスプリットアンドマージアルゴリズムが導入されています。

于 2009-12-07T06:38:17.663 に答える
0

小さな改善。コンパイルされた C コードで実行される組み込みの numpy.sinc(x) 関数を使用します。

より大きな改善の可能性：オンザフライで補間できますか（プロットが発生すると）？それとも、マトリックスのみを受け入れるプロット ライブラリに縛られていますか?

于 2009-12-05T06:55:06.047 に答える