高速フーリエ変換概要を理解する

高速フーリエ変換のポイントをまとめてみました.

# 動画

【数分解説|眺めて理解】高速フーリエ変換 FFT: 離散フーリエ変換を回転因子の特性を活かしてコンピュータでの計算を高速にし計算量をNlogNにする手法.【高速フーリエ変換4/4】
<iframe width="560" height="315" src="https://www.youtube.com/embed/GseqjWgmw6Q" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

# 高速フーリエ変換3行概要

- 高速フーリエ変換は離散フーリエ変換をプログラム時に高速に求める計算の工夫
- できることは離散フーリエ変換と変わらないが、早い
- 回転因子の特性を活かして、行列計算の計算を省けるところ省く

# まずは回転因子

回転因子の幾らかの特性を使って行列を簡単にするため、先にそれを把握.

tは整数とする.
回転因子Wは
```math
W_N^{t} = e^{2\pi i \frac{t}{N}}
```
と定義されており、これはtが増えると単位円上に移動していくことを表す.
また$t=N$の時を考えるとこれは一周してゼロに戻ることになる.
離散フーリエ変換で使うケースは上記のWの表現で全て表現できるため、この性質をそのまま活用することができます.

## 周期性 

単位円上にあるため、$t+N$を行うと一周して元の値と同じになります.
```math
W_N^{t+N} = e^{2\pi i \frac{t+N}{N}}
= e^{2\pi i (\frac{t}{N}+1}
= e^{2\pi i (\frac{t}{N})}e^{2\pi i}
= e^{2\pi i (\frac{t}{N})}(1)
```

## 対称性

単位円上にあるため、$N$が偶数で、$N/2$が定義できると
```math
W_N^{t+\frac{N}{2}} = e^{2\pi i \frac{t+\frac{N}{2}}{N}}
= e^{2\pi i (\frac{t}{N}+\frac{\frac{N}{2}}{N})}
= e^{2\pi i (\frac{t}{N}+\frac{1}{2})}
= e^{2\pi i (\frac{t}{N})}e^{2\pi i \frac{1}{2}}
= e^{2\pi i (\frac{t}{N})}e^{\pi i}
= e^{2\pi i (\frac{t}{N})}(-1)
```
となって、tに$N/2$を足すとちょうど反対側に来ることがわかります.
この特性を後で使います

## kの係数の偶奇

![fourier.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/438951b3-1aae-e733-02f9-505535390d08.png)

離散フーリエ変換の係数の行列はこのようになっています.
乗数はkとnに合わせてそれぞれで値が大きくなっています.
当然ですが、kが偶数だと係数も偶数になります.奇数なら奇数になります.
この特性と先ほどの対称性を使ってこの後、式を簡単にします.

# 高速フーリエ変換の計算

N=8の場合で確認.行列は8x8となり、入力データも周波数も8個ずつ得られる.

## 行列の並び替え

この後の処理のために、最初に行の並び替えを行う.
偶数行目を上に、奇数行目を下に移動.

![guki.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/013d65de-eb86-b56a-8296-a98d446e8dbf.png)

上半分を偶数の集まり、下半分を奇数の集まりと捉えて別々に処理します.

小さくそれぞれの回転量を表現してみました.

![small.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/07f7cf4b-397e-2282-238e-5cfb76d360eb.png)

上半分は、4つずれると元の位置に戻ってきています.
これはkが偶数のため、-1回転させたものを偶数回かけるので、結果的に何周かするに過ぎません.
特定の差分に注目すると4x2=8の回転が起きています. 8の倍数なので周期性が現れます.

下半分は、やはり半分回転したままです.
これはkが奇数のため、何周かしますが、必ず最後は半回転することになります.

以上を踏まえると、上側の左半分の要素は一周したものがN/2=4つ右隣に現れることになり、
すなわち同一の行列になることがわかります.
一方で下側の左半分の要素は半周したものがN/2=4つ右隣に現れることになり、
すなわち-1の行列になることがわかります.

![same](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/dbd4fde8-99dd-c79c-a8e3-a4a998087708.png)

## 偶数行列の計算

上側を取り出すと下記のように同じ行列が並びます.
ここで先に二つを足し合わせてから係数にかけるようにすると計算回数を減らすことができます.
実際に計算すると変換まえは64回、変換後は36回の計算になります.

![even.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/8c15c8e3-6867-0ec0-eeb0-28c8bfb6ead0.png)

上記のような工夫で計算回数を減らし高速化しています.

しかしまだ終わりません.
高速フーリエ変換ではさらにこの4x4の行列に対してこれまでと同様の処理を行います.
同様の処理を適用するには、添字が8x8の時と同様にならなくてはなりません.
これまではN=8を暗黙にして、いましたが、これをN=4にすることで以下のようにkの添字を半分にすることができます.

![yakubun.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/0afeb7cf-6bb3-d146-7d80-cb2b7ed36d2e.png)


## 奇数行列の計算

こちらも偶数行列の計算と大体同じです.
今回は-1の係数のため、足し算ではなく、引き算になることに注意してください.

![odd.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/3d3c13a9-1206-df96-81ab-35be5488b861.png)

偶数の時に添字を調整したように、奇数ではもう少し調整が必要です.
4x4にしたときに一行目が全て1になる必要がありますが、これが現状なっていません.
それをあらかじめ引いた値にかけておくことで、偶数側と全く同じ係数にしてしまいます.

![sameeven.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/da549504-1345-486d-ddb3-7f344e2b1cb4.png)

## これでN=8のサイズは完了. あとは繰り返し.

N=8の時の処理はこれで以上です.
同様の処理をN=4,N=2と繰り返すとどのような計算を行えば良いかがわかります.
上記のような理由で、Nは2,4,8,16のように2の累乗であることが求められます.

# バタフライ演算

先ほど出てきた、要素同士の足し算と、引き算の処理がクロスするような処理のためバタフライ演算と呼びます.

![butterfly.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/0af244a0-6f17-f69e-78a0-0ad6bd0c5c7c.png)

# 計算量

高速フーリエ変換の計算量は簡単に求まります.
2のx乗のデータを使って処理してくため、$log_2(2^x)=x$回行列を小さくする処理が行われます.
それぞれの小さくする処理では、N回の足し算とN回の引き算、それと$\frac{N}{2}$回転処理が行われるため、これらを掛け合わせて
全体としてのオーダーは、$O(Nlog_2N)$とわかります

# ビットリバース

ビットリバースは大文字Fの並び順がどうなるかを求める処理になります.

今回の高速フーリエ変換では、偶数番目と奇数番目を並び替える処理を繰り返していました.
初回は、その一番末尾のbitで偶数番めか奇数番目かわかります.
次は、次のbitに注目して並び替えることになります.
これを繰り返していくと1のbitが下に並んでいるはほど下側に集まるようになります.
すなわち、このbitを逆順に並び替えてみると数字順になることがわかります.

![bitreverse.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541716/478eea3c-ce43-be75-d554-3010829951a8.png)

つまり、

1. Fの添字を2進数で表現
2. bitの並び順を逆順にする
3. 求まった値がそのFが出力時にくる場所

となります.