はじめに
背景
面白そうなので離散フーリエ変換を行う変換行列をつくりたい.
目標
離散フーリエ変換の変換行列を作成し, 実際に変換とその逆変換ができることを確認する.
離散フーリエ変換
離散フーリエ変換の説明については⇒wikiに丸投げする.
またwikiの数式にそって変換するコードは以前離散フーリエ変換の実装で実装しているのでそちらを.
変換行列
まず離散フーリエ変換の式が以下
${\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}$
これから変換行列を作成する関数をつくる↓
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の例でそれぞれが何なのか表示しておく.
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になっている. ↓ここの部分
M = cmath.e**(-1j * 2 * cmath.pi * T * X / N)
実験
変換元となる波形yが必要になるので, 適当にsinの合成波として作る.
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()
こんな感じ
次に作成したdft_matrix関数にサンプル数Nを渡して変換行列Mを求める.
変換行列Mを用いた行列演算でyを変換する.行列演算はnumpyでnumpy.dotが用意されているのでそれを使う.
# 変換行列
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()
変換結果
次に逆変換を行うが, その際にMの逆行列を求める必要がある.
逆行列への変換もnumpyでnumpy.linalg.invが用意されているのでそれを使う.
後は変換後のfyを先と同じ要領で逆変換する. 結果を元波形と比較してみる.
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
最後に
dftを埋め込んですぐ動くコードを置いておく.(numpyとmatplotlibが無い場合はpipでちょいちょいっといれてください)
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()
まとめ
離散フーリエ変換の変換行列を作成し, 実際に変換とその逆変換ができることを確認した.
結果：面白かったρ(・ω・、)