離散フーリエ変換

# 概要
音などのサンプルが入ったベクターを、周波数ごとに分解することができる。

例えば、１秒間、40000回、音をサンプリングすると、40000のサンプリング結果が取得できる。
このサンプリング結果を、一つのベクターに入れて、40000要素のベクターを作る。

このベクターを下のような、Standard basisで表されたベクターと考えて、

```math

[1,0,0,0,.....] \\
[0,1,0,0,.....] \\
[0,0,1,0,.....]
```

これを、サイン波のbasisで構成されるDiscrete Fourier basisで表す。
そうすると、各basisの係数を見ることで、音がどの周波数をどれだけの度合いで含んでいるのかがわかる。

ここで、係数をいじって、Standard basisに戻すことで、周波数ごとに音をいじることもできる。

# 詳細
## Discrete Fourier basis

$ \omega = e^{(2\pi/n)i} $ として、各ベーシスを、

$$ F_k(t) = \frac{1}{\sqrt{n}}(\omega^{k})^{t}$$

と定義する。サンプル数が$n$ とすると、$0\leq k \leq n-1$ のベーシスが必要。

### 注意点
- $F_0$は、$ \frac{n}{0} $で周波数が無限大になるので、常に $ \frac{1}{\sqrt{n}}+0i $
- オイラーの公式 $e^{\theta i} = \cos{\theta} + i\cdot \sin{\theta}\;$から、$t=n,2n,3n,...$ごとに、$F_0(t) = 1+0i$ 
- $F_0$の周期は$n$で、周波数は$\frac{1}{n}$になる
- $F_0,F_1,F_2,...$は、$t$がゼロの時には全て$\frac{1}{\sqrt{n}}+0i$を返す。周期はそれそれ違って$\frac{n}{k}$（$k$は$F_k$の$k$）。例えば、$F_2(t)$は、$t=n/2$ で初めて$\frac{1}{\sqrt{n}}+0i$に戻る。
- Normを1にするために、$F_k$に、$\frac{1}{\sqrt{n}}$を掛けている

### 離散フーリエ空間
Discrete Fourier basisで構成されるベクター空間

### サンプルベクター
時間$0, 1, .., n-1$ にサンプリングした$n$個の要素が含まれたベクター

```math

\begin{bmatrix}
S(0) \\
S(1) \\
...  \\
S({n-1}) 
\end{bmatrix}

```

### 離散フーリエMatrix
サンプルベクターに対して、同じ時間を使って、Discrete Fourier basisの列を作り、並べていくと、離散フーリエMatrixができる。

```math

\begin{bmatrix}
F_0(0) & F_1(0) & ... & F_{n-1}(0) \\
F_0(1) & F_1(1) & ... & F_{n-1}(1) \\
... & ... & ... & ... \\
F_0(n-1) & F_1(n-1) & ... & F_{n-1}(0) 
\end{bmatrix}

```

### サンプルベクターの離散フーリエ空間での表現


```math

\begin{bmatrix}
F_0(0) & F_1(0) & ... & F_{n-1}(0) \\
F_0(1) & F_1(1) & ... & F_{n-1}(1) \\
... & ... & ... & ... \\
F_0(n-1) & F_1(n-1) & ... & F_{n-1}(0) 
\end{bmatrix}

\begin{bmatrix}
\phi_0 \\
\phi_1 \\
... \\
\phi_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
S(0) \\
S(1) \\
... \\
S(n-1)
\end{bmatrix}

```

これを解くと、$(\phi_0,\phi_1,...,\phi_{n-1})$が、離散フーリエ空間での、サンプルベクターの表現になる。

### 離散フーリエMatrixのインバース

```math

W(\omega, n)=
\begin{bmatrix}
\omega^{0\cdot 0} & \omega^{0\cdot 1} & ... & \omega^{0\cdot (n-1)} \\
\omega^{0\cdot 0} & \omega^{0\cdot 1} & ... & \omega^{0\cdot (n-1)} \\
... & ... & ... & ... \\
\omega^{(n-1)\cdot 0} & \omega^{(n-1)\cdot 1} & ... & \omega^{(n-1)\cdot (n-1)} 

\end{bmatrix}

```

と定義すると、サンプル数$n$の離散フーリエMatrixは、$F_n = \frac{1}{\sqrt{n}}W(e^{2\pi i}/n, n)$と表せる。このインバースは、$F_n = \frac{1}{\sqrt{n}}W(e^{-2\pi i}/n, n)$になる。

### 高速フーリエ変換
フーリエ変換の行列式をインバースを使って解くと$n^2$の計算量になる。変換の対象になるベクターは通常サイズが大きいため、以下の条件が満たせる場合には、高速フーリエ変換を使うことで、計算量を$O(n \log n)$に減らす。

- $n$が2の累乗
- $\omega^{n}=1$

#### 計算方法

##### 関数FFT

関数$FFT$を以下のように定義する

- 入力: $(\omega, S)$
  - $n$は2の累乗
  - 複素数 $\omega$ 。
  - 複素数のリスト$S=[s_0, s_1, ..., s_{n-1}]$
- 出力: 複素数のリスト$[z_0, z_1, ..., z_{n-1}]$ で、要素は以下を満たす。
$$
z_k = s_0(\omega^k)^0+s_1(\omega^k)^1+s_2(\omega^k)^2+...s_{n-1}(\omega^k)^{n-1}
$$

##### 関数FFTの実装
簡単にするために、

$$
s(x)=s_0+s_1x+s_2{x^2}+...+s_{n-1}x^{n-1}
$$

として、$z_k=s(\omega^k)$とする。

次に、$s$を、添字が偶数のものと奇数のものに分ける。

```math
s_{even} = [s_0, s_2, ..., s_{n-2}] \\
s_{odd} = [s_1, s_3, ..., s_{n-1}] 
```

$s(x)$の形で書くと、

```math
s_{even}(x) = s_0+s_2x+s_4x^2+...+s_{n-2}x^{\frac{n-2}{2}} \\
s_{odd}(x) = s_1+s_3x+s_5x^2+...+s_{n-1}x^{\frac{n-2}{2}}
```

これらを使って

```math
s(x) = (s_0+s_2(x^2)+s_4(x^2)^2+...+s_{n-2}(x^2)^{\frac{n-2}{2}}) + x(s_1+s_3(x^2)+s_5(x^2)^2+...+s_{n-1}(x^2)^{\frac{n-2}{2}}) \\

s(x) = (s_0+s_2x^2+s_4x^4+...+s_{n-2}x^{n-2}) + (s_1x+s_3x^3+s_5x^5+...+s_{n-1}x^{n-1}) \\

s(x) = s_{even}(x^2) + x \cdot s_{odd}(x^2) \\
```

この定義を使って、$s(x)$を、$\omega^0,\omega^1,...,\omega^{n-1},$について書くと、

```math
s(\omega^0) = s_{even}((\omega^0)^2) + \omega^0 \cdot s_{odd}((\omega^0)^2) \\
s(\omega^1) = s_{even}((\omega^1)^2) + \omega^0 \cdot s_{odd}((\omega^1)^2) \\
... \\
s(\omega^{n-1}) = s_{even}((\omega^{n-1})^2) + \omega^{n-1} \cdot s_{odd}((\omega^{n-1})^2)
```

ここで、
$$
(\omega^0)^2,(\omega^1)^2,...,(\omega^{n-1})^2
$$

は、FFTの条件である$\omega^n=1$を使うと、

```math
(\omega^0)^2 = (\omega^0)^2(\omega^{n/2})^2 = (\omega^{0+n/2})^2 \\
(\omega^1)^2 = (\omega^1)^2(\omega^{n/2})^2 = (\omega^{1+n/2})^2 \\
(\omega^2)^2 = (\omega^2)^2(\omega^{n/2})^2 = (\omega^{2+n/2})^2 \\
... \\
(\omega^{\frac{n}{2}-1})^2 = (\omega^{\frac{n}{2}})^2(\omega^{n/2})^2 = (\omega^{\frac{n}{2}-1+\frac{n}{2}})^2 \\
```

となるので、$(\omega^0)^2,(\omega^1)^2,...,(\omega^{n-1})^2$を計算すれば、残り半分にも使える。

##### 関数FFTの実行

```math

f_0 = FFT(\omega^2, [s_0, s_2, s_4, ..., S_{n-2}]) \\
f_1 = FFT(\omega^2, [s_1, s_3, s_5, ..., S_{n-1}])

```

を計算して、計算結果のリスト内の多項式を $s_{even}$と$s_{odd}$ で置き換えると

```math

f_0 = [s_{even}((\omega^2)^0), s_{even}((\omega^2)^1), s_{even}((\omega^2)^2), s_{even}((\omega^2)^{\frac{n}{2}-1}))] \\
f_1 = [s_{odd}((\omega^2)^0), s_{odd}((\omega^2)^1), s_{odd}((\omega^2)^2), s_{odd}((\omega^2)^{\frac{n}{2}-1}))]

```

となるので、これらの要素を使うと、$s(\omega^k)$を以下のように計算できる。


```math
s(\omega^0) = s_{even}((\omega^0)^2) + \omega^0 \cdot s_{odd}((\omega^0)^2) \\
s(\omega^1) = s_{even}((\omega^1)^2) + \omega^0 \cdot s_{odd}((\omega^1)^2) \\
... \\
s(\omega^{n-1}) = s_{even}((\omega^{n-1})^2) + \omega^{n-1} \cdot s_{odd}((\omega^{n-1})^2)
```


## 参考文献
Coding the matrix / Philip N. Klein