畳み込みと相関の基礎:理論・畳み込み定理・Python実装

信号処理の基礎である畳み込みと相関を、連続/離散の定義から畳み込み定理、FFTを用いた高速畳み込み、移動平均との関係まで体系的に解説し、NumPy/SciPyによるPython実装で実証します。

はじめに

信号処理の世界では、**畳み込み(convolution)相関(correlation)**が根本的な役割を担っています。フィルタリング、特徴抽出、パターンマッチング、システム同定など、応用範囲は極めて広く、FFTウェーブレット変換といった周波数解析手法もこれらの概念を土台にしています。

本記事では、連続/離散の畳み込みと相関の定義を確認し、両者の数学的関係を整理します。さらに畳み込み定理を通じて時間領域と周波数領域がどのように接続されるかを示し、np.convolvescipy.signal.correlate および FFT による高速畳み込みを Python で実装・比較します。最後に移動平均フィルタが「ボックス窓との畳み込み」として表現されることを確認します。

畳み込みの定義

連続時間の畳み込み

2つの連続関数 \(f(t)\) と \(g(t)\) の畳み込み \((f * g)(t)\) は次のように定義されます。

\[(f * g)(t) = \int_{-\infty}^{\infty} f(\tau) \cdot g(t - \tau) \, d\tau \tag{1}\]

直感的には「\(g\) を時間反転して \(t\) だけずらし、\(f\) と重ね合わせて積分する」操作です。線形時不変(LTI)システムの応答が入力とインパルス応答の畳み込みで表現されるため、フィルタ処理は本質的に畳み込みそのものといえます。

離散時間の畳み込み

長さが有限あるいは可算無限の離散信号 \(x[n]\) と \(h[n]\) について、畳み込みは次のように定義されます。

\[y[n] = (x * h)[n] = \sum_{k=-\infty}^{\infty} x[k] \cdot h[n - k] \tag{2}\]

実装上は \(h\) を時間反転(\(h[n-k]\) )してスライドさせ、対応する \(x[k]\) と要素積を取り、和をとります。\(x\) の長さが \(N\) 、\(h\) の長さが \(M\) の場合、出力長は \(N + M - 1\) となります。

主要な性質

畳み込みは以下の代数的性質を満たします。

  • 交換律: \(x * h = h * x\)
  • 結合律: \((x * h) * g = x * (h * g)\)
  • 分配律: \(x * (h + g) = x * h + x * g\)

これにより、複数のフィルタを直列に適用する場合、それらを事前に畳み込んで1つの等価フィルタにまとめられます。

相関の定義

相互相関(Cross-Correlation)

2つの信号 \(x[n]\) と \(y[n]\) の相互相関は、ラグ \(\ell\) を引数として次のように定義されます。

\[R_{xy}[\ell] = \sum_{n=-\infty}^{\infty} x[n] \cdot y[n + \ell] \tag{3}\]

相互相関は「\(y\) を \(\ell\) だけずらしたときに \(x\) とどれだけ似ているか」を示し、信号間の遅延推定やテンプレートマッチングに用いられます。

自己相関(Auto-Correlation)

\(y = x\) とした特殊な場合が自己相関で、信号自身の周期性や類似ラグを検出するのに使われます。

\[R_{xx}[\ell] = \sum_{n=-\infty}^{\infty} x[n] \cdot x[n + \ell] \tag{4}\]

\(R_{xx}[0]\) は信号の総エネルギー(または分散)に対応し、自己相関は \(\ell = 0\) で最大値をとります。

畳み込みとの関係

相互相関と畳み込みは第二引数を時間反転するか否かだけが異なります。実数信号では次の関係が成り立ちます。

\[R_{xy}[\ell] = (x * \tilde{y})[\ell], \quad \tilde{y}[n] = y[-n] \tag{5}\]

つまり、相関は「時間反転した信号との畳み込み」として実装できます。複素信号の場合は時間反転に加えて複素共役をとります(\(\tilde{y}[n] = y^*[-n]\) )。

畳み込み定理

時間領域畳み込み ↔ 周波数領域積

畳み込みの最も強力な性質は、周波数領域では単純な積になることです。離散フーリエ変換(DFT)を \(\mathcal{F}\) で表すと、

\[\mathcal{F}\{x * h\}[k] = X[k] \cdot H[k] \tag{6}\]

逆に、時間領域の積は周波数領域の畳み込みに対応します(窓関数とPSDの記事で扱ったスペクトル漏れの根拠)。

\[\mathcal{F}\{x \cdot h\}[k] = \frac{1}{N}(X * H)[k] \tag{7}\]

高速畳み込みの計算量

長さ \(N\) 同士の素朴な畳み込みは \(O(N^2)\) ですが、畳み込み定理と \(O(N \log N)\) のFFTを組み合わせると、

\[y = \mathcal{F}^{-1}\{ \mathcal{F}\{x\} \cdot \mathcal{F}\{h\} \} \tag{8}\]

として \(O(N \log N)\) で計算できます。長いカーネルでは数桁の高速化が得られます。なお、線形畳み込みを得るには循環畳み込みを避けるため、両信号を \(N + M - 1\) 以上の長さへゼロパディングする必要があります。

相関定理

相関にも同様の定理が成り立ちます。

\[\mathcal{F}\{R_{xy}\}[k] = X^*[k] \cdot Y[k] \tag{9}\]

特に自己相関のフーリエ変換はパワースペクトル密度に等しく(Wiener–Khinchin の定理)、PSD 推定の理論的基盤となっています。

Python 実装

np.convolvescipy.signal.correlate

最初に基本的なAPIで挙動を確認します。

import numpy as np
from scipy.signal import correlate

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
h = np.array([1.0, 0.5, 0.25])

# 線形畳み込み(mode='full' は出力長 N+M-1)
y_conv = np.convolve(x, h, mode='full')
print('convolve:', y_conv)

# 相互相関(mode='full')
y_corr = correlate(x, h, mode='full')
print('correlate:', y_corr)

# 式 (5) の確認: 相関 = 反転した h との畳み込み
y_corr_via_conv = np.convolve(x, h[::-1], mode='full')
print('correlate via flipped conv:', y_corr_via_conv)

correlate(x, h)convolve(x, h[::-1]) が一致することで、式 \((5)\) が実証されます。

FFT による高速畳み込み

長いカーネルで np.convolve と FFT ベースの実装の速度を比較します。

import numpy as np
import time
from scipy.signal import fftconvolve

np.random.seed(0)
N = 1 << 16   # 65536
M = 1 << 12   # 4096
x = np.random.randn(N)
h = np.random.randn(M)

# 直接畳み込み
t0 = time.perf_counter()
y_direct = np.convolve(x, h, mode='full')
t_direct = time.perf_counter() - t0

# FFT 畳み込み(SciPy)
t0 = time.perf_counter()
y_fft = fftconvolve(x, h, mode='full')
t_fft = time.perf_counter() - t0

# 自前 FFT 畳み込み(式 (8) の実装)
L = N + M - 1
n_fft = 1 << int(np.ceil(np.log2(L)))
t0 = time.perf_counter()
X = np.fft.rfft(x, n=n_fft)
H = np.fft.rfft(h, n=n_fft)
y_manual = np.fft.irfft(X * H, n=n_fft)[:L]
t_manual = time.perf_counter() - t0

print(f'direct  : {t_direct*1e3:7.2f} ms')
print(f'fftconv : {t_fft*1e3:7.2f} ms')
print(f'manual  : {t_manual*1e3:7.2f} ms')
print(f'max err : {np.max(np.abs(y_direct - y_fft)):.2e}')

実行例では、\(N = 65536\) 、\(M = 4096\) で np.convolve が数百 ms かかるのに対し、fftconvolve は十数 ms 程度で完了し、結果の最大誤差は浮動小数点誤差レベル(\(10^{-10}\) オーダー)です。

移動平均との接続

長さ \(L\) の単純移動平均(SMA)は、係数がすべて \(1/L\) のボックス窓 \(h_{\text{box}}[n]\) との畳み込みに等価です。

\[h_{\text{box}}[n] = \frac{1}{L}, \quad 0 \leq n \leq L - 1 \tag{10}\] \[y[n] = \frac{1}{L}\sum_{k=0}^{L-1} x[n - k] = (x * h_{\text{box}})[n] \tag{11}\]
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1)
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2*np.pi*5*t) + 0.4*np.random.randn(len(t))

L = 21
h_box = np.ones(L) / L

# 移動平均 = ボックス窓との畳み込み
sma = np.convolve(signal, h_box, mode='same')

# 検証: scipy の uniform_filter1d と一致
from scipy.ndimage import uniform_filter1d
sma_ref = uniform_filter1d(signal, size=L, mode='nearest')

plt.figure(figsize=(10, 4))
plt.plot(t, signal, alpha=0.4, label='Noisy signal')
plt.plot(t, sma, label=f'Convolution with box (L={L})')
plt.plot(t, sma_ref, '--', label='uniform_filter1d (reference)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

両者の出力は端効果の処理を除いて一致し、移動平均が畳み込みの特殊例であることが確認できます。フィルタ係数を変えれば(例: 三角窓・Hann 窓)、自然と FIR フィルタへ拡張できます。

まとめ

  • 畳み込みは「片方を反転してずらし、積和をとる」操作であり、LTI システムの応答を完全に記述する
  • 相関は時間反転を伴わない積和で、相互相関 \(R_{xy}\) はテンプレートマッチング、自己相関 \(R_{xx}\) は周期検出に用いる
  • 畳み込み定理(式 \((6)\) )により時間領域畳み込みは周波数領域の積になり、FFT で \(O(N \log N)\) の高速畳み込みが可能になる
  • 移動平均はボックス窓との畳み込みであり、FIR フィルタの最も単純なケースとして位置付けられる
  • Wiener–Khinchin の定理により、自己相関と PSD は同一概念の異なる表現である

これらの基礎を押さえることで、フィルタ設計・スペクトル解析・時間-周波数解析が一貫した数学体系として理解できるようになります。

関連記事

参考文献

  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • Smith, S. W. (1997). The Scientist and Engineer’s Guide to Digital Signal Processing. California Technical Publishing.
  • NumPy convolve documentation
  • SciPy signal documentation