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

はじめに

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

本記事では、連続/離散の畳み込みと相関の定義を確認し、両者の数学的関係を整理します。さらに畳み込み定理を通じて時間領域と周波数領域がどのように接続されるかを示し、np.convolve・scipy.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]\) ）。

マッチドフィルタとの関係

なぜ相関演算が「マッチドフィルタ（matched filter）」と呼ばれる検出器の実体になるのかを、雑音下の最適検出問題として導出します。既知の波形（テンプレート）\(s[n]\) （\(n = 0, \dots, K-1\) ）が加法性白色雑音 \(w[n]\) （平均0、分散 \(\sigma^2\) 、各サンプル独立）に埋もれて受信されるとします。

\[ r[n] = s[n] + w[n] \]

この \(r[n]\) を長さ \(K\) のFIRフィルタ \(h[n]\) に通した出力 \(y[n] = (r * h)[n]\) のうち、テンプレートがちょうど窓内に収まる時刻 \(n = K-1\) でのサンプルに着目すると、

\[ y[K-1] = \sum_{k=0}^{K-1} r[k]\, h[K-1-k] = \underbrace{\sum_{k=0}^{K-1} s[k]\, h[K-1-k]}_{\text{信号成分 } S} + \underbrace{\sum_{k=0}^{K-1} w[k]\, h[K-1-k]}_{\text{雑音成分 } W} \]

信号成分 \(S\) は決定論的な定数であり、雑音成分 \(W\) は分散 \(\text{Var}[W] = \sigma^2 \sum_{k=0}^{K-1} h[K-1-k]^2\) を持つ確率変数です（\(w[k]\) が独立なため分散が線形に加算されます）。出力のSN比は

\[ \text{SNR}(h) = \frac{S^2}{\text{Var}[W]} = \frac{\left(\sum_{k=0}^{K-1} s[k]\, h[K-1-k]\right)^{\!2}}{\sigma^2 \sum_{k=0}^{K-1} h[K-1-k]^2} \]

で与えられます。分子にコーシー・シュワルツの不等式を適用すると、

\[ S^2 = \left(\sum_{k=0}^{K-1} s[k]\, h[K-1-k]\right)^{\!2} \le \left(\sum_{k=0}^{K-1} s[k]^2\right)\left(\sum_{k=0}^{K-1} h[K-1-k]^2\right) \]

であり、等号は \(h[K-1-k] \propto s[k]\) 、すなわち

\[ h[n] = c \cdot s[K-1-n] \tag{6} \]

（\(c\) は任意定数）のときにのみ成立します。このとき \(\text{SNR} = \sum_k s[k]^2 / \sigma^2\) （信号エネルギー÷雑音分散）が達成され、これが線形フィルタの中で到達可能な最大値です。式\((6)\) は「テンプレートを時間反転したフィルタ」そのものであり、式\((5)\) の関係に照らせば、この \(h\) を用いた畳み込み \((r * h)\) は \(r\) とテンプレート \(s\) の相互相関 \(R_{rs}\) に一致します。つまり「未知の到来時刻を持つ既知波形をSN比最大で検出する最適フィルタ」は、テンプレートとの相関演算そのものであり、これが相関がマッチドフィルタとして機能する数学的根拠です。レーダー・ソナー・GPS信号捕捉など「既知波形の検出」を要するあらゆる場面でこの結果が使われています。

畳み込み定理

畳み込み定理の証明（連続時間）

畳み込みの最も強力な性質は、時間領域の畳み込みが周波数領域では単純な積になることです。まず連続時間で、この事実を定義から直接証明します。フーリエ変換を

\[ \mathcal{F}\{f\}(\omega) = F(\omega) = \int_{-\infty}^{\infty} f(t)\, e^{-j\omega t}\, dt \]

と定義します。式\((1)\) の畳み込みにこの定義を適用すると、

\[ \mathcal{F}\{f * g\}(\omega) = \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} f(\tau)\, g(t - \tau)\, d\tau \right) e^{-j\omega t}\, dt \]

被積分関数が可積分であれば（\(f, g \in L^1(\mathbb{R})\) ）フビニの定理により積分順序を交換でき、

\[ \mathcal{F}\{f * g\}(\omega) = \int_{-\infty}^{\infty} f(\tau) \left( \int_{-\infty}^{\infty} g(t - \tau)\, e^{-j\omega t}\, dt \right) d\tau \]

を得ます。内側の積分に変数変換 \(u = t - \tau\) （すなわち \(t = u + \tau\) 、\(dt = du\) 、積分範囲は \(\tau\) に依らず \((-\infty, \infty)\) のまま）を施すと、

\[ \int_{-\infty}^{\infty} g(t - \tau)\, e^{-j\omega t}\, dt = \int_{-\infty}^{\infty} g(u)\, e^{-j\omega(u + \tau)}\, du = e^{-j\omega \tau} \int_{-\infty}^{\infty} g(u)\, e^{-j\omega u}\, du = e^{-j\omega \tau}\, G(\omega) \]

となります（\(e^{-j\omega\tau}\) は \(u\) に依存しないため積分の外に出せます）。これを代入すれば、

\[ \mathcal{F}\{f * g\}(\omega) = \int_{-\infty}^{\infty} f(\tau)\, e^{-j\omega \tau}\, G(\omega)\, d\tau = G(\omega) \int_{-\infty}^{\infty} f(\tau)\, e^{-j\omega \tau}\, d\tau = F(\omega)\, G(\omega) \]

と、目的の関係 \(\mathcal{F}\{f * g\} = F \cdot G\) が畳み込みとフーリエ変換の定義のみから導かれます。鍵となるのは「時間反転・シフトされた \(g(t-\tau)\) 」を変数変換で通常の \(g(u)\) に戻す操作であり、これによりフーリエ核 \(e^{-j\omega\tau}\) を積分の外に括り出せる点にあります。

離散フーリエ変換における畳み込み定理と巡回畳み込み

離散信号でも全く同じ議論が成り立ちますが、変数変換が「\(\bmod N\) の巡回シフト」になる点が本質的に異なります。長さ \(N\) の離散フーリエ変換を

\[ X[k] = \sum_{n=0}^{N-1} x[n]\, e^{-j2\pi kn/N} \]

と定義し、長さ \(N\) の2信号 \(x[n], h[n]\) に対する**巡回畳み込み（circular convolution）**を

\[ (x \circledast h)[n] = \sum_{m=0}^{N-1} x[m]\, h[(n - m) \bmod N] \tag{7} \]

と定義します（線形畳み込み式\((2)\) とは異なり、添字が \(N\) を法として折り返される点に注意）。このDFTを取ると、

\[ \mathcal{F}\{x \circledast h\}[k] = \sum_{n=0}^{N-1} \left( \sum_{m=0}^{N-1} x[m]\, h[(n-m) \bmod N] \right) e^{-j2\pi kn/N} \]

和の順序を交換し、内側の和に変数変換 \(p = (n - m) \bmod N\) （したがって \(n \equiv m + p \pmod N\) ）を施します。複素指数 \(e^{-j2\pi kn/N}\) は \(n\) について周期 \(N\) であるため、\(n\) を \(m+p\) の代表元のどれで取っても値は変わらず、

\[ \mathcal{F}\{x \circledast h\}[k] = \sum_{m=0}^{N-1} x[m]\, e^{-j2\pi km/N} \sum_{p=0}^{N-1} h[p]\, e^{-j2\pi kp/N} = X[k]\, H[k] \tag{8} \]

が得られます。連続時間の変数変換 \(u=t-\tau\) が「並進」だったのに対し、離散DFTでは「\(\bmod N\) の巡回シフト」に置き換わる点が唯一の違いです。重要なのは、この定理が成り立つのは巡回畳み込み \(\circledast\) に対してであり、式\((2)\) の線形畳み込みそのものではない、という点です。両者の違いと、線形畳み込みをFFTで計算する際にゼロパディングが必要な理由を後の節で見ます。

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

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

直接畳み込みの計算量：\(O(NM)\) の導出

長さ \(N\) の \(x\) と長さ \(M\) の \(h\) の線形畳み込み \(y[n] = \sum_k x[k] h[n-k]\) （式\((2)\) ）を素朴に計算する場合に必要な乗算回数を正確に数えます。出力の各点 \(y[n]\) の計算で有効な項は、\(k \in [0, N-1]\) かつ \(n-k \in [0, M-1]\) を満たす \(k\) についてのみ存在します。\(m = n - k\) とおくと、これは「\(k \in [0,N-1]\) 、\(m \in [0,M-1]\) の組 \((k,m)\) すべてに対し、出力インデックス \(n = k+m\) に \(x[k] h[m]\) を1回加算する」という操作に対応します。\((k,m)\) の組は重複なくちょうど \(N \times M\) 通り存在するため、乗算回数は厳密に \(NM\) 回であり、これが \(O(NM)\) （\(N \sim M\) のとき \(O(N^2)\) ）という計算量の由来です。

FFTによる高速化：\(O((N+M)\log(N+M))\)

畳み込み定理（式\((8)\) ）と、 Cooley–Tukeyアルゴリズムによる \(O(N\log N)\) のFFT を組み合わせると、線形畳み込みは次の手順で計算できます。

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

必要なのは長さ \(N_{\text{fft}}\) の順方向FFT2回・要素積1回（\(O(N_{\text{fft}})\) ）・逆FFT1回であり、合計で \(O(N_{\text{fft}} \log N_{\text{fft}})\) です。\(N_{\text{fft}}\) は次節で導出するように \(N+M-1\) 程度に取れば十分なので、全体の計算量は \(O((N+M)\log(N+M))\) となります。カーネルが長い（\(M\) が \(N\) に対して無視できない）場合、直接法の \(O(NM)\) に比べて数桁の高速化が得られます。

線形畳み込みと巡回畳み込みの違い、ゼロパディング

式\((8)\) の畳み込み定理はDFT長 \(N_{\text{fft}}\) 上の巡回畳み込みに対して成り立つため、FFTでそのまま計算すると得られるのは巡回畳み込みであって、式\((2)\) の線形畳み込みではありません。両者を一致させる条件を導出します。

長さ \(N\) の \(x\) と長さ \(M\) の \(h\) の線形畳み込み \(y = x * h\) は、非零区間が \(n = 0, \dots, N+M-2\) （長さ \(L = N+M-1\) ）に限られます。\(x, h\) をどちらも長さ \(N_{\text{fft}} \ge L\) になるよう末尾にゼロを詰めてから式\((7)\) の巡回畳み込みを取ると、

\[ (x \circledast h)[n] = \sum_{m=0}^{N_{\text{fft}}-1} x[m]\, h[(n-m) \bmod N_{\text{fft}}] \]

において、\(x[m] \neq 0\) となるのは \(m \in [0, N-1]\) 、\(h[\cdot] \neq 0\) となるのは引数が \([0, M-1]\) のときに限られます。\(N_{\text{fft}} \ge L = N+M-1\) であれば、\(n \in [0, N_{\text{fft}}-1]\) の範囲で \(n - m\) が負になったり \(N_{\text{fft}}\) を超えて折り返したりすることがなく（\(0 \le n-m \le N_{\text{fft}}-1\) が保証される）、巡回畳み込みは線形畳み込みと完全に一致します（\(n \ge L\) の部分はゼロ）。

一方、\(N_{\text{fft}} < L\) の場合、線形畳み込みの本来 \(n \ge N_{\text{fft}}\) に現れるはずの値が \(n \bmod N_{\text{fft}}\) の位置に折り返して加算され、エイリアシング（折り返し誤差）として観測されます。特に影響を受けるのは先頭 \(L - N_{\text{fft}}\) サンプルであり、\(N_{\text{fft}} = N\) （ゼロパディングなし）の場合は先頭 \(M-1\) サンプルが汚染されます。したがって、FFTで線形畳み込みを正しく計算するには、\(N_{\text{fft}} \ge N + M - 1\) を満たすようゼロパディングする必要があります。実装上は基数2 FFTの都合上 \(N_{\text{fft}}\) を2のべき乗に切り上げることが多いですが、scipy.fft.next_fast_len を使えば2・3・5・7の積で表せる「高速な合成数長」の中から最小のものを選べるため、常に2のべき乗にするより無駄なパディングを減らせます。この現象は後のPython実装で数値的に確認します。

相関定理

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

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

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

Python 実装

`np.convolve` と `scipy.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 畳み込み（式 (10) の実装）
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}\) オーダー）です。この誤差の小ささ自体が、式\((8)\) の畳み込み定理（時間領域畳み込み＝周波数領域の積の逆変換）が数値的にも成立していることの直接的な検証になっています。

巡回畳み込みとエイリアシングの数値実験

前節で導出した「ゼロパディングが不十分だと巡回畳み込みは線形畳み込みとずれる」という主張を、実際にFFTで確認します。

import numpy as np

x = np.array([1., 2., 3., 4., 5., 6., 7., 8.])   # N=8
h = np.array([1., 1., 1., 1.])                    # M=4

# 真の線形畳み込み（長さ N+M-1=11）
y_linear = np.convolve(x, h, mode='full')

# 巡回畳み込み：ゼロパディングなし（FFT長=N=8、不十分）
Xf = np.fft.fft(x, n=8)
Hf = np.fft.fft(h, n=8)
y_circular_noP = np.real(np.fft.ifft(Xf * Hf))

# 巡回畳み込み：十分な長さへゼロパディング（N+M-1=11 以上、16に切り上げ）
n_fft = 16
Xf2 = np.fft.fft(x, n=n_fft)
Hf2 = np.fft.fft(h, n=n_fft)
y_circular_pad = np.real(np.fft.ifft(Xf2 * Hf2))[:11]

print('linear conv        :', np.round(y_linear, 3))
print('circular (n_fft=8) :', np.round(y_circular_noP, 3))
print('circular (n_fft=16):', np.round(y_circular_pad, 3))
print('max diff, no padding:', np.max(np.abs(y_circular_noP - y_linear[:8])))
print('max diff, padded    :', np.max(np.abs(y_circular_pad - y_linear)))

実行結果は次の通りです。

linear conv        : [ 1.  3.  6. 10. 14. 18. 22. 26. 21. 15.  8.]
circular (n_fft=8) : [22. 18. 14. 10. 14. 18. 22. 26.]
circular (n_fft=16): [ 1.  3.  6. 10. 14. 18. 22. 26. 21. 15.  8.]
max diff, no padding: 21.0
max diff, padded    : 1.7763568394002505e-15

ゼロパディングなし（n_fft=8）では、線形畳み込みの末尾3サンプル（インデックス8, 9, 10の値 21, 15, 8）が先頭3サンプル（インデックス0, 1, 2）に折り返して加算されており（例えば 22 = 1 + 21、18 = 3 + 15、14 = 6 + 8）、最大誤差は21.0にも達します。一方、\(N_{\text{fft}}=16 \ge L=11\) を満たすようゼロパディングすると、巡回畳み込みは浮動小数点誤差レベル（\(\sim 10^{-15}\) ）で線形畳み込みと一致します。この結果を図示すると次のようになります。

線形畳み込みと巡回畳み込みの比較。左：ゼロパディングなし（FFT長8）では末尾3サンプルが先頭に折り返してエイリアシングが生じる。右：十分な長さ（FFT長16）にゼロパディングすると巡回畳み込みは線形畳み込みと完全に一致する

相関とマッチドフィルタの数値実験

式\((6)\) の結論——「未知の到来時刻の既知波形を検出する最適フィルタはテンプレートとの相関」——を、雑音に埋もれたパルスの検出で確認します。

import numpy as np
from scipy.signal import correlate

np.random.seed(7)
N = 2000
pulse = np.exp(-0.5 * (np.arange(-15, 16) / 4.0)**2)  # ガウスパルス（長さ31）
true_delay = 823

signal = 0.6 * np.random.randn(N)
signal[true_delay:true_delay + len(pulse)] += pulse

# マッチドフィルタ = テンプレートとの相関
mf_output = correlate(signal, pulse, mode='same')
raw_argmax = np.argmax(mf_output)

print('true pulse center index :', true_delay + len(pulse)//2)
print('matched filter argmax   :', raw_argmax)
print('SNR before (peak/std)   :', pulse.max() / signal.std())
print('SNR after MF (peak/std) :', mf_output.max() / mf_output.std())

true pulse center index : 838
matched filter argmax   : 839
SNR before (peak/std)   : 1.675499480245094
SNR after MF (peak/std) : 4.836095661428964

真のパルス中心（インデックス838）に対し、マッチドフィルタ出力が最大値をとる位置（839）はわずか1サンプルのずれで一致しています。さらに、生信号でのSN比（ピーク振幅÷信号全体の標準偏差）が約1.68であるのに対し、マッチドフィルタ通過後のSN比は約4.84まで改善しており、約2.9倍のSN比向上が得られています。これは式\((6)\) で導いた「コーシー・シュワルツの等号成立条件（\(h[n] \propto s[K-1-n]\) ）」がまさにこの相関演算に対応するためです。

直接畳み込みとFFT畳み込みの計算時間比較

\(O(NM)\) の直接法と \(O((N+M)\log(N+M))\) のFFT法の実測クロスオーバー点を、信号長を変えながら計測します。

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

np.random.seed(42)
Ns = [8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536]

for N in Ns:
    x = np.random.randn(N)
    h = np.random.randn(N)  # N=M の正方ケース
    reps = max(1, int(2000 / N))

    t0 = time.perf_counter()
    for _ in range(reps):
        y_direct = np.convolve(x, h, mode='full')
    t_direct = (time.perf_counter() - t0) / reps

    t0 = time.perf_counter()
    for _ in range(reps):
        y_fft = fftconvolve(x, h, mode='full')
    t_fft = (time.perf_counter() - t0) / reps

    print(f'{N:6d} | direct {t_direct*1e3:9.4f} ms | fft {t_fft*1e3:8.4f} ms | speedup {t_direct/t_fft:6.2f}x')

実測結果（NumPy 2.4 / SciPy 1.18 環境）は以下の通りです。

     N=M | direct [ms] |  fft [ms] | speedup
       8 |      0.0010 |    0.0357 |   0.03x
      16 |      0.0011 |    0.0270 |   0.04x
      32 |      0.0014 |    0.0266 |   0.05x
      64 |      0.0022 |    0.0282 |   0.08x
     128 |      0.0045 |    0.0292 |   0.15x
     256 |      0.0128 |    0.0335 |   0.38x
     512 |      0.0422 |    0.1506 |   0.28x
    1024 |      0.1620 |    0.0808 |   2.00x
    2048 |      0.6563 |    0.1754 |   3.74x
    4096 |      2.6386 |    0.1340 |  19.70x
    8192 |     10.7775 |    0.3876 |  27.80x
   16384 |     50.5015 |    0.7445 |  67.84x
   32768 |    207.4742 |    1.5122 | 137.20x
   65536 |    852.2690 |    2.8904 | 294.86x

\(N=M \le 512\) では直接法の方が速く（FFTの固定オーバーヘッドが支配的）、\(N=M=1024\) 付近でクロスオーバーが生じ、それ以降はFFT法が急速に優位になります。\(N=M=65536\) では約295倍の高速化が得られており、\(O(NM)\) と \(O(N\log N)\) の理論通りの傾向が実測でも確認できます。カーネルが短い（\(M \ll N\) ）実用上のFIRフィルタ設計ではこの分岐点はカーネル長 \(M\) にも依存するため、SciPyは scipy.signal.choose_conv_method により信号長・カーネル長・データ型から最適な方式を自動選択する機能を提供しています。

オーバーラップアド法・オーバーラップセーブ法

FFT畳み込みは高速ですが、信号長 \(N\) が非常に大きい（ストリーミング処理やメモリに乗り切らない長さ）場合、\(N_{\text{fft}} \approx N\) の巨大な1回のFFTを計算するのは非効率・非現実的です。そこで長い信号をブロックに分割し、短いカーネル \(h\) （長さ \(M\) ）との畳み込みをブロック単位でFFT計算し、あとから正しくつなぎ合わせる手法が使われます。

オーバーラップアド法（Overlap-Add, OLA）：入力を長さ \(L\) の重ならないブロックに分割し、各ブロックと \(h\) の畳み込み（長さ \(L+M-1\) ）をFFTで計算します。隣接ブロックの出力は \(M-1\) サンプルだけ重なるため、この重なり部分を加算して連結します。
オーバーラップセーブ法（Overlap-Save, OLS）：各ブロックの先頭に直前ブロックの末尾 \(M-1\) サンプルを付加してから長さ \(L\) の巡回畳み込みを計算し、先頭 \(M-1\) サンプル（エイリアシングで汚染される部分）を破棄して残りだけを採用します。加算が不要な代わりに、出力の一部を捨てる点が異なります。

どちらも1回あたりのFFT長を \(M\) 程度に抑えたまま、全体としては単一の巨大FFTと同じ結果が得られます。SciPyには scipy.signal.oaconvolve としてOLA法が実装されています。

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

np.random.seed(3)
N = 2_000_000   # 長い信号（ストリーミングを想定）
M = 65          # 短いFIRカーネル
x = np.random.randn(N)
h = np.random.randn(M)

t0 = time.perf_counter()
y_full_fft = fftconvolve(x, h, mode='full')
t_full = time.perf_counter() - t0

t0 = time.perf_counter()
y_oa = oaconvolve(x, h, mode='full')
t_oa = time.perf_counter() - t0

print('max abs diff:', np.max(np.abs(y_full_fft - y_oa)))
print(f'single-shot fftconvolve: {t_full*1e3:.1f} ms')
print(f'oaconvolve (block OLA) : {t_oa*1e3:.1f} ms')

max abs diff: 3.019806626980426e-14
single-shot fftconvolve: 66.5 ms
oaconvolve (block OLA) : 16.7 ms

両者の出力は浮動小数点誤差レベル（\(3.0 \times 10^{-14}\) ）で完全に一致しており、ブロック処理でも正しく線形畳み込みが再現できることが確認できます。さらに、\(N=2{,}000{,}000\) 、\(M=65\) というカーネルが信号に比べて極端に短いケースでは、単一の巨大FFT（66.5 ms）よりもブロック分割したOLA（16.7 ms）の方がおよそ4倍高速でした。これは、\(N_{\text{fft}} \approx N\) の1回の巨大FFTがメモリ帯域に律速されるのに対し、ブロックサイズを適切に選んだOLAはキャッシュに収まる小さなFFTを繰り返すことで実効効率を高められるためです。

移動平均との接続

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

\[h_{\text{box}}[n] = \frac{1}{L}, \quad 0 \leq n \leq L - 1 \tag{12}\] \[ y[n] = \frac{1}{L}\sum_{k=0}^{L-1} x[n - k] = (x * h_{\text{box}})[n] \tag{13} \]

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 フィルタへ拡張できます。

発展：深層学習における畳み込み定理の再評価（2023年以降）

畳み込み定理は100年近い歴史を持つ古典的な結果ですが、2023年以降のディープラーニング研究において、Transformerの自己注意機構が持つ \(O(L^2)\) （系列長 \(L\) の2乗）の計算量を回避する手段として再評価されています。

Gu らのS4（Structured State Space Sequence model）（2022, ICLR）は、状態空間モデルから導かれる大域的な畳み込みカーネルをFFTで適用することで、系列全体にわたる依存関係を \(O(L \log L)\) で捉えます。これを継承したHyena Hierarchy（Poli et al., ICML 2023）は、自己注意を「陰的にパラメータ化された長い畳み込み」と乗法的ゲーティングの組み合わせで完全に置き換えるアーキテクチャで、最適化済みの注意機構に対し系列長8,192で同等、系列長65,536では約100倍高速という結果を報告しています。数万〜数十万トークンの推論タスクでは既存の劣線形手法（状態空間モデル等）に対して50ポイント以上の精度改善を達成し、標準的な言語モデリングベンチマーク（WikiText103, The Pile）でも密な注意機構を使わないアーキテクチャとして最高性能を示しました。

もっとも、式\((10)\) のFFT畳み込みは3回のFFT（順方向2回・逆1回）と要素積という複数のメモリ律速なカーネル呼び出しを要するため、GPU上では理論上の計算量優位性がそのまま実行時間に反映されにくいという実装上の課題がありました。FlashFFTConv（Fu et al., 2023）は、FFTをモナーク分解（Monarch decomposition）により行列積の連続として書き換え、GPUのテンソルコアで実行できるようにカーネルを融合する手法で、PyTorch標準のFFT畳み込みに対し最大7.93倍の高速化、エンドツーエンドで最大4.4倍の高速化を報告しています。これにより、Hyena-GPT-sのPerplexityが2.3ポイント改善し、M2-BERT-baseのGLUEスコアが3.3ポイント向上し、高解像度画像タスクPath-512では従来50%を超えられなかった精度が96.1%に到達したと報告されています。

本記事で導出した「畳み込み定理によりFFTで高速化できる」という古典的な事実、そして「線形畳み込みを正しく計算するにはゼロパディングが必要」という注意点は、こうした最新のシーケンスモデリング研究の実装においても本質的に同じ形で登場する基礎理論です。

まとめ

畳み込みは「片方を反転してずらし、積和をとる」操作であり、LTI システムの応答を完全に記述する
畳み込み定理 \(\mathcal{F}\{x \circledast h\} = X \cdot H\) （式\((8)\) ）は、連続時間では変数変換 \(u=t-\tau\) 、離散DFTでは \(\bmod N\) の巡回シフトという変数変換によって、定義から直接証明できる
DFTの畳み込み定理が成り立つのは巡回畳み込みに対してであり、線形畳み込みと一致させるには \(N_{\text{fft}} \ge N+M-1\) を満たすゼロパディングが必須。不十分だと先頭サンプルにエイリアシングが生じる（実験では最大誤差21.0 → ゼロパディング後は \(10^{-15}\) オーダーに収束）
相関は時間反転を伴わない積和で、相互相関 \(R_{xy}\) はテンプレートマッチング、自己相関 \(R_{xx}\) は周期検出に用いる。コーシー・シュワルツの不等式から「既知波形を最大SN比で検出する最適フィルタ＝相関演算（マッチドフィルタ）」であることも厳密に導ける
直接畳み込み \(O(NM)\) とFFT畳み込み \(O((N+M)\log(N+M))\) の実測クロスオーバーは \(N=M\approx 1024\) 付近にあり、\(N=M=65536\) では約295倍の高速化が得られる
信号が極端に長い場合はオーバーラップアド／オーバーラップセーブ法でブロック単位のFFT畳み込みに分割でき、実験では単一の巨大FFTより約4倍高速だった
移動平均はボックス窓との畳み込みであり、FIR フィルタの最も単純なケースとして位置付けられる
Wiener–Khinchin の定理により、自己相関と PSD は同一概念の異なる表現である
2023年以降、畳み込み定理とFFTは Hyena Hierarchy や FlashFFTConv など、Transformerの自己注意に代わる劣二次時間の系列モデリング手法の中核技術として再評価されている

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

高速フーリエ変換（FFT）の仕組みとPython実装 - 本記事の畳み込み定理の前提となる DFT/FFT のアルゴリズムを解説しています。
窓関数とパワースペクトル密度（PSD）の理論とPython実装 - 自己相関と PSD の関係（Wiener–Khinchin 定理）の応用として PSD 推定法を扱っています。
ウェーブレット変換の理論とPython実装 - 畳み込みを基盤とする時間-周波数解析の発展形を解説しています。
移動平均フィルタの種類と比較 - 移動平均がボックス窓との畳み込みであることを踏まえ、SMA・WMA・EMA を比較しています。
自己相関・相互相関の理論とPython実装 - 本記事の相関の定義をベースに、ACF/PACFを使った時系列分析や遅延推定（GCC-PHAT）まで詳しく扱います。
ケプストラム分析の理論とPython実装 - 畳み込みを対数スペクトルで分解する代表的応用。F0推定・スペクトル包絡推定の枠組みを提供します。
FIRフィルタとIIRフィルタの比較 - 畳み込みの定式化を直接的に実装するFIRフィルタと、再帰的構造を取るIIRの設計思想の違いを解説しています。
DTFT・DFT・FFTの関係を体系的に理解する - 畳み込み定理が成立する周波数変換（DTFT・DFT）の階層関係を整理し、FFT畳み込みの理論的背景を補完します。

参考文献

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.
Gu, A., Goel, K., & Ré, C. (2022). Efficiently Modeling Long Sequences with Structured State Spaces. ICLR 2022.
Poli, M., Massaroli, S., Nguyen, E., et al. (2023). Hyena Hierarchy: Towards Larger Convolutional Language Models. ICML 2023. arXiv:2302.10866
Fu, D. Y., Kumbong, H., Nguyen, E., & Ré, C. (2023). FlashFFTConv: Efficient Convolutions for Long Sequences with Tensor Cores. arXiv:2311.05908
NumPy convolve documentation
SciPy signal documentation