MATLAB,Python,Scilab,Julia比較 第5章 その89【複素フーリエ係数⑭】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia5-backnumber/

はじめに
登場人物
【再掲】プログラムフロー
Pythonコード
処理結果
まとめ

はじめに

任意波形から複素フーリエ係数抽出。
それを元に元波形を複素フーリエ級数で再現。
これを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ファイルはこれ。

wave.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数分の等差数列

C = np.zeros(N*2+1, dtype=complex) # C係数群格納用(±両方)

for n in range(-N,N+1):
    # 係数C_n算出
    # C[N]がC_0
    # C_n = (1/2π)∫f(x)exp(-inx)dx
    C[n+N] = fx@np.exp(-1j*n*x).T*dx/(2*np.pi)


Ns = [10,50,200]
fig = plt.figure()

for i in range(0,len(Ns)):
    NN = Ns[i]; # 今回のC_n項数

    # f(x)=Σ(C_n exp(inx))
    Fourier_series = np.zeros(points) # 合成関数を0初期化
    for n in range(-NN,NN+1):
        Fourier_series = Fourier_series+(C[n+N]*np.exp(1j*n*x))
    
    ax = fig.add_subplot(len(Ns), 1, i+1)
    
    # 元波形と複素フーリエ級数波形の表示
    ax.plot(x, fx,'-b',lw=3)
    ax.plot(x, Fourier_series.real,'-r',lw=2)
    ax.set_title('n=±{:d}'.format(NN),fontsize=22,fontweight="bold")
    ax.set_ylim(-0.1,1.1)
    ax.set_xlim(-np.pi,np.pi)
    ax.grid(linestyle='dotted')

plt.show()