円周率.jp (http://xn--w6q13e505b.jp/method/fft/dwt.html)

離散荷重変換

DFT を行う前後で、データ列と同じ長さの適当な係数列 $\{a_i\}$ を付ける変換を考えてみよう。

\[ X_k = \sum_{j=0}^{n-1}a_jx_j\omega^{jk}, \quad x_j = a_j^{-1} n^{-1}\sum_{k=0}^{n-1}X_k\omega^{-jk} \]

この変換を畳み込み乗算のルーチンに組み込むと

\[ \begin{eqnarray} z_l &=& a_l^{-1}n^{-1} \sum_{k=0}^{n-1} \left(\sum_{j=0}^{n-1} a_j x_j \omega^{jk} \right) \left(\sum_{j'=0}^{n-1} a_{j'} y_{j'} \omega^{j'k} \right) \omega^{-kl}\\ &=& \sum_{(j+j') \bmod n = l} (a_l^{-1} a_j a_{j'}) (x_j y_{j'}) \end{eqnarray} \]

という形になり、$\{a_j\}$ の定義次第で計算のバリエーションを増やすことができる。 この計算方法を離散荷重変換 (Discrete Weighted Transform) という。 以下に使い勝手の良い例をいくつか挙げる。

負巡回畳み込み乗算

$\{a_j\} = \{\omega^{j/2}\}$ とすると \[ \begin{eqnarray} z_l &=& \sum_{(j+j') \bmod n = l} \omega^{j+j'-l} (x_j y_{j'})\\ &=& \sum_{j+j'=l} x_j y_{j'} + \omega^{n/2} \sum_{j+j'=l+n} x_j y_{j'}\\ &=& \sum_{j+j'=l} x_j y_{j'} - \sum_{j+j'=l+n} x_j y_{j'} \end{eqnarray} \]

となり負巡回畳み込み乗算になる。

複素巡回畳み込み乗算

計算が複素数体でなされている時、$\{a_j\} = \{\omega^{j/4}\}$ とすると $\omega^{n/4} = -i$ となるので \[ \begin{eqnarray} z_l &=& \sum_{(j+j')\bmod n=l} \omega^{j+j'-l} (x_j y_{j'})\\ &=& \sum_{j+j'=l} x_j y_{j'} + \omega^{n/4} \sum_{j+j'=l+n} x_j y_{j'}\\ &=& \sum_{j+j'=l} x_j y_{j'} - i \sum_{j+j'=l+n} x_j y_{j'} \end{eqnarray} \]

となり、巡回部分が虚数部へ巡回する畳み込み乗算になる。 この DWT は実数 DFT ではないが、 畳み込み乗算の無駄を省く手法の 1 つである。

なお大前提をひっくり返す話ではあるが、FFT の実装の際に DFT の定義式として $\omega=\exp(-2\pi i/n)$ ではなく $\omega=\exp(2\pi i/n)$ を採用すると \[ z_l = \sum_{j+j'=l} x_j y_{j'} + i \sum_{j+j'=l+n} x_j y_{j'} \]

と分かりやすい形になる。

無理数基底離散荷重変換 (IBDWT)

多倍長整数の計算をするとき、通常なら $2^{64}$ や $10^{19}$ の様な定数を基底 $R$ として \[ X = \sum_{j=0}^{n-1} x_jR^j\quad (0\leq x_j \lt R) \]

というように数列 $\{a_j\}$ に分割定義して計算する。

ここではこれを少し変えて $R$ にバリエーションを持たせてみる。 具体的には $r$ を 2 や 10 の小さな単位での基数として、

\[ X = \sum_{j=0}^{n-1} x_j\prod_{k=0}^j R_k\quad (0\leq x_k \lt R_k = r^{d_k}) \]

のように各要素に詰める桁数をバラけさせる。 このとき

\[ D = \sum_{j=0}^{n-1} d_j \]

とすれば多倍長整数 $X$ は $r$ 進数では $D$ 桁の数である。 もしこの計算空間で巡回畳み込み乗算ができれば $\pmod{r^D-1}$ の計算が容易になるが、それを実現する方法がここで紹介する 無理数基底離散荷重変換 (IBDWT ; Irrational Base DWT) である。 ちなみに長々と一般的な話で紹介しているが、 現在、メルセンヌ素数($r=2$、$D\in P$)の確認以外に使いみちは知られてない。

計算方法

$\{d_j\}$、$\{a_j\}$ を以下のように定義する。

\[ d_j = \left \lceil \frac{D(j+1)}{n} \right \rceil - \left \lceil \frac{Dj}{n} \right \rceil \] \[ a_j = r^{\left\lceil \frac{Dj}{n} \right\rceil - \frac{Dj}{n}} \Leftrightarrow \log_r a_j = \left\lceil \frac{Dj}{n} \right\rceil - \frac{Dj}{n} \]

わかりにくいが、 $a_j$ は $\dfrac{Dj}{n}$ の小数部分 $f_j$ を使って $r$ の $(1-f_j)$ 乗を取った数である。厳密には $\dfrac{Dj}{n}$ が整数になる場合は $a_j=r^0=1$ である。 この定義を使って DWT を行うと

\[ z_l = \sum a_{j}a_{j'}a_l^{-1} \cdot x_j y_{j'} = \sum r^{E(j,j',l)} (x_j y_{j'}) \] \[ \left( E(j,j',l) = \left\lceil \dfrac{Dj}{n} \right\rceil + \left\lceil \dfrac{Dj'}{n} \right\rceil - \left\lceil \dfrac{Dl}{n} \right\rceil - \dfrac{D(j+j'-l)}{n} \right) \]

となる。 ちなみに $j+j'-l=\{0,n\}$ ということを考えると $E(j,j',l)=\{0,1\}$ となり、計算結果が整数になることが保証される。