ディジタル信号処理を勉強するひとは、すべからくマスターすべき演算が畳み込みである（私見）。 畳み込み演算の重要性は論を俟たない。いわゆる線形時不変システムはインパルス応答と入力信号との畳み込みにより記述される。音声のディジタル信号処理の文脈でいえば、分析合成フィルタ（いわゆるボコーダー）の実現に畳み込み演算が用いられる。
Pythonライブラリのnumpyには畳み込みを実行する関数convolveがあり、またscipyにもconvolveがあるため、実用上はそれらを使えば良い（というかぜひ使うべきである）。
ただスクラッチで畳み込みを実装せよと言われたとき、はたと手が止まる人も多いのではなかろうか。
本記事ではディジタル信号処理を対象とし、FIRフィルタ（インパルス応答が有限長）としての畳み込み演算の実装について説明する。
説明
インパルス応答とはデジタルシステムにインパルス信号を入力したときの出力を指す。ここで，入力信号を$x[n]$，インパルス応答を$h[n]$，出力信号を$y[n]$とする。$n$は時間インデックスである。出力信号$y[n]$は次式で示される。
$$ y[n] = \sum\limits_{m = -\infty}^{\infty} x[n-m]h[m] $$
ただし$x[n]$の長さを$N$とする。
$$ x[0], x[1], \ldots, x[N-1] $$
さらに$h[n]$の長さを${M}$とする。
$$ h[0], h[1], \ldots, h[M -1] $$
すると上式の総和は$m=0$から$M - 1$までとなり、$y[n]$の長さは$N + M -1$となる。
上記の説明より、$n = 0, 1, ..., N + M - 1$について
$$ y[n] = \sum\limits_{m = 0}^{M -1} x[n-m]h[m] $$
である。
$x[n]$は長さ$N$の条件を持っている。しかしながら、例えば上記の畳み込みの式において、$n - m $が $N$以上になることが起こり得る。一方で、実際のデータとしては
$$ x[0], x[1], \ldots, x[N-1] $$
しか手元にない。そのとき、$x[n - m] = 0$としたいので、事前に$x[n]$を拡張し、長さ$N + M -1$の系列となるよう、$x[n]$の右側に0を詰めておくのがスマートである。
\begin{align} \tilde{x}[n] = \begin{cases} x[n] & (0 \leq n \leq N-1)\\ 0 & (N\leq n \leq N + M - 2) \end{cases} \end{align}
この結果として、$\tilde{x}[n]$が長さ$N + M -1$の系列になる。 さらに$\tilde{x}[n-m]$のインデックスが非負になる$n \geq m $の範囲で和を取ればよい。 つまり、 $$y[n] = \sum\limits_{m = 0}^{n} \tilde{x}[n-m]h[m]$$ である。 このとき、和のインデックス上限$n$と$h[m]$のインデックス上限$M - 1$の兼ね合いで、$x[n]$と類似の問題が起きる。この問題の煩わしさを避けるため、$h[n]$についても同様の拡張を行う。
\begin{align} \tilde{h}[n] = \begin{cases} h[n] & (0 \leq n \leq N-1)\\ 0 & (N\leq n \leq N + M - 2) \end{cases} \end{align}
最終的に畳み込みは以下のように書ける。 $$ y[n] = \sum\limits_{m = 0}^{n} \tilde{x}[n-m]\tilde{h}[m] \; \; \; (0 \leq n \leq N + M -1) $$
実装
素直にPythonコードに落とし込めばよい。入力の拡張としてのゼロ詰めにはnumpyのzeros関数を使う。
地のfor文を使うという、ナイーブな実装なのでPython的には（とても）遅く、あまり実用的ではないのは承知のうえで。
# 前提 # 入力 x[0] ... x[N-1] # インパルス応答 h[0], ... h[M-1] import numpy as np # 出力 y[0] ... y[M+N-2] y = np.zeros(len(h) + len(x) - 1, dtype=np.float32) # ゼロづめによる拡張 hzero = np.hstack([h, np.zeros(len(x) - 1)]) xzero = np.hstack([x, np.zeros(len(h) - 1)]) for n in range(0, len(y)): for k in range(0, n + 1): y[n] = y[n] + hzero[k] * xzero[n - k]
プロットの実例はそのうち追記予定（ほんまか）。
参考
これを読もう。 www.ic.is.tohoku.ac.jp
畳み込みと相互相関の違いとか。 qiita.com
