이산 부립엽 변환

참조²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)

테스트

실장 코드를 기준으로 테스트를 진행하다.

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

변환과 역변환을 확인합니다.

추기

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으로 나누는 부분은 조정이 필요하다.
A:가능하면앞으로는"변환행렬의이산부립엽변환"같은글을 쓰려고합니다.(ノ・ω･)승낙

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