重畳加算法

重畳加算法 (ちょうじょうかさんほう、英語: 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,Nx$、フィルタ $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 として効率的に計算できる事にある:

$y_k[n] = \textrm{IFFT}\left(\textrm{FFT}\left(\mathbf{x}_k\right)\cdot\textrm{FFT}\left(\mathbf{h}\right)\right)[n]$

(式1)

ここで $\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: 重畳加算法

図1に、重畳加算法のアイデアを示す。

1. 信号$\mathbf{x}$ (長さNx) を区間長$L$で区切り、オーバーラップの無いk個の断片の列$\mathbf{x}_k$を得る。
2. 各区間について、断片 $\mathbf{x}_k$ (長さ≦L) と フィルタ $\mathbf{h}$ (長さM) のFFT結果の積をとり逆FFTして、区間毎の畳み込み結果 $\mathbf{y}_k$ (長さ≦L+M-1) を得る。
   (なお畳み込み結果は、元信号(区間長L)と比べ (フィルタの長さ-1) だけ長くなり、オーバーラップが生じる)
3. 最終的な出力信号は、図に示すように、各区間の結果$\mathbf{y}_k$ のオーバーラップ(重畳)部を加算しながら連結して得られる。

区間長$L$は、FFT開発初期にはしばしば効率上の理由でFFTのブロック長$N$が2の冪乗をとるよう調整されたが、更なる研究開発の結果、Nを大きな素数で素因数分解する効率的変換方法が明らかにされ、このパラメータに対する計算上の敏感さは低減された。 (訳注: FFTのブロック長$N$が2の冪乗の場合の計算量は$O(N\log_2{N})$、一般にブロック長が $N=\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. 1周期分の$y[n]$は、ちょうど1周期分の信号 $x[n]$ (長さNx) と フィルタ$h[n]$ (長さM) の畳み込みで得られる。ここでは前述のアルゴリズム 1を使う。
   畳み込み結果$y[n]$は、区間 [M, Nx] で正しい。
2. 先頭のM-1個(区間[1, M-1])の値に、末尾のM-1個(区間[Nx+1, Nx+M-1])の値を加算すれば、区間 [1, Nx] が正しい畳み込み結果を表すようになる。

変更した疑似コードは以下の通り:

   アルゴリズム 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
   end

計算量

畳み込みの計算コストは、操作に関わる複素数乗算の回数と関連付ける事ができる。主要な計算量は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)\,$

なお巡回行列版の$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)\,$

したがって重畳加算法の計算量はほぼ $O(N_x\log_2{N})$でスケールし、他方、普通の巡回畳み込みの計算量はほぼ $O(N_x\log_2{N_x})$である。(訳注: $N>N_x$なので重畳加算法の方が計算量が多い事になってしまう。2つの計算量は要検証) しかしながら、この種の推計は複素数乗算の計算量だけ考慮しており、アルゴリズムに関わる他の処理は度外視している。各アルゴリズムに要する計算時間の直接測定こそ、より大きな関心事である。

図2に、式1による標準的な巡回畳み込みの計算時間と、アルゴリズム 2の形の重畳加算法による同様な畳み込みの計算時間の比を示す; 縦軸は信号長$N_x$の対数表示、横軸はフィルタ長$N_h=M$の対数表示で、計算時間の比を等高線で示している。両アルゴリズムはMatlabで実装した。青い太線は、重畳加算法の方が標準巡回畳み込みより速い(比>1)領域の境界線である。このケースでは重畳加算法は標準的手法より最大約3倍高速だった。

図2: 式1の計算時間と、重畳加算法アルゴリズム 2の計算時間の比; 縦軸は信号長$N_x$の対数表示、横軸はフィルタ長$N_h=M$の対数表示

関連項目

• 重畳保留法 (Overlap-save method)

参考文献

• Rabiner, Lawrence R.; Gold, Bernard (1975). Theory and application of digital signal processing. Englewood Cliffs, N.J.: Prentice-Hall. pp. 63–67. ISBN 0-13-914101-4. 
• Oppenheim, Alan V.; Schafer, Ronald W. (1975). Digital signal processing. Englewood Cliffs, N.J.: Prentice-Hall. ISBN 0-13-214635-5. 
• Hayes, M. Horace (1999). Digital Signal Processing. Schaum's Outline Series. New York: McGraw Hill. ISBN 0-07-027389-8.