バックナンバーはこちら。
https://www.simulationroom999.com/blog/stock-predict-matlabpython-backnumber/
はじめに
前回は、フーリエ変換、逆フーリエ変換のMATLABコードを動作させてみた。
まずはFFTと同等の整数倍の周波数特性を取って、
その後、最大周波数を調整すると細かい周波数特性が取れることを確認。
今回は同様の処理をPython(Numpy)でもやれるようコードを作成していく。
登場人物
博識フクロウのフクさん
イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1
エンジニア歴8年の太郎くん
イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1
Python(Numpy)版のコード作成はなんかメンドクサイ
今回は、前回動かしたMATLABのフーリエ変換、逆フーリエ変換のPython(Numpy)版だね。
ぶっちゃけPython(Numpy)版を作るのはメンドクサクなってきたな。
でも、MATLABが常に使えるとは限らないし、
大体似たようなコードになるんでしょ？
まぁそうだね。
似たようなコードになるというより、似たコードに寄せてるってイメージだけど。
とりあえず、見てみよう。
フーリエ変換、逆フーリエ変換のPython(Numpy)版のコード
コードはこれになる。
import numpy as np
import matplotlib.pyplot as plt
L=np.pi # 波形の期間(-L～L)
w_max = 20 # 取りたい最大周波数
x=np.linspace(-L,L,255)
ft=np.sin(x)+np.sin(3*x)+np.sin(7*x) # 変換用波形読み込み
N=len(ft) # 波形のplot数取得
dt=2*L/N # dt
dw=w_max/(N/2) # dω 取りたい最大周波数から逆算
# 無限長、無限次元の関数同士の内積を疑似的に実現するため、
# 時間領域関数と周波数領域関数は同じ要素にするに必要あり。
w=np.linspace(-dw*N/2,dw*N/2,N) # ω_n
t=np.linspace(-L,L,N) # t_n
Fw = np.zeros(len(ft), dtype=complex) # F(ω) フーリエ変換後の関数格納用
# F(ω)=∫f(t)e^(-iωt)dt
cnt=0
for tn in t:
Fw[cnt]=ft@np.exp(-1j*w*tn).T*dt
cnt = cnt+1
fx = np.zeros(len(Fw), dtype=complex) # f(x) 逆フーリエ変換後の関数格納用
# f(x)=(1/2π)∫F(ω)e^(iωt)dω
cnt=0
x=t
for wn in w:
fx[cnt]=Fw@np.exp(1j*wn*x).T*dw/(2*np.pi)
cnt = cnt+1
fig = plt.figure()
ax = fig.add_subplot(3, 1, 1)
ax.plot(w,np.abs(Fw),'b',lw=3)
ax.set_xlim(w[0],w[-1])
ax.grid(linestyle='dotted')
c=int(len(w)/2)
ax = fig.add_subplot(3, 1, 2)
mask = (0<=w) & (w<=12);
ax.plot(w[mask],np.abs(Fw[mask]),'b',lw=3)
ax.grid(linestyle='dotted')
ax = fig.add_subplot(3, 1, 3)
ax.plot(t,ft,'b','b',lw=5)
ax.plot(t,fx.real,'r',lw=2)
ax.set_xlim(t[0],t[-1])
ax.grid(linestyle='dotted')
plt.show()コード見た感想
やっぱりMATLABの流れとかわらないね。
(変わらないようにしたんだよ・・・。)
演算の中の「@」って何だっけ？
「@」は内積だな。
普通の積は「*」なのだが、これをベクトルに対して行うとアダマール積になる。
アダマール積？
こんな感じだな。
>>> import numpy as np
>>> A=np.array([1,2,3])
>>> B=np.array([4,5,6])
>>> A@B
32
>>> A*B
array([ 4, 10, 18])
ほう！
アダマール積の方は各ベクトルの要素単位の掛け算になるのか！
というわけで意味が全く異なるので注意が必要だな。
まとめ
まとめだよ。
- フーリエ変換、逆フーリエ変換のPython(Numpy)版のコードを作成。
- 基本的にはMATLABと一緒。
- というか、MATLABに寄せた。
- 内積の演算子は「＠」。
- 「*」だとアダマール積になり、結果が全く異なる。
バックナンバーはこちら。
コメント