3
応用数学解説
文献あり

離散フーリエ変換

245
0

はじめに

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

離散フーリエ変換(DFT)

定義

複素数のデータ列xiがあり、データ数はnとします。xiのフーリエ変換x^i1n乗根 ω=e2πi/nを使って
x^j=k=0n1xkωjk
と表せます。

逆変換

k=0n1wk=0
また
..w2n+k=wn+k=wk=wn+k=w2n+k=..+
なので
xj=1nk=0n1x^kωjk
となります。jkjkとし離散フーリエ変換したものを要素数で割ることで元に戻ります。

データ列としてa,b,cを使うとω=1+i32として

順変換

A=a+b+cB=a+bω2+cωC=a+bω+cω2

逆変換

1+w+w2=0なので
3a=A+B+C3b=A+Bω+Cω23c=A+Bω2+Cω

畳み込み定理

(fg)i=k=0n1fkgik
とすると
(fg^)j=f^jg^j

f^jg^j=(f0+f1ωj+f2ω2j+ ..+fn1ω(n1)j) ×(g0+g1ωj+g2ω2j+ ..gn1ω(n1)j)=(f0g0+f0g1ωj+f0g2ω2j+ .. +f0gn1ω(n1)jf1gn1+f1g0ωj+f1g1ω2j+ .. +f1gn2ω(n1)jf2gn2+f2gn1ωj+f2g0ω2j+ .. +f2gn3ω(n1)j..)=(fg)0+(fg)1ωj+(fg)2ω2j+ ..+(fg)n1ωn1=(fg^)j

A=a+b+cB=a+bω2+cωC=a+bω+cω2X=x+y+zY=x+yω2+zωZ=x+yω+zω2

AX=(ax+by+cz)+(ay+bz+cx)+(az+bx+cy)BY=(ax+byω+czω2)+(ayω2+bz+cxω)+(azω+bxω2+cy)CZ=(ax+byω2+czω)+(ayω+bz+cxω2)+(azω2+bxω+cy)

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

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

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

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

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

C=0にして高周波成分カットしてみます。
A=a+b+cB=a+bω2+cωC=0
A+B+C=2abωcω2A+Bω+Cω2=2baω2cωA+Bω2+Cω=2caωbω2
複素数になってしまいました。複素共役のペアができていないと、実数のデータには戻りません。BのペアはCなので片方だけを削るとそうなります。なので実軸に対称な形の周波数フィルタにします。
A=a+b+cB=a+12(bω2+cω)C=a+12(bω+cω2)

A+B+C=3a+12(b+c)A+Bω+Cω2=2b+12cA+Bω2+Cω=2c+12b
データが3点だと端の影響がでて分かりにくいですが、移動平均のような形になっています。ローパスフィルタによりデータが"ならされた"わけです。

多項式の求解との関係

複素共役が対になる条件

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

解の形

解の形も離散フーリエ変換と類似しています。例えば子葉さんの記事siyoに先の例と同じ式が現れています。次数nの解ける多項式の解は、ν,a,b,c, ..en乗根より大きい冪根を含まない式として
x1=a+bνn+cν2n+dν3n+ ..+eνn1nx2=a+ωbνn+ω2cν2n+ω3dν3n+ ..+ωn1eνn1nx3=a+ω2bνn+ω4cν2n+ω6dν3n+ ..+ω2(n1)eνn1n..
という形をしておりe282ωの代わりにωνn、つまりνn乗根を使った離散フーリエ変換、もしくはa,bνn,cν2n,..eνn1nの離散フーリエ変換とみなせます。

高速フーリエ変換(FFT)

異なる(j,k)に対してxjwjkが同じになることがあれば、そうなるケースを検知して、同じ計算を防ぐことで、計算量の削減ができます。もちろん、ループを回して計算の重複を検知するなど本末転倒ですから、数学の問題としてこの条件をあらかじめ求めておきたいわけです。これは、あるωjに対して、ωjk1=ωjk2=ωjk3 ..となるようなk1,k2,..の組を見つけるという課題になります。

FFTのキモ

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

Cooley-Turkey法

n=n1n2のときj=j1+j2n1,k=k1n2+k2とおきます
x^j1+j2n1=k2=0n21k1=0n11xk1n2+k2ω(j1+j2n1)(k1n2+k2)=k2=0n21(k1=0n11xk1n2+k2ω(j1k1n2))ω(j1+j2n1)k2
このように変形され、カッコの中が小さいサイズのDFTになります。DFTの計算のωjkを要素とする行列をDFT行列といいます。DFT行列に対するいまの操作を図示すると以下のようになります。
Cooley-Turkey法 Cooley-Turkey法

図と同じことをしてみます。
ω=eπi/3とします。ω3=1,ω4=ω,ω5=ω2に留意して
A=a+b+c+d+e+fB=abω2cωd+eω2+fωC=abω+cω2+deω+fω2D=ab+cd+efE=a+bω2cω+d+eω2fωF=a+bω+cω2deωfω2
並べ替え
A=a+c+e+b+d+fB=acω+eω2bω2d+fωC=a+cω2eωbω+d+fω2D=a+c+ebdfE=acω+eω2+bω2+dfωF=a+cω2eω+bωdfω2
サイズ3のDFTを以下のように求めて
A=a+c+eC=acω+eω2E=a+cω2eωB=b+d+fD=bdω+fω2F=b+dω2fω
A=A+BB=CDω2C=EFωD=ABE=C+Dω2F=E+Fω
このように計算量が削減できます。

Prime-Factor FFT

互いに素な整数n1,n2によりn=n1n2と分解します。
(j,k)を以下のようなルールで(j1,k1),(j2,k2)の組に分解します。
j=j1n2+j2n1modnk=k1n21n2+j2n11n1modn
ここでn11n11n1=1modn2をみたす整数です。
jk=(j1n2+j2n1)(k1n21n2+k2n11n1)modn=j1k1n2+j2k2n1modn
なので
x^j=k=0n1xkωj1k1n2ωj2k2n1
これよりxkを上記のルールで並べ変えることで、二重DFTに変形できます。

並べ替えルール

(j,k)から(j1,k1),(j2,k2)への対応は数式で書くと分かりにくいかもしれませんが、表を使うと分かりやすいです。まずj=j1n2+j2n1modnという対応はn1=3,n2=5のとき以下のように表せます。
!FORMULA[78][2102910366][0] j=j1n2+j2n1modn
またk=j1n21n2+j2n11n1modn という対応は下図のようになります。
!FORMULA[80][-2031344760][0] k=j1n21n2+j2n11n1modn
これは中国の剰余定理に基づく対応です。中国の剰余定理は互いに素なn1,n2よりn=n1n2のときZ/nZZ/n1Z×Z/n2Zであるという定理です。n1=3,n2=5のとき、以下の表のように mod 15 の整数j は mod 3 の整数j1 と mod 5 の整数 j2 の組(j1,j2)と1対1対応します。さきほどの図ではこの(j1,j2)が座標になっていることが分かります。
mod 1501234567891011121314mod 3012012012012012mod 5012340123401234
これらの並び換えを使うとDFT行列は下図のようになります。
Prime Factor FFT Prime Factor FFT
これは二重DFTになっていることが分かります。
Prime Factor FFT Prime Factor FFT

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

Raderのアルゴリズム

nが素数のとき、ωjk1=ωjk2=ωjk3 ..となるk1,k2,..は存在しません。
サイズが素数のDFT行列 サイズが素数のDFT行列
このためこれまでの原理で計算量は削減できません。しかし別の原理によりこの場合も計算量を削減することができます。まずZ/pZCp1という対応に基づいた並べ替えを行います。nの原始根をgとし、j=gp,k=gqとすると、jk=gpqとなります。ωgpq=fpqと考えれば、DFTの計算がxfの畳み込みになります。ただしこれはj=0,k=0を除いた部分になります。図で表すと Raderのアルゴリズム Raderのアルゴリズム
巡回畳み込みは、畳み込み原理によりFFTにより計算できるので、これも計算量の削減になります。

簡単なn=3の場合に結果だけ示すと
A=a+b+cB=a+12((b+c)+(ω2ω)(bc))C=a+12((b+c)(ω2ω)(bc))
a,b+c,bc,ω2ωを使い回すことができています。

私の仮説

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

剰余群のDFT

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

方針

まず0,1,..m1に含まれる数値kωmk=e2πik/mと符号化します。n個の配列(ωms1,ωms2,..ωmsn)のDFTは
x^j=k=0n1ωnjkωmsk
ここでN=mnとします。
ωnjkωmsk=ωNnskmjk
となります。

問題

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

ωNakの形の足し算

!FORMULA[118][2078508914][0]に対する!FORMULA[119][1168296218][0] N=3,4,5,6,7,8に対するωNak
この格子は原点からN角形の頂点に線を引き、同様のことを各頂点を中心として行うということを繰り返すことでできます。これをl回繰り返したとき全てが異なる点になれば、Nl個の点が生まれるのですが、実際には重なる点があるため、それより少なくなります。
n個のZ/mZの元 a,b, ..,zから決まるωma+ωmb+ ..+ωmz が取りうる独立の点を数えると以下のようになります。縦はmで横はnになります。
111111234567361015212849162536495153570126210619376191127
例えばn=2,m=2とすると、11=0,1+1=2,11=2の3点なので、2行2列目は3になっています。これはMathematicaで作ったのですが6行目だけおかしい気がします。コードは以下にあります。 https://www.wolframcloud.com/obj/star-yoshi/Published/AnglePath.nb

計算機のパラドックス

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

問題

n個のZ/mZの元 a,b, ..,zから決まる ωma+ωmb+ ..+ωmz が取りうる独立の点の数をf(n,m)とするとき、f(n,m)はどのような式になるか?

n個のZ/mZの元 a,b, ..,zから決まる ωma+ωmb+ ..+ωmz が取りうる独立の点の集合をSとするとき、a,b, ..,cSにマップする式はどのようになるか?

与えられたjに対しn個のZ/mZの元 a,b, ..,zから
ωNnamj+ωNnb2mj+ωNnc3mj+ ..+ωNnznmj
を計算するうまい方法はあるか?

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

参考文献

投稿日:2024515
更新日:2024518
OptHub AI Competition

この記事を高評価した人

高評価したユーザはいません

この記事に送られたバッジ

バッジはありません。
バッチを贈って投稿者を応援しよう

バッチを贈ると投稿者に現金やAmazonのギフトカードが還元されます。

投稿者

17世紀の数学を学び始めました。 https://www.17centurymaths.com/ このサイト素晴らしい。

コメント

他の人のコメント

コメントはありません。
読み込み中...
読み込み中
  1. はじめに
  2. 離散フーリエ変換(DFT)
  3. 定義
  4. 逆変換
  5. 畳み込み定理
  6. 多項式の求解との関係
  7. 高速フーリエ変換(FFT)
  8. 私の仮説
  9. 剰余群のDFT
  10. 方針
  11. 問題
  12. 問題
  13. 参考文献