私の畳み込みルーチンはnumpyとscipyの違いとは何ですか？

Question

私は時系列分類のためにカーネルで遊んでいたので、手動で1次元畳み込みをコーディングしたいと思っていました。ここで見られるように、有名なウィキペディア畳み込みイメージを作ることにしました。

ここに私のスクリプトがあります。私はstandard formula for convolution for a digital signalを使用しています。

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.ndimage 

plt.style.use('ggplot') 

def convolve1d(signal, ir): 
    """ 
    we use the 'same'/'constant' method for zero padding. 
    """ 
    n = len(signal) 
    m = len(ir) 
    output = np.zeros(n) 

    for i in range(n): 
     for j in range(m): 
      if i - j < 0: continue 
      output[i] += signal[i - j] * ir[j] 

    return output 

def make_square_and_saw_waves(height, start, end, n): 
    single_square_wave = [] 
    single_saw_wave = [] 
    for i in range(n): 
     if start <= i < end: 
      single_square_wave.append(height) 
      single_saw_wave.append(height * (end-i)/(end-start)) 
     else: 
      single_square_wave.append(0) 
      single_saw_wave.append(0) 

    return single_square_wave, single_saw_wave 

# create signal and IR 
start = 40 
end = 60 
single_square_wave, single_saw_wave = make_square_and_saw_waves(
    height=10, start=start, end=end, n=100) 

# convolve, compare different methods 
np_conv = np.convolve(
    single_square_wave, single_saw_wave, mode='same') 

convolution1d = convolve1d(
    single_square_wave, single_saw_wave) 

sconv = scipy.ndimage.convolve1d(
    single_square_wave, single_saw_wave, mode='constant') 

# plot them, scaling by the height 
plt.clf() 
fig, axs = plt.subplots(5, 1, figsize=(12, 6), sharey=True, sharex=True) 

axs[0].plot(single_square_wave/np.max(single_square_wave), c='r') 
axs[0].set_title('Single Square') 
axs[0].set_ylim(-.1, 1.1) 

axs[1].plot(single_saw_wave/np.max(single_saw_wave), c='b') 
axs[1].set_title('Single Saw') 
axs[2].set_ylim(-.1, 1.1) 

axs[2].plot(convolution1d/np.max(convolution1d), c='g') 
axs[2].set_title('Our Convolution') 
axs[2].set_ylim(-.1, 1.1) 

axs[3].plot(np_conv/np.max(np_conv), c='g') 
axs[3].set_title('Numpy Convolution') 
axs[3].set_ylim(-.1, 1.1) 

axs[4].plot(sconv/np.max(sconv), c='purple') 
axs[4].set_title('Scipy Convolution') 
axs[4].set_ylim(-.1, 1.1) 

plt.show() 

そして、ここでは、私が手にプロットです：

あなたが見ることができるように、何らかの理由で、私の畳み込みがシフトしています。曲線の数値（y値）は同じですが、フィルタ自体のサイズの約半分にシフトしています。

誰でもここで何が起こっているのか分かりますか？

Answer 1

あなたがリンクしている式のように、畳み込みマイナスからプラスの無限大まで。有限シーケンスの場合は、どうやら不可避のオーケールになる境界効果に対処する必要があります。numpyのとそうするさまざまな方法を提供しますscipyのダウンロード：

numpy convolve：

モード：{ 'フル'、 '有効'、 '同じ'}、オプションの

scipy convolve：

モード：{'反映'、 '定数'、 '最も近い'、 'ミラー'、 'ラップ'}、オプション

次のポイントは、あなたの起源を置く場所です。指定した実装では、信号をt=0から開始し、負の値の加算を破棄します。t。 Scipyはそのためにパラメータoriginを提供しています。

origin：array_like、オプション originパラメータは、フィルタの配置を制御します。デフォルトは0です。

であるあなたは、実際にscipyのダウンロードconvolveを使用して実装の動作を模倣することができます

あなたはどのように決定する必要があり畳み込みをやって、それをまとめる

from scipy.ndimage.filters import convolve as convolve_sci 
from pylab import * 

N = 100 
start=N//8 
end = N-start 
A = zeros(N) 
A[start:end] = 1 
B = zeros(N) 
B[start:end] = linspace(1,0,end-start) 

figure(figsize=(6,7)) 
subplot(411); grid(); title('Signals') 
plot(A) 
plot(B) 
subplot(412); grid(); title('A*B numpy') 
plot(convolve(A,B, mode='same')) 
subplot(413); grid(); title('A*B scipy (zero padding and moved origin)') 
plot(convolve_sci(A,B, mode='constant', origin=-N//2)) 
tight_layout() 
show() 

シーケンスの外側の値を処理します（ゼロに設定する（numpy）、リフレクション、ラッピングなど）。シグナルの起点はどこに置くか。

もnumpyののデフォルトをそのに気づくと、彼らは境界効果（ゼロパディング対反映）をどのように扱うか異なりscipyのダウンロードしてください。

Answer 2

output[i] += signal[i - j] * ir[j]が、それは単純に理解できるようにすることドキュメントの変数名を使用してoutput[i] += signal[j] * ir[i - j]

する必要があり、ドキュメントの表記に合わせて、まず：あなたはこのためにコンボリューションで逃げる

i = len(signal) 
for n in range(i): 
    for k in range(n): 
     output[n] += signal[k] * ir[n - k] 

したがって、f*g == g*f（あなたのダイアグラムを参照）

主な違いは、長さがmシグナルと長さが0の「基本」コンボリューションインパルスは長さm + n -1（np.convolve docsを参照）ですが、np.convolve(. . . , mode = 'same')とscipy.ndimage.convolve1dの両方は、信号の両端から要素をトリミングすることによって長さを返します（m）。

だからあなたの問題はあなただけnp.convolve(..., mode = 'same')と同じトリミングを行うために、なぜ

np.all(
     np.convolve(single_square_wave, single_saw_wave)[:len(single_square_wave)]\ 
     ==\ 
     convolve1d(single_square_wave, single_saw_wave) 
     ) 

True 

をある右からあなたの信号をトリミングしているである、あなたが必要があるだろう：

def convolve1d_(signal, ir): 
    """ 
    we use the 'same'/'constant' method for zero padding. 
    """ 
    pad = len(ir)//2 - 1 
    n_ = range(pad, pad + len(signal)) 
    output = np.zeros(pad + len(signal)) 

    for n in n_: 
     kmin = max(0, n - len(ir) + 1) 
     kmax = min(len(ir), n) 
     for k in range(kmin, kmax): 
      output[n] += signal[k] * ir[n - k] 

    return output[pad:] 

テスト：

np.all(
     np.convolve(single_square_wave, single_saw_wave, mode = 'same')\ 
     ==\ 
     convolve1d_(single_square_wave,single_saw_wave) 
     ) 

True