numpyでconvolveしてみる

numpyを使って数値計算で畳み込みをしてみたのでメモしておきます。

numpyで畳み込みするにはnumpy.convolveという関数を使用します。
\[f(x) = x^2(0\le x \le 2)\]と幅2の長方形の畳み込みをしてみます。

import numpy as np
import matplotlib.pyplot as plt
from math import *

x = np.arange(0,2.01,0.01) 
f = x**2  //f(x)の配列を作る
g = np.ones(200)/(len(np.ones(200))) //g(x)の配列を作る
con = np.convolve(f,g) //畳み込み
conx = np.arange(-1,3,0.01)
plt.plot(conx,con)
plt.show()

まずf(x)の配列を作ります。

x = np.arange(0,2.01,0.01) 
f = x**2  //f(x)の配列を作る

ここでは0~2まで0.01刻みの配列にしています。
次にg(x)の配列を作ります。

g = np.ones(200)/(len(np.ones(200))) //g(x)の配列を作る

np.convolve関数は配列同士を片方逆順にして、かけ算して足し合わせていくというもの(([1,2,3],[0,1,0.5])の配列をnp.convolveすると([1*0,2*0+1*1,3*0+2*1+1*0.5,1*3+2*0.5,3*0.5]) = ([0,1,2.5,4,1.5])という配列ができる)なので、畳み込みにするにはどちらかを規格化しておく必要があります。
ここではg(x)を規格化して配列の和が1になるようにしています。

con = np.convolve(f,g) //畳み込み
conx = np.arange(-1,3,0.01)
plt.plot(conx,con)
plt.show()

こうして作った配列をnp.convolveで畳み込みます。
np.convolve(f,g)の配列の長さはlen(f)+len(g)-1となります。また0≦x≦2のf(x)に幅2の長方形g(x)を畳み込むので出てくる関数は-1≦x≦3です。したがってそのような配列np.arange(-1,3,0.01)を作っています。これですべて材料が揃ったので後はプロットするだけです。

出てくるグラフは次のようになります。