畳み込みに最適な FFT 点数

はじめに

FFT を使って長い入力信号 \(\boldsymbol{x}\) に有限長の固定インパルス応答 \(\boldsymbol{h}\) を畳み込む場合、 計算の主要部分は、図1 の「FFT → → FFT\({}^{-1}\)」の上側経路である (左下の点線で囲まれた FFT は最初に 1 回計算するだけでよいので無視する)。

主要部分の演算量

シンプルに実装した場合の演算量

よく使われる基数 2 の \(N\) 点 FFT の場合、処理は \(\log_2 N\) 段から構成され、各段は \(N/2\) 個のバタフライ演算からなり、 各バタフライ演算には複素乗算 1 回 (4 回の実数乗算) が含まれるので、1回の FFT に必要な実数乗算回数は次のように表せる。 \begin{eqnarray} 4\cdot\frac{N}{2}\log_2 N &=& 2N\log_2 N \label{CostFFT} \end{eqnarray} 逆 FFT の演算量も同じである (正変換または逆変換時に \(1/N\) の乗算が必要だが、事前に \(\boldsymbol{h}\) または \(\boldsymbol{H}\) を \(1/N\) 倍しておけばよいので無視できる)。 周波数特性 \(\boldsymbol{H}\) と入力スペクトル \(\boldsymbol{X}\) の要素ごとの複素乗算 は、 非負周波数のビンに対してのみ行えばよく、1回の複素乗算は4回の実数乗算で構成されるので、この部分の実数乗算回数は次のように表せる。 \begin{eqnarray} 4\cdot\left(\frac{N}{2}+1\right) &=& 2N+4 \label{CostMult} \end{eqnarray} 式(\ref{CostFFT})と式(\ref{CostMult})から、シンプルに実装した場合の 1 回の畳み込みに必要な実数乗算回数は、次のようになる。 \begin{eqnarray} 2\times \left( 2N\log_2 N \right) + (2N+4) &=& 4N\log_2 N + 2N+4 \label{simpleConv} \end{eqnarray}

一般化した演算量モデル

実際の実装では

• DC と Nyquist 周波数の虚部は 0 なので、虚部の乗算を省略する
• 回転子 \(W\) に対して \(W^0=1\) の乗算は省略、\(W^{N/2}=-1\) の乗算は符号反転、\(W^{N/4}=-i,\ W^{3/4N}=i\) の乗算は実部・虚部の入れ替えと符号反転で済ませる
• 2 つの実数波形を一度に FFT する
• 基数 2 以外の FFT を使用する

サンプル当たりの演算量

無限に続く入力信号に対して、入力 1 サンプル当たりの畳み込みに要する演算量を最小化したい。 長さ \(L\) サンプルのインパルス応答を \(N\) 点 FFT を用いて畳み込みを行う場合、 一度に処理する信号の長さ \(M\) は次式を満たさねばならない。 \begin{eqnarray} M &\leq& N-L+1 \end{eqnarray} 線形畳み込みが成立する限界値 \begin{eqnarray} M &=& N-L+1 \end{eqnarray} を採用した場合、入力 1 サンプル当たりの畳み込み演算量 \(c(N)\) は次式で表せる。 \begin{eqnarray} c(N) &=& \frac{C(N)}{M}\ =\ \frac{a N\log N + b N + c}{N-L+1},\ a\neq 0 \label{cN} \end{eqnarray} 【 計算例 】 演算量が式(\ref{CNmodel})で表せ、\(a=4/\log 2,\ b=2,\ c=4\) の場合に、インパルス応答長 \(L=50\) に対する \(c(N)\) を図2に例示する (\(a\) を \(\log 2\) で割っているのは底変換 \(\log_2 \to \log_e\) のため)。

サンプル当たりの演算量を最小にする\(N\)

\(N\) を実変数、\(L,a,b,c\) を定数と考えた場合、\(c(N)\) が最小になるのは \begin{eqnarray} \frac{d}{dN}c(N) &=& \frac{d}{dN} \frac{a N\log N + b N + c}{N-L+1} \\ % &=& \frac{(N-L+1)\displaystyle\frac{d}{dN}(a N\log N + b N + c) - (a N\log N + b N + c)\displaystyle\frac{d}{dN}(N-L+1)}{(N-L+1)^2} \nonumber\\ % &=& \frac{\left\{a(1+\log N) + b\right\}(N-L+1) - (a N\log N + b N + c)}{(N-L+1)^2} \\ &=& \frac{a N +a(1-L) \log N + a (1-L) +b(1-L)-c}{(N-L+1)^2} \ =\ 0\label{dc0} % &=& 0 \end{eqnarray} が成立する時、すなわち式(\ref{dc0}) の分子が \begin{eqnarray} a N +a(1-L) \log N + a (1-L) +b(1-L)-c &=& 0 \end{eqnarray} になる時である。仮定より \(a\neq 0\) だから、両辺を \(a\) で割って定数項を移項し、次のように書ける。 \begin{eqnarray} N +(1-L) \log N &=& (L-1) + \frac{b(L-1)+c}{a} \label{N1L} \end{eqnarray} 式(\ref{N1L})の左辺は \begin{eqnarray} N +(1-L) \log N &=& (1-L) \frac{N}{1-L}+(1-L) \log N \\ &=& (1-L) \log e^{\frac{N}{1-L}}+(1-L) \log N \\ % &=& (1-L) \left(\log e^{\frac{N}{1-L}}+ \log N\right) \\ &=& (1-L) \log N e^{\frac{N}{1-L}} \end{eqnarray} 式(\ref{N1L})の右辺は \begin{eqnarray} (L-1) + \frac{b(L-1)+c}{a} &=& \frac{(a+b)(L-1)+c}{a} \end{eqnarray} と変形できるので、式(\ref{N1L})は次のように書き換えられる。 \begin{eqnarray} (1-L) \log N e^{\frac{N}{1-L}} &=& \frac{(a+b)(L-1)+c}{a} \end{eqnarray} 両辺を \(1-L\) で割って \begin{eqnarray} \log N e^{\frac{N}{1-L}} &=& \frac{(a+b)(L-1)+c}{a(1-L)} \end{eqnarray} 両辺の指数を取ると \begin{eqnarray} N e^{\frac{N}{1-L}} &=& e^{\frac{(a+b)(L-1)+c}{a(1-L)}} \label{Ne} \end{eqnarray} ここで \begin{eqnarray} x &=& \frac{N}{1-L} \label{x} \end{eqnarray} と置くと、式(\ref{Ne})は \begin{eqnarray} (1-L)x e^x &=& e^{\frac{(a+b)(L-1)+c}{a(1-L)}} \end{eqnarray} と書き換えることができ、両辺を \(1-L\) で割って次式を得る。 \begin{eqnarray} x e^x &=& \frac{1}{1-L}e^{\frac{(a+b)(L-1)+c}{a(1-L)}} \end{eqnarray} \(y=x e^x\) の解は Lambert の \(W\) 関数\({}^{[1]}\)を用いて \(x=W(y)\) と表せる。

実験

測定条件

計算例と同じ条件で、長さ \(L=50\) サンプルのインパルス応答を FFT で畳み込むプログラムを作成し、 \(N\) を \(2^6=64\) から \(2^{16}=65536\) の範囲で変えて、入力 1 サンプル当たりの処理時間を計測した。

実験結果

計測した入力 1 サンプル当たりの畳み込み所要時間を図4に示す。 \(c(N)\) は演算量(今回は乗算回数)の推定値であって処理時間の推定値ではないが、当てはまりの度合いを見るため、スケーリングした \(c(N)\) も実線で重ね描きした。

考察

\(N=2^{11}\)辺りからサンプル当たりの演算量 \(c(N)\) と実際の計算時間が比例しなくなるのは、\(N\) が大きいとデータ全体が CPU のキャッシュに入り切らず、遠く離れたアドレス間のバタフライ演算の増加に伴い、キャッシュのヒット率が下がってゆくためと思われる。

まとめ

\(N\) 点 FFT を用い、長さ \(L\) のインパルス応答に対して線形畳み込みが成立する限界サンプル数 \(M=N-L+1\) で入力信号をフレーム分割し、 1回の畳み込みに必要な演算量が \begin{eqnarray} C(N) &=& a N\log N + b N + c \nonumber \end{eqnarray} でモデル化される場合、入力 1 サンプル当たりの畳み込み演算量 \begin{eqnarray} c(N) &=& \frac{C(N)}{M}\ =\ \frac{a N\log N + b N + c}{N-L+1},\ a\neq 0 \nonumber \end{eqnarray} が最小になる FFT 点数 \(N_{opt}\) は Lambert の \(W\) 関数を使って次のように表せる。 \begin{eqnarray} N_{opt} &=& (1-L)W_{-1}\left(\frac{1}{1-L}e^{\frac{(a+b)(L-1)+c}{a(1-L)}}\right) \nonumber \end{eqnarray}

参考文献