5
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

Python: 離散フーリエ変換(実装コードとテスト)

Last updated at Posted at 2019-02-22

離散フーリエ変換

参考 ⇒ wiki

離散フーリエ変換の式

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

逆離散フーリエ変換の式

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

実装したdftのコード

dft.py
import cmath

import numpy as np


def dft(f):
    N = len(f)
    A = np.arange(N)
    T = A.reshape(1, -1)
    X = A.reshape(-1, 1)
    # 行列演算形式になっているがwikiの数式そのまま
    M = f * cmath.e**(-1j * 2 * cmath.pi * T * X  / N)
    return np.sum(M, axis=1)


def idft(f):
    N = len(f)
    A = np.arange(N)
    X = A.reshape(1, -1)
    T = A.reshape(-1, 1)
    # 同上
    M = f * cmath.e**(1j * 2 * cmath.pi * X * T  / N)
    return np.sum(M, axis=1) / N


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

テスト

実装コードをimportしてテスト.

%matplotlib inline

import dft
import numpy as np
import cmath
import matplotlib.pyplot as plt

# ステップ
s = 0.1
# 横軸(0~2πまで)
x = np.arange(0, 2 * cmath.pi, s)
y = np.sin(x) + 2 * np.sin(2 * x + cmath.pi)
fy = dft.dft(y)
yd = dft.idft(fy)
print('元波形')
plt.plot(x, y)
plt.show()
print('変換後:実部')
plt.plot(fy.real)
plt.show()
print('変換後:虚部')
plt.plot(fy.imag)
plt.show()
print('逆変換後')
plt.plot(yd.real)
plt.show()

変換と逆変換ができていることを確認.
実行結果.png

追記

x,tの積部分はNxNパターンあるので, 以下のような行列演算で一気にやることができる.

A = np.arange(4)
A.reshape(1, -1) * A.reshape(-1, 1)

実行結果

array([[0, 0, 0, 0],
       [0, 1, 2, 3],
       [0, 2, 4, 6],
       [0, 3, 6, 9]])

追追記

よく見るとよくある行列用いた空間系の変換みたいにできそうなのでやってみた.
逆変換行列の作成方法が変なのと, 最後Nで割っているところは調整がいる.
⇒できたら今後「変換行列による離散フーリエ変換」みたいな記事書きたい.+.(ノ*・ω・)ノ*.

%matplotlib inline

import dft
import numpy as np
import cmath
import matplotlib.pyplot as plt

# ステップ
s = 0.1
# 横軸(0~2πまで)
x = np.arange(0, 2 * cmath.pi, s)
y = np.sin(x) + 2 * np.sin(2 * x + cmath.pi)
# 変換行列
M = dft_matrix(y)
# 波形 ⇒ M ⇒ 周波数表現
fy = np.dot(M, y)
# 周波数表現 ⇒ -M ⇒ 波形
yd = np.dot(-M, fy) / len(fy)
print('元波形')
plt.plot(x, y)
plt.show()
print('変換後:実部')
plt.plot(fy.real)
plt.show()
print('変換後:虚部')
plt.plot(fy.imag)
plt.show()
print('逆変換後')
plt.plot(yd.real)
plt.show()

実行結果.png

5
6
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
5
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?