16

私が行っているいくつかの計算作業をプロファイリングすると、私のプログラムの1つのボトルネックは、基本的にこれを実行する関数であることがわかりました(npis numpyspis scipy):

def mix1(signal1, signal2):
    spec1 = np.fft.fft(signal1, axis=1)
    spec2 = np.fft.fft(signal2, axis=1)
    return np.fft.ifft(spec1*spec2, axis=1)

両方の信号の形状(C, N)C、データセットの数(通常は20未満)とN各セットのサンプル数(約5000)です。各セット(行)の計算は、他のセットから完全に独立しています。

これは単なる畳み込みであると考えたので、次のように置き換えようとしました。

def mix2(signal1, signal2):
    outputs = np.empty_like(signal1)

    for idx, row in enumerate(outputs):
        outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same')

    return outputs

...同じ結果が得られたかどうかを確認するためだけに。しかし、私はしませんでした、そして私の質問は次のとおりです:

  1. なぜだめですか?
  2. に相当するものを計算するためのより良い方法はありmix1()ますか?

mix2おそらく、そのままでは高速ではなかったと思いますが、並列化の開始点としては適切だったかもしれません。)

これをすばやく確認するために使用した完全なスクリプトは次のとおりです。

import numpy as np
import scipy as sp
import scipy.signal

N = 4680
C = 6

def mix1(signal1, signal2):
    spec1 = np.fft.fft(signal1, axis=1)
    spec2 = np.fft.fft(signal2, axis=1)
    return np.fft.ifft(spec1*spec2, axis=1)

def mix2(signal1, signal2):
    outputs = np.empty_like(signal1)

    for idx, row in enumerate(outputs):
        outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same')

    return outputs

def test(num, chans):
    sig1 = np.random.randn(chans, num)
    sig2 = np.random.randn(chans, num)
    res1 = mix1(sig1, sig2)
    res2 = mix2(sig1, sig2)

    np.testing.assert_almost_equal(res1, res2)

if __name__ == "__main__":
    np.random.seed(0x1234ABCD)
    test(N, C)
4

3 に答える 3

11

だから私はこれをテストし、今いくつかのことを確認することができます:

1)numpy.convolveは循環的ではありません。これは、fftコードが提供するものです。

2)FFTは内部で2の累乗にパディングしません。次の操作の大きく異なる速度を比較します。

x1 = np.random.uniform(size=2**17-1)
x2 = np.random.uniform(size=2**17)

np.fft.fft(x1)
np.fft.fft(x2)

3)正規化は違いではありません。a(k)* b(ik)を合計して単純な巡回畳み込みを行うと、FFTコードの結果が得られます。

2の累乗にパディングすると、答えが変わります。長さの素因数(数値レシピで言及されているがコード化されていない)を巧みに使用することでこれに対処する方法があるという話を聞いたことがありますが、実際にそうする人を見たことがありません。

于 2011-07-28T07:27:04.520 に答える
2

scipy.signal.fftconvolveはFFTによって畳み込みます。これはPythonコードです。ソースコードを調べて、mix1関数を修正することができます。

于 2011-07-29T00:40:31.547 に答える
1

前述のように、scipy.signal.convolve関数は巡回畳み込みを実行しません。(fftを使用するのとは対照的に)実空間で巡回畳み込みを実行する場合は、scipy.ndimage.convolve関数を使用することをお勧めします。'wrap'に設定できるモードパラメータがあり、巡回畳み込みになります。

for idx, row in enumerate(outputs):
    outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode='wrap')
于 2013-08-20T15:21:25.750 に答える