重畳加算法
重畳加算法 (ちょうじょうかさんほう、OverLap-Add method, OA, OLA) とは、非常に長い信号 \mathbf{x} と FIRフィルタ \mathbf{h} の 離散畳み込み \mathbf{h} * \mathbf{x} を分割して(効率的に)処理する手法である。
\begin{align}
y[n] = (\mathbf{h} * \mathbf{x})[n]
=  \sum_{m=0}^{M-1} h[m] \cdot x[n-m]
\end{align}
ここで 信号 x[n]は n=1,\ldots,N_x、フィルタ h[m]は m=0,\ldots,M-1 以外で0、またM<N_xである。
重畳加算法の基本アイデアは、長い信号 \mathbf{x} を区間長Lで区切り、複数の短い断片 \mathbf{x}_k と フィルタ \mathbf{h} に関する複数の畳み込みに分割して、適切なブロック長のFFTの積として効率的に計算することにある。
\begin{align}
x_k[n]\ & \stackrel{\mathrm{def}{=} \
\begin{cases}
x[n+kL],& n=1,2,\ldots,L\\
0, & \textrm{otherwise}
\end{cases}\\
x[n] & = \sum_{k} x_k[n-kL], \\
\end{align}
離散畳み込み(\mathbf{h} * \mathbf{x})[n]は、各区間の畳み込みの総和で表せる:
\begin{align}
y[n]   = (\mathbf{h} * \mathbf{x})[n]
& = \sum_{m=0}^{M-1} \left(h[m] \cdot \sum_{k} x_k[n-m-kL]\right) \\
& = \sum_{k} (\mathbf{h} * \mathbf{x}_k)[n - kL] \\
\end{align}
ここで各区間の畳み込み y_k[n] \ \stackrel{\mathrm{def}{=} \ (\mathbf{h} * \mathbf{x}_k)[n] は 領域 [1, L+M-1] 以外で 0 であり、任意の N \ge (L+M-1) に関し、領域 [1, N] における \mathbf{x}_k と \mathbf{h} の N点巡回畳み込みと等価である。
重畳加算法のアドバンテージは、この巡回畳み込みが畳み込み定理により、次のような FFTの積の逆FFT として効率的に計算できる事にある:
ここで \textrm{FFT}と\textrm{IFFT}は(ブロック長Nの)高速フーリエ変換と逆高速フーリエ変換である。
FFTアルゴリズムによっては、(巡回畳み込み計算のために)FFTブロック長Nを調整する事が理に適っている。例えばCooley-Tukey型FFTアルゴリズム (radix-2アルゴリズム) を使う場合、2の冪乗のブロック長を選ぶと有利である:
N = L + M - 1 = 2^p, \quad p \in \mathbb{N}
図1に、重畳加算法のアイデアを示す。
信号\mathbf{x} (長さNx) を区間長Lで区切り、オーバーラップの無いk個の断片の列\mathbf{x}_kを得る。
各区間について、断片 \mathbf{x}_k (長さ≦L) と フィルタ \mathbf{h} (長さM) のFFT結果の積をとり逆FFTして、区間毎の畳み込み結果 \mathbf{y}_k (長さ≦L+M-1) を得る。 (なお畳み込み結果は、元信号(区間長L)と比べ (フィルタの長さ-1) だけ長くなり、オーバーラップが生じる)
最終的な出力信号は、図に示すように、各区間の結果\mathbf{y}_k のオーバーラップ(重畳)部を加算しながら連結して得られる。
区間長Lは、FFT開発初期にはしばしば効率上の理由でFFTのブロック長Nが2の冪乗をとるよう調整されたが、更なる研究開発の結果、Nを大きな素数で素因数分解する効率的変換方法が明らかにされ、このパラメータに対する計算上の敏感さは低減された。2013年1月 (\Pi{n_i} と素因数分解できる場合の計算量はO(N\sum{n_i})であり、ブロック長Nをゼロ・パディングで調整して計算量を最適化できる。)
アルゴリズムの疑似コードは以下の通り:
アルゴリズム 1 (重畳加算法(OLA)による線形畳み込み)
使用するFFTアルゴリズムに応じてFFTブロック長 N と 分割区間長 L に最適値を設定
H = FFT(h,N)       (ゼロ・パディングしたFFT )
i = 1
while i <= Nx
il = min(i+L-1,Nx)
yt = IFFT( FFT(x(i:il),N) * H, N)
k  = min(i+M-1,Nx)
y(i:k) = y(i:k) + yt    (オーバーラップ区間を加算 )
i = i+L
end
一般に信号\mathbf{x}が周期的でその周期がN_xの場合、畳み込み結果y[n]も周期的で同じ周期をとる。
1周期分のy[n]は、ちょうど1周期分の信号 x[n] (長さNx) と フィルタh[n] (長さM) の畳み込みで得られる。ここでは前述のアルゴリズム 1を使う。 畳み込み結果y[n]は、区間 [M, Nx] で正しい。
先頭のM-1個(区間[1, M-1])の値に、末尾のM-1個(区間[Nx+1, Nx+M-1])の値を加算すれば、区間 [1, Nx] が正しい畳み込み結果を表すようになる。
変更した疑似コードは以下の通り:2013年1月
アルゴリズム 2 (重畳加算法(OLA)による巡回畳み込み)
アルゴリズム 1 を評価
y(1:M-1) = y(1:M-1) + y(Nx+1:Nx+M-1)
y = y(1:Nx)
end
【参考】英語版記事初版のアルゴリズム 2
(アルゴリズム 1 の必要範囲を引用 )
使用するFFTアルゴリズムに応じてFFTブロック長 N と 分割区間長 L に最適値を設定
H = FFT(h,N)       (ゼロ・パディングしたFFT )
(ここまで引用 )
ML = floor((N-1)/L)    (未使用 )
i = Nx-L+1
k = N - L
while k >= 1
il = i+L-1
yt = IFFT( FFT(x(i:il,N)) * H )
y(1:k) = y(1:k) + yt(N-k+1:N)
k  = k-L
i  = i-L
訳注: 英語版記事 Overlap–add method 04:36, 10 December 2012 (UTC) 版のこの章にはいくつか問題があります:
計算量C_{OA}, C_sが根拠不明 — 特にC_sの表式は、radix-2アルゴリズムの計算量に過ぎない.
計算量の比較意図が不明 — 重畳加算法のブロック畳み込みは巡回畳み込みと同じ形なのでFFTで高速化できるというストーリーなのに、重畳加算法を変形した巡回畳み込みAlgorithm 2と、高速化の鍵であるFFT処理(Eq.1)による(巡回畳み込み)アルゴリズム(詳細不明)を比較しようとしている.
評価基準の説明に不整合があり、"standard ..."の指す具体的アルゴリズムが不明.
本文記述: "standard circular convolution using  Eq.1" — Algorithm 2とは別の、Eq.1に基づくFFTベースの循環畳み込みの標準アルゴリズム
Figure 2キャプション: "time required by  Eq.1" — Eq.1 (="Algorithm 1"相当): 循環畳み込みアルゴリズムではない
図面の解説: "standard FFT method" — 標準FFT: 循環畳み込みアルゴリズムではない
本文中の計算量C_s: N_xが2の冪乗の場合の radix-2アルゴリズムの形 — 重畳加算法より計算量が小さいアルゴリズムが"standard"と呼ばれている
評価対象"Algorithm 2"の計測と図面作成後に、アルゴリズムが書き換えられている.
記事初版のFigure 2キャプションによれば、20:11, 9 June 2007 版Algorithm 2の計算時間を実測(→図面に反映)
図面は02-Jun-2007 20:20:47作成
その後、12:42, 12 July 2008 (UTC)にAlgorithm 2が書き換えられている
畳み込みの計算コストは、操作に関わる複素数乗算の回数と関連付ける事ができる。主要な計算量はFFT演算によるもので、Radix-2アルゴリズム(ブロック長Nが2の冪乗のアルゴリズム)を長さN_x=Nの信号に適用する場合およそ C = (N\log_2{N})/2回の複素数乗算が行われる。重畳加算法における複素数乗算の回数は、FFTとフィルタ乗算とIFFTを考慮して:
C_{OA}=\left\lceil \frac{N_x}{N-M+1}\right\rceil
N\left(\log_2 N+1\right)\,2013年1月
なお巡回行列版のM_Lセクション (先頭と末尾の長さ(M-1)の区間?) の追加コストは、通常とても小さいので単純化のために無視できる。
FFTのブロック長Nの最適値は、\log_2{M}\le n\le \log_2{N_x}の範囲で動かしてC_{OA}\left(N=2^n\right)の最小値を整数nを(数値的に)探す事で得られる。Nが2の冪乗であれば、FFTを効率的に計算できる。Nの値が定まれば、信号\mathbf{x}を最適に区切る区間長L=N-M+1が定まる。
比較のため、普通の巡回畳み込みの計算量も示しておく:
C_S=N_x\left(\log_2 N_x+1\right)\,2013年1月
したがって重畳加算法の計算量はほぼ O(N_x\log_2{N})でスケールし、他方、普通の巡回畳み込みの計算量はほぼ O(N_x\log_2{N_x})である。(N>N_xなので重畳加算法の方が計算量が多い事になってしまう。2つの計算量は要検証) しかしながら、この種の推計は複素数乗算の計算量だけ考慮しており、アルゴリズムに関わる他の処理は度外視している。各アルゴリズムに要する計算時間の直接測定こそ、より大きな関心事である。
図2に、式1による標準的な巡回畳み込み2013年1月の計算時間と、アルゴリズム 2の形の重畳加算法による同様な畳み込みの計算時間の比を示す; 縦軸は信号長N_xの対数表示、横軸はフィルタ長N_h=Mの対数表示で、計算時間の比を等高線で示している。両アルゴリズムはMatlabで実装した。青い太線は、重畳加算法の方が標準巡回畳み込みより速い(比>1)領域の境界線である。このケースでは重畳加算法は標準的手法より最大約3倍高速だった。2013年1月