バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia5-backnumber/
はじめに
フーリエ係数に至る道。
今回はフーリエ係数を求めるプログラムをPythonで実現。
登場人物
博識フクロウのフクさん
イラスト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
【再掲】フーリエ係数を求めるプログラムフロー
太郎くん
まずは、プログラムフローを再掲。
- csvファイル読み込み
- 各種変数初期化
- フーリエ係数算出
- n=10,50,200のパターンでフーリエ級数で波形を合成
- グラフにプロット
フクさん
今回はフーリエ係数を求めるプログラムをPythonで実現。
Pythonコード
フクさん
まず、使用する波形を取り込んだcsvファイル
フクさん
Pythonコードは以下になる。
import numpy as np
import matplotlib.pyplot as plt
N=1000 # 係数算出項数(同定元波形のplotよりも少なく)
wave=np.loadtxt('wave.csv',delimiter=',') # 同定波形読み込み
points=len(wave) # 波形のplot数取得
fx=np.array(wave) # 波形を行ベクトルへ
dx=2*np.pi/points # 1plotあたりのx軸幅
x=np.linspace(-np.pi,np.pi,points); # -π~+πの範囲で波形plot数分の等差数列
a = np.zeros(N) # a係数群格納用
b = np.zeros(N) # b係数群格納用
for n in range(1,N+1):
# 係数a_n算出
# a_n = (1/π)∫f(fx)cos(nx)dx
a[n-1] = fx@np.cos(n*x).T*dx/np.pi;
# 係数b_n算出
# a_n = (1/π)∫f(fx)cos(nx)dx
b[n-1] = fx@np.sin(n*x).T*dx/np.pi;
# 係数a_0算出
a0=np.sum(fx)*dx/np.pi
Ns = [10,50,200]
fig = plt.figure()
for i in range(0,len(Ns)):
NN = Ns[i]; # 今回のa_n,b_n項数
# f(x)=a_0+Σ(a_n cos(nx)+ b_n sin(nx))
Fourier_series=np.ones(points)*a0/2
for n in range(1,NN+1):
Fourier_series = Fourier_series+(a[n-1]*np.cos(n*x)+b[n-1]*np.sin(n*x))
ax = fig.add_subplot(len(Ns), 1, i+1)
# 元波形とフーリエ級数波形の表示
ax.plot(x, fx,'-b',lw=3)
ax.plot(x, Fourier_series,'-r',lw=2)
ax.set_title('n={:d}'.format(NN))
ax.set_ylim(-0.1,1.1)
ax.set_xlim(-np.pi,np.pi)
ax.grid(linestyle='dotted')
plt.show()
処理結果
フクさん
処理結果は以下。
考察
太郎くん
結果としてはMATLABと一緒かな。
不連続点で振動するのも一緒だ。
フクさん
まぁ、同じ数式を元にしているからね。
Pythonの0オリジンな振る舞いには注意だな。
まとめ
フクさん
まとめだよ。
- フーリエ係数を求めるプログラムをPythonで実現。
- おおよそ元の波形を再現できる係数が算出できている。
- 不連続点では流石に振動している。
バックナンバーはこちら。
マンガでわかるフーリエ解析
https://amzn.to/4azyPDj
手を動かしてまなぶ フーリエ解析・ラプラス変換
https://amzn.to/4bDH3vu
物理数学 量子力学のためのフーリエ解析・特殊関数
https://amzn.to/4bsygfU
単位が取れるフーリエ解析ノート
https://amzn.to/3V83fIl
今日から使えるフーリエ変換 普及版 式の意味を理解し、使いこなす
https://amzn.to/3ysbfvf
コメント