Pythonで離散フーリエ変換(DFT) - 元理系院生の新入社員がPythonとJavaで色々頑張るブログ

この記事に書かれている事

Pythonによる離散フーリエ変換(DFT)と逆変換の実装

離散フーリエ変換

離散フーリエ変換は、次の式で信号 $f(x)$ を周波数スペクトル $F(t)$ に変換します。

$F(t) = \sum_{x=1}^{N} f(x)exp(-i\frac{2\pi t x}{N})$
ここで、Nを解析する信号長さ、 $i$ を虚数単位とします。

逆変換は、次の式で周波数スペクトル $F(t)$ から信号 $f(x)$ に変換します。
$f(x) = \sum_{t=1}^{N} F(t)exp(i\frac{2\pi t x}{N})$

Pythonによる実装

__author__ = 'emoson'
import math


def dft(f):
    return [sum([f[x]*math.e**(-1j*2*math.pi*t*x/len(f)) for x in range(len(f))]) for t in range(len(f))]


def idft(F):
    return [sum([(F[t]*math.e**(1j*2*math.pi*t*x/len(F))).real for t in range(len(F))])/len(F) for x in range(len(F))]

この関数は次の様に使います。

#信号
A = [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]

#DFT(周波数スペクトルに変換)
_A = dft(A)

#表示
print("元信号:"+"\t".join([str("%2.3f" % a) for a in A]))
print("周波数スペクトル:\n"+"\n".join([str(a) for a in _A]))
print("-------------------------------------------------------")
print("元信号:"+"\t".join([str("%2.3f" % a) for a in A]))
print("復元した信号:"+"\t".join([str("%2.3f"% a) for a in idft(_A) ]))

このコードを実行すると次の様になります

元信号:1.000	2.000	3.000	4.000	5.000	6.000	7.000	8.000	9.000	0.000
周波数スペクトル:
(45+0j)
(-13.090169943749476+9.510565162951535j)
(-8.090169943749473-2.6286555605956714j)
(-1.9098300562505237-5.877852522924732j)
(3.0901699437494785-4.253254041760197j)
(5+4.898587196589413e-15j)
(3.0901699437494696+4.253254041760204j)
(-1.9098300562505495+5.877852522924751j)
(-8.090169943749476+2.6286555605956554j)
(-13.090169943749475-9.510565162951544j)
-------------------------------------------------------
元信号:1.000	2.000	3.000	4.000	5.000	6.000	7.000	8.000	9.000	0.000
復元した信号:1.000	2.000	3.000	4.000	5.000	6.000	7.000	8.000	9.000	-0.000

ちゃんと元信号に復元されていますね。