人気ブログランキング | 話題のタグを見る

畳み込み(convolution)

Haskell で畳み込み(convolution) の計算をしてみた。これはWikipediaでは次のように定義されている。

畳み込み(たたみこみ、convolution)とは関数fを平行移動しながら関数gを重ね足し合わせる二項演算である。畳み込み積分・合成積・重畳積分とも呼ばれる。(Wikipedia より)

これを読んでも分かりにくいが、たとえば、電気回路のアンプに様々な波形の信号を入れると、アンプの特性によっていろいろな波形が出力に現れてくる。このとき、このアンプにインパルス信号をを入力したときの波形(インパルス応答)が分かっていれば、あとは入力波形にどんなものを持ってきても、インパルス応答と入力波形のたたみ込みを計算すると、出力の波形が計算できてしまうというすぐれものだ。

説明はともあれ、プログラムは次のようになる。

畳み込み計算プログラム(convolution.hs)
conv :: Num a => [a] -> [a] -> [a]
conv hs xs =
  let hs' = reverse hs
      hn = length hs
      xn = length xs
      xs' = replicate (hn-1) 0 ++ xs
  in take xn $ conv' hs' xs'

conv' :: Num a => [a] -> [a] -> [a]
conv' hs zs@(x:xs) = (sum $ zipWith (*) hs zs) : conv' hs xs

たったこれだけだ。これを Hugs でテストしてみよう。(上のプログラムをHugsで実行する手順の紹介は省略する。)
Main> conv [1,2] [9,2,3,6]
[9,20,7,12]

インパルス応答 [1,2] と入力信号[9,2,3,6] との畳み込みは次のようになる。
入力9の前の入力は0なので最初の出力は9*1+0*2=9になる。
次の出力は二番目の入力2にインパルス応答の1をかけたものと最初の入力9にインパルス応答2の重みをかけたものの和になるので、2*1+9*2=20 だ、同様に3番目の出力は3*1+2*2=7... と計算できる。

それでは電気回路のシミュレーションをやってみよう。インパルス応答は指数関数的に減少するハイパスフィルターのものを用いる。入力としては、矩形波のパルスを入力する。次のプログラムを作成し、Hugs のコンソールから main と入力すると、sample.csv というファイルができるので、それをExcel で取り込んでグラフを作ると、立ち上がりと、立下りにカーブのある出力信号の形を見ることができる。

main = writeData "sample.csv" (conv impulseResp stepInput)

conv :: Num a => [a] -> [a] -> [a]
conv hs xs =
  let hs' = reverse hs
      hn = length hs
      xn = length xs
      xs' = replicate (hn-1) 0 ++ xs
  in take xn $ conv' hs' xs'

conv' :: Num a => [a] -> [a] -> [a]
conv' hs zs@(x:xs) = (sum $ zipWith (*) hs zs) : conv' hs xs

writeData fname xs = writeFile fname datastr
  where datastr = unlines $ map show xs

impulseResp = map (\x -> exp (-1*x)) [0..20]

stepInput = replicate 10 1.0 ++ replicate 30 0.0

予備知識が少し必要な話になってしまったので、Haskell のプログラミングの紹介から外れてしまったかもしれないが、結構面倒なことを簡潔に表現できる Haskell プログラムの表現力には驚かされる。

ちょっとテクニカルな話になるが、convolution を計算していく時に計算を終了する場所をどこにするかということは結構悩ましい問題だったりする。Haskell では遅延評価ができるので、終了条件のことは放っておいて take 関数で必要な所だけ取り出せるのでプログラムの煩雑さが激減する。速度の問題はどうあれ、ちょっとした数値計算のシミュレーションには Haskell は非常に便利に感じる。
by tnomura9 | 2009-09-22 18:46 | Haskell | Comments(0)
<< 改行文字変換プログラム Haskell でサインカーブ... >>