Python: 離散フーリエ変換の実装：変換行列編

# はじめに

## 背景

面白そうなので離散フーリエ変換を行う変換行列をつくりたい.

## 目標

離散フーリエ変換の変換行列を作成し, 実際に変換とその逆変換ができることを確認する.

# 離散フーリエ変換

離散フーリエ変換の説明については⇒[wiki](https://ja.wikipedia.org/wiki/%E9%9B%A2%E6%95%A3%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B)に丸投げする.
またwikiの数式にそって変換するコードは以前[離散フーリエ変換の実装](https://qiita.com/sai-sui/items/9650ff14af54dd147a39)で実装しているのでそちらを.

# 変換行列

まず離散フーリエ変換の式が以下

${\displaystyle F(t)= \sum_{x=0}^{N-1} f(x) e^{-i\frac{2 \pi t x}{N}} \quad \quad }$

tとxをT, Xにおきかえて雰囲気下のような感じ(XもTも0~N-1まで変化する1次元配列)

${\displaystyle F(T)= f(X) e^{-i\frac{2 \pi T X}{N}}=f(X)M}$

これから変換行列を作成する関数をつくる↓

```dft.py
import numpy as np
import cmath

def dft_matrix(N):
    A = np.arange(N)
    T = A.reshape(1, -1)
    X = A.reshape(-1, 1)
    M = cmath.e**(-1j * 2 * cmath.pi * T * X  / N)
    return M
```

`N = 3`の例でそれぞれが何なのか表示しておく.

```python
N = 3
A = np.arange(N)
T = A.reshape(1, -1)
X = A.reshape(-1, 1)
M = cmath.e**(-1j * 2 * cmath.pi * T * X  / N)
print(A)
print('------------------')
print(T * X)
print('------------------')
print(M)
```

実行結果がこんな感じ

```
[0 1 2]
------------------
[[0 0 0]
 [0 1 2]
 [0 2 4]]
------------------
[[ 1. +0.j         1. +0.j         1. +0.j       ]
 [ 1. +0.j        -0.5-0.8660254j -0.5+0.8660254j]
 [ 1. +0.j        -0.5+0.8660254j -0.5-0.8660254j]]
```

`T * X`はN×Nの配列になっていて, その`(i, j)`の要素`(T * X)[i][j]`は　${2^i×2^j}$　となっている.
後は　${M=e^{-i\frac{2 \pi T X}{N}}}$　の残りの部分を計算して`M`になっている.　↓ここの部分

```python
M = cmath.e**(-1j * 2 * cmath.pi * T * X  / N)
```

# 実験

変換元となる波形`y`が必要になるので, 適当にsinの合成波として作る.

```python
import numpy as np
import matplotlib.pyplot as plt

import dft

x = np.linspace(-1, 1, 200)
y  = np.sin(2 * np.pi * x) + np.sin(2 * np.pi * 2.5 * x)

print('元波形')
plt.plot(x, y)
plt.show()
```

こんな感じ

![fx.png](https://qiita-image-store.s3.amazonaws.com/0/365961/fb2fd52d-66c3-5210-a1bc-18393ce6f3d0.png)

次に作成した`dft_matrix`関数にサンプル数`N`を渡して変換行列`M`を求める.
変換行列`M`を用いた行列演算で`y`を変換する.行列演算はnumpyで`numpy.dot`が用意されているのでそれを使う.

```python
# 変換行列
N = len(y)
M = dft.dft_matrix(N)
# 波形 ⇒　M　⇒ 周波数表現
fy = np.dot(y, M)
print('変換後:実部')
plt.plot(fy.real)
plt.show()
print('変換後:虚部')
plt.plot(fy.imag)
plt.show()
```

変換結果

![dftm_fy.png](https://qiita-image-store.s3.amazonaws.com/0/365961/96326c82-4df4-c042-f5fd-fa3cd64a4626.png)

次に逆変換を行うが, その際にMの逆行列を求める必要がある.
逆行列への変換もnumpyで`numpy.linalg.inv`が用意されているのでそれを使う.
後は変換後の`fy`を先と同じ要領で逆変換する. 結果を元波形と比較してみる.

```python
inv_M =  np.linalg.inv(M)
# 周波数表現 ⇒　-M　⇒ 波形
yd = np.dot(fy, inv_M)
print('元波形')
plt.plot(x, y)
plt.show()
print('逆変換後')
plt.plot(yd.real)
plt.show()
```

完璧ですね(・ω・)b

![dftm_fx.png](https://qiita-image-store.s3.amazonaws.com/0/365961/071068fb-4bb4-f927-9234-07967aff6224.png)

# 最後に

dftを埋め込んですぐ動くコードを置いておく.(numpyとmatplotlibが無い場合はpipでちょいちょいっといれてください)

```python
import cmath

import numpy as np
import matplotlib.pyplot as plt


def dft_matrix(N):
    A = np.arange(N)
    T = A.reshape(1, -1)
    X = A.reshape(-1, 1)
    M = cmath.e**(-1j * 2 * cmath.pi * T * X  / N)
    return M

# 元波形生成部
x = np.linspace(-1, 1, 200)
y  = np.sin(2 * np.pi * x) + np.sin(2 * np.pi * 2.5 * x)
print('元波形')
plt.plot(x, y)
plt.show()

# 変換行列の作成と変換
M = dft_matrix(len(y))
fy = np.dot(y, M)
print('変換後:実部')
plt.plot(fy.real)
plt.show()
print('変換後:虚部')
plt.plot(fy.imag)
plt.show()

# 逆変換行列の作成と逆変換
inv_M =  np.linalg.inv(M)
yd = np.dot(fy, inv_M)
print('逆変換後')
plt.plot(yd.real)
plt.show()
```

# まとめ

離散フーリエ変換の変換行列を作成し, 実際に変換とその逆変換ができることを確認した.
結果：面白かったρ(・ω・、)