離散フーリエ変換

はじめに

この記事は離散フーリエ変換(DFT)の紹介と、多項式の求解との関連、剰余群に対するDFTについて考えた記事です。

離散フーリエ変換(DFT)

定義

複素数のデータ列$x_i$があり、データ数は$n$とします。$x_i$のフーリエ変換$\hat{x}_i$は$1$の$n$乗根 $\omega = e^{2\pi i/n }$を使って
$$ \hat{x}_j = \sum_{k=0}^{n-1} x_k \omega^{-jk} $$
と表せます。

逆変換

$$ \sum_{k=0}^{n-1} w^k = 0 $$
また
$..w^{-2n+k}=w^{-n+k} = w^k = w^{n+k} = w^{2n+k}=..+$
なので
$$ x_j = \frac{1}{n}\sum_{k=0}^{n-1} \hat{x}_k \omega^{jk} $$
となります。$-jk \to jk$とし離散フーリエ変換したものを要素数で割ることで元に戻ります。

例

データ列として$a,b,c$を使うと$\omega=\frac{1+i\sqrt{3}}{2}$として

順変換

\begin{eqnarray} A &=& a + b + c\\ B &=& a + b \omega^2 + c\omega\\ C &=& a + b \omega + c\omega^2\\ \end{eqnarray}

逆変換

$1+w+w^2=0$なので
\begin{eqnarray} 3a &=& A + B + C\\ 3b &=& A + B \omega + C\omega^2\\ 3c &=& A + B \omega^2 + C\omega\\ \end{eqnarray}

畳み込み定理

$$ (f*g)_i= \sum_{k=0}^{n-1} f_k g_{i-k} $$
とすると
$$ \left(\widehat{f*g}\right)_j = \hat{f}_j\hat{g}_j $$

例

\begin{eqnarray} A &=& a + b + c\\ B &=& a + b \omega^2 + c\omega\\ C &=& a + b \omega + c\omega^2\\ X &=& x + y + z\\ Y &=& x + y \omega^2 + z\omega\\ Z &=& x + y \omega + z\omega^2\\ \end{eqnarray}

\begin{eqnarray} AX &=&(ax+by+cz)+(ay+bz+cx)+(az+bx+cy)\\ BY &=&(ax+by\omega+cz\omega^2)+(ay\omega^2+bz+cx\omega)+(az\omega+bx\omega^2+cy)\\ CZ &=&(ax+by\omega^2+cz\omega)+(ay\omega+bz+cx\omega^2)+(az\omega^2+bx\omega+cy)\\ \end{eqnarray}

\begin{eqnarray} AX+BY+CZ &=& 3(ax+bz+cy)\\ AX+BY\omega+CZ\omega^2 &=& 3(cz+ay+bx)\\ AX+BY\omega^2+CZ\omega &=& 3(by+cx+az) \end{eqnarray}
畳み込みの定義では、$(x,y,z)$,もしくは$(a,b,c)$の進む方向を逆にするので、これであってるのですが、素直な$(ax+by+cz)$という見え方にするには、片方のフーリエ変換を逆位相にするか、$(x,z,y)$という列を使えばよいです。

フーリエ変換の解釈と応用

フーリエ変換はよく、時間のデータを周波数のデータにする変換、また空間のデータを波数のデータにするという説明がされます。そしてフーリエ変換の実用例として、フーリエ変換により周波数のデータにし、これに高周波成分をカットするローパスフィルタを作用させ、ノイズ除去するという例があります。

畳み込みと周波数フィルタ

ノイズ除去の方法として、移動平均を取るという方法があります。それはデータ点をその近傍でとった平均で置き換えるという操作を全ての点について行うことです。そしてこれは畳み込みという演算になります。さてここで、なにか気づかないでしょうか？そうです！！！周波数フィルタを作用させることと、移動平均をとることは、同じことなのです。これが、　畳み込み定理なのです！！！

例

$C=0$にして高周波成分カットしてみます。
\begin{eqnarray} A &=& a + b + c\\ B &=& a + b \omega^2 + c\omega\\ C &=& 0 \end{eqnarray}
\begin{eqnarray} A+B+C &=& 2a-b\omega-c\omega^2\\ A+B\omega+C\omega^2 &=& 2b - a \omega^2 - c \omega \\ A+B\omega^2+C\omega &=& 2c - a \omega - b\omega^2 \end{eqnarray}
複素数になってしまいました。複素共役のペアができていないと、実数のデータには戻りません。$B$のペアは$C$なので片方だけを削るとそうなります。なので実軸に対称な形の周波数フィルタにします。
\begin{eqnarray} A &=& a + b + c\\ B &=& a + \frac{1}{2}(b \omega^2 + c\omega)\\ C &=& a + \frac{1}{2}(b \omega + c\omega^2) \end{eqnarray}

\begin{eqnarray} A + B + C &=& 3a + \frac{1}{2}(b + c) \\ A + B \omega + C \omega^2 &=& 2b+\frac{1}{2}c \\ A + B \omega^2 + C \omega &=& 2c+\frac{1}{2}b \end{eqnarray}
データが３点だと端の影響がでて分かりにくいですが、移動平均のような形になっています。ローパスフィルタによりデータが"ならされた"わけです。

多項式の求解との関係

複素共役が対になる条件

現実には多くの場合$x_i$は実数であり、この場合$\hat{x}_i$の複素共役は必ず対になり、データ数を実質半分にできます。これは、多項式の係数が実数のとき、複素数解は必ず複素共役と対をなすことと似ています。

解の形

解の形も離散フーリエ変換と類似しています。例えば子葉さんの記事siyoに先の例と同じ式が現れています。次数$n$の解ける多項式の解は、$\nu,a,b,c,~..e$ を$n$乗根より大きい冪根を含まない式として
\begin{eqnarray} x_1 &=& a + b \sqrt[n]{\nu}+c\sqrt[n]{\nu^2}+d\sqrt[n]{\nu^3}+~..+e\sqrt[n]{\nu^{n-1}}\\ x_2 &=& a + \omega b \sqrt[n]{\nu}+\omega^2 c\sqrt[n]{\nu^2}+\omega^{3}d\sqrt[n]{\nu^3}+~..+\omega^{n-1}e\sqrt[n]{\nu^{n-1}}\\ x_3 &=& a + \omega^2b \sqrt[n]{\nu}+\omega^4c\sqrt[n]{\nu^2}+\omega^6d\sqrt[n]{\nu^3}+~..+\omega^{2(n-1)}e\sqrt[n]{\nu^{n-1}}\\ .. \end{eqnarray}
という形をしておりe282、$\omega$の代わりに$\omega\sqrt[n]{\nu}$、つまり$\nu$の$n$乗根を使った離散フーリエ変換、もしくは$a,b\sqrt[n]{\nu},c\sqrt[n]{\nu^2},..e\sqrt[n]{\nu^{n-1}}$の離散フーリエ変換とみなせます。

高速フーリエ変換(FFT)

異なる$(j,k)$に対して$x_jw^{-jk}$が同じになることがあれば、そうなるケースを検知して、同じ計算を防ぐことで、計算量の削減ができます。もちろん、ループを回して計算の重複を検知するなど本末転倒ですから、数学の問題としてこの条件をあらかじめ求めておきたいわけです。これは、ある$\omega^j$に対して、$\omega^{jk_1}=\omega^{jk_2}=\omega^{jk_3}~..$となるような$k_1,k_2,..$の組を見つけるという課題になります。

FFTのキモ

アルゴリズムの詳細は置いといて、FFTのキモをずばり言います。それは並べ替えです。$(j,k)$の表を考えると、上に述べた条件を満たす$(j,k_1),(j,k_2),~..$はこの表の中にあるルールで散らばっています。$(j,k)$を一貫した方法で並べ換え、条件をみたす$(j,k)$の組を局在化させるのがFFTのキモです。

Cooley-Turkey法

$n=n_1n_2$のとき$j=j_1+j_2n_1,k=k_1n_2+k_2$とおきます
\begin{eqnarray} \hat{x}_{j_1+j_2n_1} &=& \sum_{k_2=0}^{n_2-1}\sum_{k_1=0}^{n_1-1} x_{k_1n_2+k_2} \omega^{-(j_1+j_2n_1)(k_1n_2+k_2)}\\ &=&\sum_{k_2=0}^{n_2-1}\left(\sum_{k_1=0}^{n_1-1} x_{k_1 n_2 + k_2} \omega^{-(j_1k_1n_2)}\right)\omega^{-(j_1+j_2n_1)k_2} \end{eqnarray}
このように変形され、カッコの中が小さいサイズのDFTになります。DFTの計算の$\omega^{-jk}$を要素とする行列をDFT行列といいます。DFT行列に対するいまの操作を図示すると以下のようになります。
Cooley-Turkey法

例

図と同じことをしてみます。
$\omega=e^{\pi i/3}$とします。$\omega^3=-1,\omega^4=-\omega,\omega^5=-\omega^2$に留意して
\begin{array}{} A &=&a&+&b&+&c&+&d&+&e&+&f\\ B &=&a&-&b\omega^2&-&c\omega&-&d&+&e\omega^2&+&f\omega\\ C &=&a&-&b\omega&+&c\omega^2&+&d&-&e\omega&+&f\omega^2\\ D &=&a&-&b&+&c&-&d&+&e&-&f\\ E &=&a&+&b\omega^2&-&c\omega&+&d&+&e\omega^2&-&f\omega\\ F &=&a&+&b\omega&+&c\omega^2&-&d&-&e\omega&-&f\omega^2\\ \end{array}
並べ替え
\begin{array}{} A &=&a&+&c&+&e&+&b&+&d&+&f\\ B &=&a&-&c\omega&+&e\omega^2&-&b \omega^2&-&d&+&f\omega\\ C &=&a&+&c\omega^2&-&e\omega&-&b \omega&+&d&+&f\omega^2\\ D &=&a&+&c&+&e&-&b&-&d&-&f\\ E &=&a&-&c\omega&+&e\omega^2&+&b \omega^2&+&d&-&f\omega\\ F &=&a&+&c\omega^2&-&e\omega&+&b \omega&-&d&-&f\omega^2\\ \end{array}
サイズ3のDFTを以下のように求めて
\begin{array}{} A' &=&a&+&c&+&e\\ C' &=&a&-&c\omega&+&e\omega^2\\ E' &=&a&+&c\omega^2&-&e\omega\\ B' &=&b&+&d&+&f\\ D' &=&b&-&d\omega&+&f\omega^2\\ F' &=&b&+&d\omega^2&-&f\omega\\ \end{array}
\begin{array}{} A &=&A'&+&B'\\ B &=&C'&-&D'\omega^2\\ C &=&E'&-&F'\omega\\ D &=&A'&-&B'\\ E &=&C'&+&D'\omega^2\\ F &=&E'&+&F'\omega\\ \end{array}
このように計算量が削減できます。

Prime-Factor FFT

互いに素な整数$n_1,n_2$により$n=n_1n_2$と分解します。
$(j,k)$を以下のようなルールで$(j_1,k_1),(j_2,k_2)$の組に分解します。
\begin{eqnarray} j&=&j_1 n_2 + j_2 n_1 \mod n\\ k&=&k_1 n_2^{-1} n_2 + j_2 n_1^{-1}n_1 \mod n \end{eqnarray}
ここで$n_1^{-1}$は$n_1^{-1}n_1 = 1 \mod n_2$をみたす整数です。
\begin{eqnarray} jk &=& (j_1 n_2 + j_2 n_1)(k_1 n_2^{-1}n_2 + k_2 n_1^{-1}n_1) \mod n\\ &=& j_1 k_1 n_2 + j_2 k_2 n_1 \mod n \end{eqnarray}
なので
$$ \hat{x}_j = \sum_{k=0}^{n-1} x_k \omega^{-j_1 k_1 n_2} \omega^{-j_2 k_2 n_1} $$
これより$x_k$を上記のルールで並べ変えることで、二重DFTに変形できます。

並べ替えルール

$(j,k)$から$(j_1,k_1),(j_2,k_2)$への対応は数式で書くと分かりにくいかもしれませんが、表を使うと分かりやすいです。まず$j=j_1 n_2 + j_2 n_1 \mod n$という対応は$n_1=3,n_2=5$のとき以下のように表せます。
$j=j_1 n_2 + j_2 n_1 \mod n$
また$k = j_1 n_2^{-1} n_2 + j_2 n_1^{-1} n_1 \mod n$ という対応は下図のようになります。
$k = j_1 n_2^{-1} n_2 + j_2 n_1^{-1} n_1 \mod n$
これは中国の剰余定理に基づく対応です。中国の剰余定理は互いに素な$n_1,n_2$より$n = n_1n_2 $のとき$\mathbb{Z}/n\mathbb{Z} \cong \mathbb{Z}/n_1\mathbb{Z} \times \mathbb{Z}/n_2\mathbb{Z}$であるという定理です。$n_1=3,n_2=5$のとき、以下の表のように mod 15 の整数$j$ は mod 3 の整数$j_1$ と mod 5 の整数 $j_2$ の組$(j_1,j_2)$と1対1対応します。さきほどの図ではこの$(j_1,j_2)$が座標になっていることが分かります。
\begin{array}{lcccccccccccccccc} \text{mod 15} & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 \\ \text{mod 3} & 0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2 \\ \text{mod 5} & 0 & 1 & 2 & 3 & 4 & 0 & 1 & 2 & 3 & 4 & 0 & 1 & 2 & 3 & 4 \\ \end{array}
これらの並び換えを使うとDFT行列は下図のようになります。
Prime Factor FFT
これは二重DFTになっていることが分かります。
Prime Factor FFT

例

Cooley - Turkey 法の最後の式を並び変えると
$$ \begin{array}{} A &=&A'&+&B'\\ D &=&A'&-&B'\\ E &=&C'&+&D'\omega^2\\ B &=&C'&-&D'\omega^2\\ F &=&E'&+&F'\omega\\ C &=&E'&-&F'\omega\\ \end{array} $$
3つのサイズ2のDFTの計算になることが分かります。Prime-Factor法はCooley Turkey法の出力配列も並び変えすることで、数学的に純粋にしたものという感じに見えます。

Raderのアルゴリズム

$n$が素数のとき、$\omega^{jk_1}=\omega^{jk_2}=\omega^{jk_3}~..$となる$k_1,k_2,..$は存在しません。
サイズが素数のDFT行列
このためこれまでの原理で計算量は削減できません。しかし別の原理によりこの場合も計算量を削減することができます。まず$\mathbb{Z}/p\mathbb{Z} \cong C_{p-1}$という対応に基づいた並べ替えを行います。$n$の原始根を$g$とし、$j=g^p,k=g^{-q}$とすると、$jk=g^{p-q}$となります。$\omega^{g^{p-q}}=f_{p-q}$と考えれば、DFTの計算が$x$と$f$の畳み込みになります。ただしこれは$j=0,k=0$を除いた部分になります。図で表すと Raderのアルゴリズム
巡回畳み込みは、畳み込み原理によりFFTにより計算できるので、これも計算量の削減になります。

例

簡単な$n=3$の場合に結果だけ示すと
\begin{eqnarray} A&=&a+b+c\\ B&=&a+\frac{1}{2}\left(-(b+c)+(\omega^2-\omega)(b-c)\right)\\ C&=&a+\frac{1}{2}\left(-(b+c)-(\omega^2-\omega)(b-c)\right) \end{eqnarray}
$a,b+c,b-c,\omega^2-\omega$を使い回すことができています。

私の仮説

多項式の求解とDFTの類似性、FFTが小さなサイズのDFTに帰着させる方法であることから、多項式の可解性がFFTの計算と関係しているのではないかと思いました。しかしサイズ6のDFTがサイズ3とサイズ2のDFTに分解でき、３次方程式と２次方程式が解けることから、６次方程式も解けることになりそうですが、そうではないので、まだわからないところです。

剰余群のDFT

いままで実数データに対するDFTの話でしたが、例えば文字列など、値の範囲が限られているデータに対してDFTを行う場合どうなるのでしょうか？例えば円周率の数列の連続する10個の数字の平均が大きくなる位置を探すとき、畳み込みの計算にFFTが使えます。このとき0から9まで範囲しかとらない数値を実数として扱ってしまう無駄が生じるはずです。

方針

まず$0,1,..m-1$に含まれる数値$k$を$\omega_m^k=e^{2\pi i k/m} $と符号化します。$n$個の配列$(\omega_m^{s_1},\omega_m^{s_2},..\omega_m^{s_n})$のDFTは
$$ \hat{x}_j=\sum_{k=0}^{n-1} \omega_n^{-jk} \omega_m^{s_k} $$
ここで$N=mn$とします。
$$ \omega_n^{-jk} \omega_m^{s_k} = \omega_{N}^{ns_k-mjk} $$
となります。

問題

剰余群の入力配列のデータの節約は即座にできますが、DFTの出力は$\sum \omega_{N}^{a_k}$の形の足し算になります。例えば$N=3$のときこの和は六角形の格子のどこかなので、この格子を符号化すれば出力のデータも節約できるはずです。

$\omega_{N}^{a_k}$の形の足し算

$N=3,4,5,6,7,8$に対する$\sum \omega_{N}^{a_k}$
この格子は原点から$N$角形の頂点に線を引き、同様のことを各頂点を中心として行うということを繰り返すことでできます。これを$l$回繰り返したとき全てが異なる点になれば、$N^l$個の点が生まれるのですが、実際には重なる点があるため、それより少なくなります。
$n$個の$\mathbb{Z}/m\mathbb{Z}$の元 $a,b,~..,z$から決まる$\omega_m^a+\omega_m^b+~..+\omega_m^z$ が取りうる独立の点を数えると以下のようになります。縦は$m$で横は$n$になります。
\begin{array}{cccccc}1 & 1 & 1 & 1 & 1 & 1 \\2 & 3 & 4 & 5 & 6 & 7 \\3 & 6 & 10 & 15 & 21 & 28 \\4 & 9 & 16 & 25 & 36 & 49 \\5 & 15 & 35 & 70 & 126 & 210 \\6 & 19 & 37 & 61 & 91 & 127 \\\end{array}
例えば$n=2,m=2$とすると、$1-1=0,1+1=2,-1-1=-2$の3点なので、2行2列目は3になっています。これはMathematicaで作ったのですが6行目だけおかしい気がします。コードは以下にあります。 https://www.wolframcloud.com/obj/star-yoshi/Published/AnglePath.nb

計算機のパラドックス

二組の実数で表しているところを格子点の座標として表現できれば、容量の節約が
できそうですが、計算機では場合によって加算有限 > 不可算無限となることがありそう簡単ではありません。というのも実数は適当な粗視化で64 bit にしているのに対し、格子点を整数の配列として表すと配列長によっては莫大なデータ量になる可能性があるからです。なので十分少ない加算有限にできなければ意味がありません。

問題

剰余群に対するDFTの効率化、高速化を考えるにあたって、FFTのときのように
なんとか数学の問題に翻訳できないかという試みですが、いかがだったでしょうか！！！