はじめに

背景

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

目標

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

離散フーリエ変換

離散フーリエ変換の説明については⇒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}$

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

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の例でそれぞれが何なのか表示しておく.

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()

まとめ

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