Python で 離散フーリエ変換
概要
離散フーリエ変換(以降 DFT: Descrete Fourier Transform) を Python でベタ書きしてみた.
動機
FFTではなく,あくまでもDFTを計算することで,どうやって周波数変換しているのか追いたかったから.
バタフライ演算という発想に至る過程を手で確認したい.
説明
離散フーリエ変換の式は以下
$ F[k]= \sum_{n=0}^{N-1}f[n]\exp^{-j\frac{2\pi nk}{N}}
$ f[n] が$ n番目のサンプリングデータで,最大値$ Nを標本点という.つまりは$ N個の配列である.
$ F[k] が周波数成分で,$ k=0,1,2,\cdots {N-1}となる.こちらも$ N個の配列になる.
例えば$ N=4の時,$ k=3について以下のようになる.
$ F[3]= f[0]\exp^{-j\frac{2\pi\times 0\times 3}{4}} + f[1]\exp^{-j\frac{2\pi\times 1\times 3}{4}} + f[2]\exp^{-j\frac{2\pi\times 2\times 3}{4}} + f[3]\exp^{-j\frac{2\pi\times 3\times 3}{4}}
計算すると以下のようになる
$ F[3]= f[0]\exp^{-j2\pi\frac{0}{4}} + f[1]\exp^{-j2\pi\frac{3}{4}} + f[2]\exp^{-j2\pi\frac{6}{4}} + f[3]\exp^{-j2\pi\frac{9}{4}}
ここでオイラーの公式を思い出す.なんか単位円上の点を回転角で表せたな〜,みたいな.
$ \exp^{j\theta}=\cos{\theta}+j\sin{\theta}
そして,先ほどの式の$ \expの肩に乗っている数に注目する.
$ \theta=2\pi\frac{0}{4},2\pi\frac{3}{4},2\pi\frac{6}{4},2\pi\frac{9}{4}
これを見れば,単位円上を$ 2\pi\frac{3}{4}ずつクルクル回っているイメージが湧くと思う.
そしてこれは$ k=3だったためであり,$ 2\pi\frac{nk}{N}に立ち帰れば,$ kが周波数を示す変数として機能していることが見えてくる.
この$ \expで表される「クルクル」を各サンプリングデータと掛け合わせたものの総和がその周波数の成分となる.
つまるところ,$ k=0\cdots N-1のfor文のなかで$ n=0\cdots N-1の各サンプリングデータと「クルクル」を掛け合わせた総和をfor文で回してあげれば良さそうだとなる.
実装
いよいよコードを書いていく.
方針としては説明した通りである.
計算はnumpyで自然対数,複素数,三角関数の計算に適用した.
結果表示にmatplotlibを使用して,関数の波形表示と,周波数成分表示をプロットした.
以下,コードである.
square wave と sin wave のどちらかをコメントアウトして波形を見比べる.
code:dft.py
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
DFT=np.zeros(N-1,dtype=complex)
n=np.linspace(1,N,num=N-1,endpoint=False)
x=np.linspace(0,N,num=N+1,endpoint=True)
y=0
# square wave
for k in range(1,20):
y += np.sin(2*np.pi*x*(2*k-1)/N)/(2*k-1)
y=(4/np.pi)*y
# sin wave
y = np.sin(2*np.pi*x*2/N)
# DFT
for k in range(0,N-1):
for num in range(0,N-1):
DFTk += ynum*np.exp(-1j*(2*np.pi*(k+1)*num)/N) FREQ=abs(DFT)
# plot
plt.figure(1)
plt.plot(x/N,y)
plt.grid(color='gray')
plt.figure(2)
plt.plot(n,FREQ/N)
plt.xticks(n)
plt.grid(color='gray')
plt.show()
計算結果
N=50
sin波 4Hz
波形
https://gyazo.com/c7e714d40292f8cb17fef00954a4fb16
DFT結果
https://gyazo.com/c073fc51801a62d267e22d2126a1414a
4Hzが検出できていることがわかる.
46Hzが見えているが,そもそもサンプリング定理により50Hzサンプリングでナイキスト周波数25Hz以上は無意味である.所謂エイリアシングである.
グラフ左半分が真のDFT結果である.
矩形波
波形
矩形波のフーリエ級数$ \frac{4}{\pi}\sum_{k=1}^{20}{\frac{\sin{\{2\pi\frac{(2k-1)x}{N}\}}}{(2k-1)}},x=0\cdots Nとした.$ xはサンプリング点.
https://gyazo.com/8deeded9c05d34a9720ef74f2f6be595
DFT結果
$ 2k-1の周波数(つまり奇数 Hz)の周波数が周波数の逆数の比で検出できている.
https://gyazo.com/e3bc648077f93d099a68de7acd7feea5
まとめ
個人的な話になるが,学生時代,理論計算上のフーリエ変換はだいぶ使ったが,離散フーリエ変換については触れた程度であった.
趣味の企てでどうしても必要になってきたため,色々と調べて,やっと実装レベルの理解ができるようになった.
FFTではさらにDFTのexp関数に着目し,これがほとんど定数として扱えることから,積和の回数を減らす工夫をしている.
任意の波形と位相を$ 90^\circずらした正弦波との積で特定周波数成分を抜き出すことができることもわかった.これはSDR方面に活かせるだろう.
実際にグラフにプロットすると理解も深まるし,思い通りに動いた時は感動する.
なんだかんだでまともにPythonを使用した事例になったので,その点でも収穫があった.numpy便利すぎる.
scrapbox の数式表示機能を活用できて,やはり便利だと思った.
本記事に関しましてご連絡はこちらに記載ののTwitterにお願いいたします. (純粋に感想をいただけると喜びます。)
おまけ
この記事を書いて2年経ち、多くの方に閲覧していただいているようです。
どこかの学校で課題にでもなっているのですかね。
ただ、あまり良いコードと解説でないので、参考にはならないと思いますよ。
先に示したコードを書いた際は、本当に全くpythonの文法をよくわかっていませんでした。
defを使用して書き直しました。多少は見通しが良くなったかと思います。
dcostringも便利ですね。
code: dft2.py
import matplotlib.pyplot as plt
import numpy as np
def square_wave(t):
""" 矩形波
t 時間軸
"""
y = 0
for k in range(1, 20):
y += np.sin(2*np.pi*t*(2*k-1)/(t.shape0-1))/(2*k-1) y = (4/np.pi)*y
return y
def sin_wave(t):
""" 正弦波
t 時間軸
"""
y = np.sin(2*np.pi*t*2/(t.shape0-1)) return y
def dft(y):
""" 離散フーリエ変換
y データ
"""
dft = np.zeros(datanum-1, dtype=complex)
for k in range(0, datanum-1):
for num in range(0, datanum-1):
dftk += ynum*np.exp(-1j*(2*np.pi*(k+1)*num)/datanum) f = abs(dft)
return f
def graph_plot(t, y, f):
""" グラフプロット
t 時間軸
y 波形
f DFT結果
figure1 : 波形
figure2 : DFT結果
"""
n = np.linspace(1, f.shape0, num=f.shape0, endpoint=True) plt.figure(1)
plt.plot(t/(t.shape0-1), y) plt.grid(color='gray')
plt.figure(2)
plt.xticks(n)
plt.grid(color='gray')
plt.show()
N = 50 # サンプリング周波数
t = np.linspace(0, N, num=N+1, endpoint=True) # 時間軸 0~N 51点
# --------------
# square wave
# --------------
y = square_wave(t)
# --------------
# sin wave
# --------------
# y = sin_wave(t, N)
# --------------
# dft
# --------------
f = dft(y)
# --------------
# plot
# --------------
graph_plot(t, y, f)