数百タップのFIRフィルタを、数百万サンプルにわたる音声ストリームやレーダー受信信号に適用したい——こうした場面は信号処理の現場に山ほどあります。素朴に畳み込みの定義式どおり計算すると、1出力サンプルあたりタップ数 $M$ 回の積和が必要で、全体では $O(NM)$ 回の演算がかかります。$N = 10^6$、$M = 1000$ なら $10^9$ 回。これはリアルタイム処理には重すぎます。
ところがFFT（高速フーリエ変換）を使うと、畳み込みは周波数領域の単なる掛け算になり、計算量を劇的に減らせます。しかしここに落とし穴があります。FFTが暗黙に計算するのは「循環畳み込み（円状畳み込み）」であって、私たちが欲しい「線形畳み込み」ではないのです。この食い違いを正しく埋めるテクニックが、オーバーラップ加算法（overlap-add） と オーバーラップ保留法（overlap-save） です。
この手法は身近なところで広く使われています。たとえば、リバーブやイコライザといった音響エフェクトのリアルタイム処理（数千タップのインパルス応答との畳み込み）、ソフトウェア無線（SDR）における広帯域チャネルフィルタリング、そして合成開口レーダー（SAR）の整合フィルタ処理。いずれも「長い信号 × 長いフィルタ」を高速にこなす必要があり、ブロックFFTによる高速畳み込みが心臓部になっています。
本記事の内容
- 線形畳み込みと循環畳み込みの違い、そして時間領域エイリアシングの正体
- ゼロパディングによるエイリアシング回避条件 $L + M – 1$ の導出
- オーバーラップ加算法の仕組み（末尾 $M-1$ サンプルを次ブロックに足す）
- オーバーラップ保留法の仕組み（先頭 $M-1$ サンプルを捨てる）
- 両者の等価性と、ブロック長最適化による1サンプルあたり計算量の最小化
- Pythonで直接畳み込み・overlap-add・overlap-saveを実装し、出力一致と処理時間 $O(N \log N)$ を可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
なぜFFTで畳み込みが速くなるのか
まず出発点を確認します。FIRフィルタによる出力は、入力 $x[n]$ とインパルス応答（フィルタ係数）$h[n]$ の線形畳み込みで与えられます。
$$ y[n] = (x * h)[n] = \sum_{m=0}^{M-1} h[m]\, x[n-m] $$
ここで $h[n]$ は長さ $M$（$M$ タップ）、$x[n]$ は長さ $N$ とします。この式をそのまま計算すると、出力1サンプルにつき $M$ 回の乗算と加算が必要です。出力の長さは $N + M – 1$ なので、総演算量はおよそ $N M$ 回、すなわち $O(NM)$ です。フィルタが長い（$M$ が大きい）ほど、この負担は線形に重くなっていきます。
ここで思い出したいのが、フーリエ変換の 畳み込み定理 です。「時間領域の畳み込みは、周波数領域では単なる積になる」というあの定理です。離散フーリエ変換（DFT）でこれを書くと、ある条件のもとで次が成り立ちます。
$$ \mathrm{DFT}\{x * h\} = \mathrm{DFT}\{x\} \cdot \mathrm{DFT}\{h\} $$
つまり、$x$ と $h$ をそれぞれFFTして周波数領域に移し、各周波数成分を掛け算し、逆FFT（IFFT）で時間領域に戻せば畳み込みが完成する——という筋書きです。FFTの計算量は長さ $L$ の信号に対して $O(L \log L)$。掛け算は $O(L)$ なので、全体も $O(L \log L)$ で済みます。直接畳み込みの $O(NM)$ と比べると、$M$ が大きいほど圧倒的に有利になることが予想できます。
しかし、上の畳み込み定理には「ある条件のもとで」という但し書きを付けました。実は、FFTが計算する畳み込みは私たちが欲しいものとは微妙に違うのです。次のセクションで、この「違い」の正体を突き止めていきます。
線形畳み込みと循環畳み込みのずれ
二つの畳み込み
長さ $L$ の二つの信号 $a[n], b[n]$（$n = 0, \dots, L-1$）に対して、DFTの掛け算とIDFTを通して得られる畳み込みは、循環畳み込み（circular convolution、円状畳み込み） と呼ばれ、次で定義されます。
$$ (a \circledast b)[n] = \sum_{m=0}^{L-1} a[m]\, b[\,(n-m) \bmod L\,], \quad n = 0, 1, \dots, L-1 $$
記号 $\circledast$ は循環畳み込み、$(n-m) \bmod L$ はインデックスを $L$ で割った余り、つまり「$L$ で一周してくるインデックス」を意味します。一方、私たちが本来欲しいのは線形畳み込みです。
$$ (a * b)[n] = \sum_{m} a[m]\, b[n-m] $$
この二つの決定的な違いは、インデックスがはみ出したときの扱いです。線形畳み込みでは $b[n-m]$ のインデックスが範囲外なら単にゼロとみなします（信号は外側でゼロと考える）。ところが循環畳み込みでは、はみ出した分が反対側から「巻き戻って」入り込んできます。$b[-1]$ は $b[L-1]$ に、$b[-2]$ は $b[L-2]$ に化けるのです。
直感的には、線形畳み込みは「信号を一直線のテープと見なす」のに対し、循環畳み込みは「信号を輪っかにつなげて回す」イメージです。輪っかにすると、テープの後ろが前に回り込んで重なり合ってしまう。これがまさにこれから問題にする現象です。
DFTは循環畳み込みを計算している
なぜFFTの掛け算が循環畳み込みになるのか、定義から確かめておきましょう。長さ $L$ のDFTは
$$ A[k] = \sum_{n=0}^{L-1} a[n]\, e^{-j 2\pi k n / L} $$
です。$C[k] = A[k] B[k]$ を逆変換した $c[n]$ を計算します。逆DFTの定義 $c[n] = \frac{1}{L}\sum_{k=0}^{L-1} C[k] e^{j 2\pi k n / L}$ に $C[k] = A[k]B[k]$ を代入すると、
$$ c[n] = \frac{1}{L}\sum_{k=0}^{L-1} \left(\sum_{p=0}^{L-1} a[p] e^{-j2\pi k p/L}\right)\left(\sum_{q=0}^{L-1} b[q] e^{-j2\pi k q/L}\right) e^{j 2\pi k n/L} $$
となります。ここで和の順序を入れ替え、$k$ についての和を先に実行します。$k$ に依存する因子だけを集めると $\sum_{k=0}^{L-1} e^{j 2\pi k (n – p – q)/L}$ という形になり、これは直交性により次の値を取ります。
$$ \sum_{k=0}^{L-1} e^{j 2\pi k (n-p-q)/L} = \begin{cases} L & (n – p – q) \equiv 0 \pmod L \\ 0 & \text{otherwise} \end{cases} $$
つまり $q \equiv n – p \pmod L$ のときだけ $L$ が残り、その他は打ち消し合ってゼロになります。この条件 $q = (n-p) \bmod L$ を $b[q]$ に代入し、係数 $1/L$ と $L$ を相殺させると、
$$ c[n] = \sum_{p=0}^{L-1} a[p]\, b[\,(n-p) \bmod L\,] = (a \circledast b)[n] $$
が得られます。狙いどおり、DFTどうしの積を逆変換すると、循環畳み込みがきっちり出てきました。FFTは単に高速なDFTですから、FFTを使った畳み込みも例外なく循環畳み込みなのです。
ここまでで「FFTの掛け算は循環畳み込みになる」とわかりました。では、循環畳み込みが線形畳み込みからどのようにずれるのか、その「ずれの量」を具体的に見ていきましょう。
時間領域エイリアシングとゼロパディング条件
ずれの正体
線形畳み込みの結果 $y_{\mathrm{lin}}[n] = (a * b)[n]$ は、長さ $L_a + L_b – 1$ にわたって値を持ちます（$L_a, L_b$ はそれぞれの信号長）。ところが循環畳み込み $y_{\mathrm{circ}}[n]$ は長さ $L$ にしか定義されません。$L < L_a + L_b - 1$ の場合、線形畳み込みのうち $n \geq L$ の部分が「行き場を失い」、$\bmod L$ の作用で先頭側（$n = 0, 1, \dots$）に巻き戻って加算されてしまいます。
数式で書くと、循環畳み込みと線形畳み込みの関係は次のようになります。
$$ y_{\mathrm{circ}}[n] = \sum_{r=-\infty}^{\infty} y_{\mathrm{lin}}[n + rL], \quad n = 0, 1, \dots, L-1 $$
これは、線形畳み込みの結果を周期 $L$ で折り返して重ね合わせた形です。$r=0$ の項が本来欲しい値、$r \neq 0$ の項が「余計に混入してくる」成分です。この巻き戻りによる混入を 時間領域エイリアシング（time-domain aliasing） と呼びます。サンプリング定理で周波数領域に現れるエイリアシングの、時間と周波数を入れ替えた双対概念です。
エイリアシングを消す条件
エイリアシングの混入をゼロにするには、$r \neq 0$ の項がそもそも存在しなければよいわけです。$y_{\mathrm{lin}}[n]$ は $n = 0$ から $n = L_a + L_b – 2$ までの $L_a + L_b – 1$ 点でしか非ゼロになりません。したがって、周期 $L$ を十分大きく取って折り返しが起きないようにすればよい。具体的には、
$$ L \geq L_a + L_b – 1 $$
を満たせば、$y_{\mathrm{lin}}[n + rL]$ の $r \neq 0$ の項はすべて非ゼロ領域の外に追いやられ、$y_{\mathrm{circ}}[n] = y_{\mathrm{lin}}[n]$ が $n = 0, \dots, L_a + L_b – 2$ で完全に成立します。
本記事の設定に当てはめましょう。長さ $L$ の入力ブロックと長さ $M$ のフィルタを畳み込む場合、線形畳み込みの長さは $L + M – 1$ です。したがって、循環畳み込みを線形畳み込みに一致させるために必要なFFTサイズ $N_{\mathrm{fft}}$ は
$$ \boxed{\,N_{\mathrm{fft}} \geq L + M – 1\,} $$
です。実際には、入力ブロック（長さ $L$）とフィルタ（長さ $M$）をともに長さ $N_{\mathrm{fft}}$ になるまでゼロパディングしてからFFTします。ゼロを足すことで「巻き戻る先」に空白地帯を作り、エイリアシングが入り込む余地をなくすのです。FFTの効率を考えると $N_{\mathrm{fft}}$ は2のべき乗に切り上げるのが定石です。
ここで根本的な問題が残っています。入力 $x[n]$ の全長 $N$ は数百万サンプルにもなりえます。これをまるごとFFTすると、FFTサイズが $N + M – 1$ となり、巨大なメモリと遅延が発生します。リアルタイム処理では「最初のサンプルが来てから全部揃うまで待つ」など論外です。そこで、入力を短いブロックに分割して処理するという発想が必要になります。これがオーバーラップ法の出発点です。次のセクションで、その分割の仕方を二通り見ていきます。
オーバーラップ加算法（overlap-add）
考え方
オーバーラップ加算法のアイデアは、線形畳み込みの線形性そのものです。入力 $x[n]$ を、重なりのない長さ $L$ のブロック $x_i[n]$ に分割します。
$$ x[n] = \sum_{i=0}^{\infty} x_i[n], \quad x_i[n] = \begin{cases} x[n] & iL \leq n < (i+1)L \\ 0 & \text{otherwise} \end{cases} $$
各ブロック $x_i$ は、$i$ 番目の区間だけ $x$ と同じ値を持ち、それ以外はゼロです。これらを足し合わせれば元の $x$ に戻ります。畳み込みは線形演算なので、和の畳み込みは畳み込みの和に分配できます。
$$ y[n] = x[n] * h[n] = \left(\sum_i x_i[n]\right) * h[n] = \sum_i \big(x_i[n] * h[n]\big) $$
つまり、各ブロックを個別にフィルタリングし、その結果を足し合わせれば全体の出力が得られる、という分解です。
重なりが生まれる理由
各ブロックの畳み込み $y_i[n] = x_i[n] * h[n]$ を考えます。$x_i$ は長さ $L$、$h$ は長さ $M$ なので、$y_i$ の長さは $L + M – 1$ になります。$y_i$ は入力区間 $[iL,\ (i+1)L)$ から始まり、末尾が $M-1$ サンプルだけ次の区間にはみ出します。これは、フィルタが過去 $M-1$ サンプルの記憶を持つため、ブロックの端の効果が境界を越えて尾を引くからです。
このはみ出した末尾 $M-1$ サンプルは、次のブロック $y_{i+1}$ の先頭 $M-1$ サンプルと時間的に重なります。線形性の式 $y = \sum_i y_i$ が要求するのは、この重なり部分を単純に足し合わせることです。これが「オーバーラップ加算（overlap-add）」という名前の由来です。
アルゴリズムの手順
各ブロックの畳み込みを高速にFFTで行うために、FFTサイズを $N_{\mathrm{fft}} \geq L + M – 1$ に取ります（前節の条件）。フィルタ $h$ は最初に一度だけFFTしておけば再利用できます。手順は次のとおりです。
- フィルタ $h[n]$（長さ $M$）を長さ $N_{\mathrm{fft}}$ までゼロパディングし、$H[k] = \mathrm{FFT}(h)$ を一度だけ計算しておく。
- 入力 $x$ を長さ $L$ の重ならないブロック $x_i$ に分割する。
- 各ブロック $x_i$ を長さ $N_{\mathrm{fft}}$ までゼロパディングし、$X_i[k] = \mathrm{FFT}(x_i)$ を計算する。
- 周波数領域で掛け算 $Y_i[k] = X_i[k] \cdot H[k]$ を行う。
- 逆変換 $y_i[n] = \mathrm{IFFT}(Y_i[k])$ で長さ $L + M – 1$ のブロック出力を得る。
- $y_i$ を出力バッファの位置 $iL$ から 加算（足し込み） する。隣接ブロックの重なり $M-1$ サンプルが自然に加算される。
ステップ6の「加算」が肝心です。各 $y_i$ は前ブロックの末尾と後ブロックの先頭で重なるので、その重なりを足すことで境界をまたぐフィルタの効果が正しく再構成されます。ゼロパディングのおかげで各ブロックの循環畳み込みは線形畳み込みと一致しているので、足し合わせの結果は理論上ぴったり $x * h$ になります。
オーバーラップ加算法では「はみ出した尾を足す」ことで境界を処理しました。一方、まったく逆の発想——「余計な部分を捨てる」ことで境界を処理する方法もあります。それがオーバーラップ保留法です。
オーバーラップ保留法（overlap-save）
考え方
オーバーラップ保留法は、ゼロパディングをやめて、循環畳み込みのエイリアシングが起きる部分を最初から捨てるという戦略を取ります。前節で見たように、長さ $N_{\mathrm{fft}}$ の循環畳み込みでは、線形畳み込みの末尾がはみ出して先頭 $M-1$ サンプルに巻き戻り、そこが汚染されます。逆に言えば、先頭 $M-1$ サンプルを除いた残りの部分は、循環畳み込みでも線形畳み込みと一致しているのです。
このことを利用します。各ブロックを処理して循環畳み込みを計算し、汚染された先頭 $M-1$ サンプルを破棄して、残りの「正しい」サンプルだけを採用（save）する。これが「保留法（overlap-save、別名 overlap-discard）」です。
入力ブロックの重ね方
破棄するサンプルがある分、入力ブロックの取り方を工夫する必要があります。FFTサイズを $N_{\mathrm{fft}}$ とすると、各ブロックの出力のうち有効なのは末尾 $N_{\mathrm{fft}} – (M-1)$ サンプルです。連続した出力を隙間なく得るには、入力ブロックを $M-1$ サンプルずつ重ねて(overlap) 取り出します。
具体的には、各入力ブロックの長さを $N_{\mathrm{fft}}$ とし、前のブロックの末尾 $M-1$ サンプルを次のブロックの先頭に再利用します。つまりブロック $i$ は入力の位置 $i \cdot (N_{\mathrm{fft}} – M + 1) – (M-1)$ から $N_{\mathrm{fft}}$ サンプルを切り出します。最初のブロックには過去がないので、先頭に $M-1$ 個のゼロを補います。1ブロックあたり新しく送り出せる有効サンプル数は
$$ L_{\mathrm{step}} = N_{\mathrm{fft}} – (M – 1) $$
です。この量がスループットを決めます。
なぜ末尾は正しいのか
長さ $N_{\mathrm{fft}}$ のブロック $x_i$（過去 $M-1$ サンプルを含む）とフィルタ $h$ の循環畳み込みを考えます。出力サンプル $n$ は
$$ y_{\mathrm{circ}}[n] = \sum_{m=0}^{M-1} h[m]\, x_i[\,(n-m) \bmod N_{\mathrm{fft}}\,] $$
です。$n \geq M-1$ のとき、$n – m$ は $m = 0, \dots, M-1$ の全範囲で $0$ 以上に収まり、$\bmod$ による巻き戻りが起こりません。したがって循環畳み込みと線形畳み込みが一致します。一方、$n < M - 1$ では $n - m < 0$ となる項があり、その項が $x_i$ の末尾から巻き戻った値（前ブロックとは無関係な値）を拾うため、汚染されます。だからこそ先頭 $M-1$ サンプルを捨て、$n = M-1, \dots, N_{\mathrm{fft}}-1$ の有効部分だけを残すのです。捨てた分を補うために入力を重ねて切り出す——これがオーバーラップ保留法の整合的な仕組みです。
アルゴリズムの手順
- フィルタ $h[n]$ を長さ $N_{\mathrm{fft}}$ にゼロパディングし、$H[k] = \mathrm{FFT}(h)$ を一度だけ計算する。
- 入力の先頭に $M-1$ 個のゼロを補う。
- 入力から長さ $N_{\mathrm{fft}}$ のブロックを、$M-1$ サンプル重ねながら（ステップ $L_{\mathrm{step}} = N_{\mathrm{fft}} – M + 1$ で）切り出す。
- 各ブロックを $X_i[k] = \mathrm{FFT}(x_i)$、$Y_i[k] = X_i[k] H[k]$、$y_i = \mathrm{IFFT}(Y_i)$ で処理する。
- $y_i$ の 先頭 $M-1$ サンプルを破棄 し、残り $L_{\mathrm{step}}$ サンプルを出力に連結する。
オーバーラップ加算が「足す」、オーバーラップ保留が「捨てる」。一見まったく違う処理ですが、実は両者は同じ線形畳み込みを別の経路で計算しているにすぎません。次に、その等価性と、どちらにも共通する計算量の最適化を見ていきます。
両者の等価性と計算量の最適化
等価性
オーバーラップ加算とオーバーラップ保留は、最終的に同じ $y[n] = x[n] * h[n]$ を出力します。両者の違いは「ブロック境界でフィルタの記憶（過去 $M-1$ サンプルの寄与）をどう扱うか」だけです。
- 加算法: ブロックをゼロパディングして独立に畳み込み、はみ出した尾を足すことで記憶を復元する。
- 保留法: ブロックを重ねて取り、汚染された先頭を捨てることで記憶の不整合を排除する。
どちらも「フィルタが境界をまたいで $M-1$ サンプル分の影響を持つ」という同じ事実に対処しているため、結果は一致します。加算法は加算バッファが要る分わずかに実装が複雑ですが、ゼロパディングのおかげで境界の特別扱いが要りません。保留法は加算が不要で実装が単純な反面、入力の重ね合わせ管理が必要です。実用上は、ストリーミング処理では保留法、ブロック単位の一括処理では加算法が好まれる傾向があります。
1サンプルあたりの計算量
高速畳み込みが本当に速いのか、計算量で確かめましょう。1ブロックの処理には、FFTが2回（入力のFFTとIFFT。フィルタのFFTは前計算で償却）と、周波数領域の複素乗算が $N_{\mathrm{fft}}$ 回必要です。FFTの演算量はおよそ $\frac{N_{\mathrm{fft}}}{2}\log_2 N_{\mathrm{fft}}$ 回の複素積として見積もれます。よって1ブロックあたりの演算量は
$$ C_{\mathrm{block}} \approx 2 \cdot \frac{N_{\mathrm{fft}}}{2}\log_2 N_{\mathrm{fft}} + N_{\mathrm{fft}} = N_{\mathrm{fft}}\big(\log_2 N_{\mathrm{fft}} + 1\big) $$
です。1ブロックで処理できる有効サンプル数は $L = N_{\mathrm{fft}} – M + 1$ なので、出力1サンプルあたりの演算量は
$$ C_{\mathrm{sample}} = \frac{C_{\mathrm{block}}}{L} = \frac{N_{\mathrm{fft}}\big(\log_2 N_{\mathrm{fft}} + 1\big)}{N_{\mathrm{fft}} – M + 1} $$
となります。この式が高速畳み込みの効率を決める中心的な指標です。
最適なFFTサイズ
$C_{\mathrm{sample}}$ を $N_{\mathrm{fft}}$ の関数として眺めると、トレードオフが見えてきます。$N_{\mathrm{fft}}$ を小さくしすぎると分母 $N_{\mathrm{fft}} – M + 1$ が小さくなり（極端には $N_{\mathrm{fft}} = M$ で有効サンプルが1になり）、1サンプルあたりのコストが跳ね上がります。逆に $N_{\mathrm{fft}}$ を大きくしすぎると、分子の $\log_2 N_{\mathrm{fft}}$ が増え、また巨大なFFTが無駄に重くなります。したがって、ある中間の $N_{\mathrm{fft}}$ で $C_{\mathrm{sample}}$ が最小になります。
経験則として、最適なFFTサイズはフィルタ長 $M$ の数倍から十数倍（$N_{\mathrm{fft}} \sim 4M$ から $16M$ 程度）で、2のべき乗に丸めた値が選ばれます。重要なのは、最適化された $C_{\mathrm{sample}}$ は $O(\log N_{\mathrm{fft}}) = O(\log M)$ のオーダーで増えるという点です。直接畳み込みの1サンプルあたり $O(M)$ と比べると、$M$ が大きいほど高速畳み込みの優位が拡大します。たとえば $M = 1000$ なら、直接法が1000回の積に対し、高速畳み込みは数十回程度で済む計算になります。
全体としては、入力長 $N$ に対してブロック数が $N/L$ 個、各ブロックが $O(N_{\mathrm{fft}}\log N_{\mathrm{fft}})$ なので、総演算量は $O(N \log M)$ となり、直接畳み込みの $O(NM)$ から大幅に改善します。一括処理（分割なしで全体を1回FFT）の $O(N\log N)$ と比較しても、ブロック化によって $\log N$ が $\log M$ に置き換わる分、長い入力ではブロック法が有利です。
理論はここまでです。「本当に出力が一致するのか」「本当に速いのか」は、実装して確かめるのが一番です。次のセクションでPythonで三つの方法を実装し、数値で検証しましょう。
Pythonでの実装
直接畳み込みと2つのオーバーラップ法
まず、直接畳み込み（参照用）、オーバーラップ加算法、オーバーラップ保留法の3つを実装します。出力が一致することを確かめるのが目的です。
import numpy as np
def conv_direct(x, h):
"""線形畳み込みの直接計算（参照用）"""
return np.convolve(x, h)
def overlap_add(x, h, L):
"""オーバーラップ加算法
x: 入力信号, h: フィルタ係数(長さM), L: ブロック長"""
M = len(h)
N = len(x)
Nfft = 1
while Nfft < L + M - 1: # FFTサイズを L+M-1 以上の2のべき乗に
Nfft *= 2
H = np.fft.rfft(h, Nfft) # フィルタは一度だけFFT（前計算）
y = np.zeros(N + M - 1) # 出力バッファ（線形畳み込みの長さ）
for start in range(0, N, L):
block = x[start:start + L] # 長さLのブロック
Y = np.fft.rfft(block, Nfft) * H # 周波数領域で積（自動ゼロパディング）
y_block = np.fft.irfft(Y, Nfft) # 逆変換 → 長さ L+M-1 相当
end = start + len(block) + M - 1
y[start:end] += y_block[:end - start] # 重なり部分を加算
return y
ここでのポイントは2つです。np.fft.rfft(block, Nfft) の第2引数 Nfft が自動的にゼロパディングを行うので、明示的にゼロを足す必要がありません。そして出力バッファへの代入が +=（加算）になっているのが加算法の本質で、隣接ブロックの末尾 $M-1$ サンプルがここで足し合わされます。
続いてオーバーラップ保留法です。こちらは「捨てる」処理が中心になります。
import numpy as np
def overlap_save(x, h, Nfft):
"""オーバーラップ保留法
x: 入力信号, h: フィルタ係数(長さM), Nfft: FFTサイズ"""
M = len(h)
N = len(x)
L_step = Nfft - (M - 1) # 1ブロックあたりの有効サンプル数
H = np.fft.rfft(h, Nfft) # フィルタは一度だけFFT
# 先頭に M-1 個のゼロを補い、末尾も切り出しに足りるよう余裕を持たせる
x_pad = np.concatenate([np.zeros(M - 1), x])
n_out = N + M - 1 # 線形畳み込みの出力長
y = np.zeros(n_out)
pos = 0 # 出力書き込み位置
start = 0
while pos < n_out:
block = x_pad[start:start + Nfft]
if len(block) < Nfft: # 末尾はゼロで埋める
block = np.concatenate([block, np.zeros(Nfft - len(block))])
Y = np.fft.rfft(block, Nfft) * H
y_block = np.fft.irfft(Y, Nfft)
valid = y_block[M - 1:] # 先頭 M-1 を破棄
take = min(L_step, n_out - pos) # 出力に必要な分だけ採用
y[pos:pos + take] = valid[:take]
pos += L_step
start += L_step
return y
保留法では valid = y_block[M - 1:] で汚染された先頭 $M-1$ サンプルを切り捨て、start += L_step で入力を $M-1$ サンプル重ねながら進めています。先頭のゼロパディング np.zeros(M - 1) が、最初のブロックの「存在しない過去」を補う役割を果たします。
出力が一致するか検証
3つの方法が同じ結果を返すか、ランダム信号で確かめます。
import numpy as np
# テスト信号とフィルタ
np.random.seed(0)
N = 5000 # 入力長
M = 128 # フィルタ長（タップ数）
x = np.random.randn(N)
h = np.random.randn(M)
# 3つの方法で計算
y_ref = conv_direct(x, h)
y_oa = overlap_add(x, h, L=512)
y_os = overlap_save(x, h, Nfft=1024)
# 最大誤差を表示
err_oa = np.max(np.abs(y_ref - y_oa))
err_os = np.max(np.abs(y_ref - y_os))
print(f"overlap-add 最大誤差: {err_oa:.2e}")
print(f"overlap-save 最大誤差: {err_os:.2e}")
print(f"出力長: ref={len(y_ref)}, oa={len(y_oa)}, os={len(y_os)}")
このコードを実行すると、両手法の最大誤差はおよそ 1e-12 程度（浮動小数点演算の丸め誤差の水準）になります。これは、ゼロパディング条件 $N_{\mathrm{fft}} \geq L + M – 1$ を守る限り、循環畳み込みが線形畳み込みと数値的に完全に一致することを意味します。理論で導いたエイリアシング回避条件が、実装でもそのまま成り立っていることが確認できました。
エイリアシングをわざと起こす
逆に、FFTサイズを条件未満にするとどうなるかを見てみましょう。理論で予言した時間領域エイリアシングが本当に現れるかの確認です。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1)
M = 64
L = 200
h = np.random.randn(M)
x_block = np.random.randn(L) # 1ブロックぶんの入力
y_lin = np.convolve(x_block, h) # 線形畳み込み（正解, 長さ L+M-1）
# 条件を満たすFFTサイズ
Nfft_ok = 1
while Nfft_ok < L + M - 1:
Nfft_ok *= 2
y_ok = np.fft.irfft(np.fft.rfft(x_block, Nfft_ok) * np.fft.rfft(h, Nfft_ok), Nfft_ok)
# 条件を満たさない（小さすぎる）FFTサイズ → エイリアシング発生
Nfft_bad = 256 # L+M-1=263 より小さい
y_bad = np.fft.irfft(np.fft.rfft(x_block, Nfft_bad) * np.fft.rfft(h, Nfft_bad), Nfft_bad)
plt.figure(figsize=(10, 5))
plt.plot(y_lin, 'k-', lw=2, label='Linear conv (true)')
plt.plot(y_ok[:len(y_lin)], 'g--', lw=1.5, label=f'Circular, Nfft={Nfft_ok} (OK)')
plt.plot(y_bad, 'r:', lw=1.8, label=f'Circular, Nfft={Nfft_bad} (aliasing)')
plt.axvline(Nfft_bad - 1, color='red', alpha=0.3, ls='--')
plt.xlabel('Sample index n')
plt.ylabel('y[n]')
plt.title('Time-domain aliasing when Nfft < L+M-1')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('aliasing_demo.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、3つの曲線の関係がはっきり読み取れます。$N_{\mathrm{fft}} = 512$（条件を満たす）の循環畳み込みは線形畳み込みと完全に重なります。一方、$N_{\mathrm{fft}} = 256$（$L+M-1=263$ より小さい）の場合、出力の先頭部分が線形畳み込みからずれています。これは、はみ出した末尾 $263 – 256 = 7$ サンプル分が $\bmod 256$ で先頭に巻き戻り、加算されているためです。理論で導いた $y_{\mathrm{circ}}[n] = \sum_r y_{\mathrm{lin}}[n+rL]$ の折り返しが、まさに目に見える形で現れました。
処理時間のスケーリング計測
最後に、入力長 $N$ を変えながら3手法の処理時間を計測し、計算量のオーダーを可視化します。
import numpy as np
import time
import matplotlib.pyplot as plt
def time_it(func, *args, repeat=3):
"""関数の最小実行時間を計測"""
best = np.inf
for _ in range(repeat):
t0 = time.perf_counter()
func(*args)
best = min(best, time.perf_counter() - t0)
return best
M = 256 # 長いFIRフィルタ（256タップ）
h = np.random.randn(M)
Ns = [2**k for k in range(10, 20)] # N = 1024 ... 524288
t_direct, t_oa, t_os = [], [], []
for N in Ns:
x = np.random.randn(N)
t_direct.append(time_it(conv_direct, x, h))
t_oa.append(time_it(lambda x, h: overlap_add(x, h, L=2048), x, h))
t_os.append(time_it(lambda x, h: overlap_save(x, h, Nfft=4096), x, h))
plt.figure(figsize=(10, 6))
plt.loglog(Ns, t_direct, 'o-', label='Direct conv O(N·M)')
plt.loglog(Ns, t_oa, 's-', label='Overlap-add O(N·log M)')
plt.loglog(Ns, t_os, '^-', label='Overlap-save O(N·log M)')
plt.xlabel('Input length N')
plt.ylabel('Time [s]')
plt.title(f'Convolution timing (M={M} taps)')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('conv_timing.png', dpi=150, bbox_inches='tight')
plt.show()
両対数プロットから、3手法のスケーリングの違いが明瞭に読み取れます。直接畳み込みの曲線が最も急な傾きを持ち、$N$ が大きくなるほど処理時間が一気に増大します（傾きは $O(NM)$ に対応）。これに対し、overlap-add と overlap-save の曲線はほぼ平行で、傾きが緩やかです。両者はほぼ重なっており、前節で論じた等価性と同程度の計算量を裏付けています。$N$ が小さい領域ではFFTのオーバーヘッドにより直接法が勝つこともありますが、$N$ が数万を超えるあたりからオーバーラップ法が逆転し、その差は $N$ とともに開いていきます。$M = 256$ という長めのフィルタでは、大きな $N$ で1桁以上の高速化が得られることが見て取れます。
1サンプルあたり計算量の最適化を確認
最後に、FFTサイズ $N_{\mathrm{fft}}$ を変えたときの1サンプルあたりコストを理論式とともにプロットし、最適点が存在することを確かめます。
import numpy as np
import matplotlib.pyplot as plt
M = 256 # フィルタ長
Nffts = np.array([2**k for k in range(9, 18)]) # 512 ... 131072
# 理論式 C_sample = Nfft(log2 Nfft + 1) / (Nfft - M + 1)
C = Nffts * (np.log2(Nffts) + 1) / (Nffts - (M - 1))
plt.figure(figsize=(9, 5))
plt.semilogx(Nffts, C, 'o-', base=2)
opt = Nffts[np.argmin(C)]
plt.axvline(opt, color='red', ls='--', alpha=0.6, label=f'optimal Nfft={opt}')
plt.xlabel('FFT size Nfft')
plt.ylabel('Cost per output sample (relative)')
plt.title(f'Per-sample cost vs FFT size (M={M})')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('per_sample_cost.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"フィルタ長 M={M} のとき最適FFTサイズ ≈ {opt}（M の約 {opt/M:.1f} 倍）")
このグラフは、前節で導いた $C_{\mathrm{sample}}$ の式がU字型の曲線を描くことを示しています。$N_{\mathrm{fft}}$ が $M$ に近い左端では、有効サンプル数 $N_{\mathrm{fft}}-M+1$ が小さくコストが急上昇します。右に行くと一旦下がって最小値を取り、その後 $\log_2 N_{\mathrm{fft}}$ の増加と巨大FFTの非効率でゆるやかに増加に転じます。最適点はおおむねフィルタ長 $M$ の数倍に位置し、理論で述べた「$N_{\mathrm{fft}} \sim 4M$ から $16M$」という経験則と整合します。実装でFFTサイズを選ぶ際は、この曲線の谷を狙うのが定石です。
まとめ
本記事では、長い信号と長いFIRフィルタの畳み込みを高速化するオーバーラップ加算法・保留法について、理論の導出から実装・検証まで解説しました。
- FFTの畳み込みは循環畳み込み: DFTの直交性から、DFTどうしの積を逆変換すると線形畳み込みではなく循環畳み込み $(a\circledast b)[n] = \sum_m a[m]b[(n-m)\bmod L]$ が得られる。
- 時間領域エイリアシング: 循環畳み込みは線形畳み込みを周期 $L$ で折り返した $y_{\mathrm{circ}}[n]=\sum_r y_{\mathrm{lin}}[n+rL]$ になる。これを防ぐ条件が $N_{\mathrm{fft}} \geq L + M – 1$ であり、ゼロパディングで実現する。
- オーバーラップ加算法: 入力を重ならないブロックに分け、ゼロパディングして個別に畳み込み、はみ出した末尾 $M-1$ サンプルを次ブロックに加算する。線形性 $y=\sum_i(x_i*h)$ がそのまま根拠になる。
- オーバーラップ保留法: 入力を $M-1$ サンプル重ねて切り出し、循環畳み込みで汚染された先頭 $M-1$ サンプルを破棄して残りを採用する。
- 等価性と最適化: 両者は同じ線形畳み込みを別経路で計算する。1サンプルあたりコスト $C_{\mathrm{sample}} = N_{\mathrm{fft}}(\log_2 N_{\mathrm{fft}}+1)/(N_{\mathrm{fft}}-M+1)$ はU字型で、$N_{\mathrm{fft}} \sim 4M$–$16M$ で最小化される。総計算量は $O(N\log M)$ で直接法の $O(NM)$ を圧倒する。
数値実験では、両手法とも直接畳み込みと丸め誤差レベルで一致し、入力長に対する処理時間が直接法よりはるかに緩やかにスケールすることを確認しました。高速畳み込みは、ディジタル信号処理における「長いフィルタを実時間でかける」あらゆる場面の基礎技術です。
次のステップとして、以下の記事も参考にしてください。