MATLAB,Python,Scilab,Julia比較第5章その98【複素フーリエ係数(周期2L)⑦】

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

はじめに
登場人物
【再掲】周期2Lの複素フーリエ
【再掲】プログラムフロー
Juliaコード
処理結果
まとめ

はじめに

前回までで複素フーリエの任意周期化が完了。
今回はこれをJuliaで動作確認する。

登場人物

博識フクロウのフクさん

イラスト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

【再掲】周期2Lの複素フーリエ

太郎くん

まずは周期2Lの複素フーリエを再掲

複素フーリエ級数

\(
\displaystyle f(x)=\sum_{n=-\infty}^{\infty}C_n e^{i\frac{n\pi x}{L}}
\)

複素フーリエ級数

\(
\displaystyle C_n=\frac{1}{2L}\int_{-L}^{L}f(x)e^{-i\frac{n\pi x}{L}}\,dx
\)

【再掲】プログラムフロー

太郎くん

そしてプログラムフローを再掲

csvファイル読み込み
各種変数初期化
フーリエ係数算出
n=10,50,200のパターンでフーリエ級数で波形を合成
グラフにプロット

フクさん

今回はこれをJuliaで実現する。

Juliaコード

フクさん

使用するcsvファイルは以下。

wave.csv ダウンロード

フクさん

Juliaコードは以下。

using PyPlot
using CSV
using DataFrames
using Printf

N=1000;  # 係数算出項数(同定元波形のplotよりも少なく)
L=10;    # 周期/2
csv_data=CSV.read("wave.csv", DataFrame); # 同定波形読み込み
wave = csv_data[:,1];
points=length(wave); # 波形のplot数取得
fx=wave'; # 波形を行ベクトルへ
dx=2*L/points; # 1plotあたりのx軸幅
x=range(-L,L,length=points); # -π～+πの範囲で波形plot数分の等差数列

C = zeros(Complex{Float64},1,N*2+1); # C係数群格納用(±両方)

for n = -N:N
    # 係数C_n算出
    # C[N+1]がC_0
    # C_n = (1/2π)∫f(x)exp(-inxπ/L)dx
    C[n+N+1] = fx*exp.(-1im*n*x*pi/L)*dx/(2*L);
end

Ns = [10,50,200];

for i = 1:length(Ns)
    NN = Ns[i]; # 今回のC_n項数
    
    # f(x)=Σ(C_n exp(inxπ/L))
    Fourier_series = zeros(1,points); # 合成関数を0初期化
    for n = -NN:NN
        Fourier_series = Fourier_series+(C[n+N+1]*transpose(exp.(1im*n*x*pi/L)));
    end
    # 元波形と複素フーリエ級数波形の表示
    subplot(length(Ns),1,i);
    plot(x, fx,linewidth = 3);
    plot(x, real(Fourier_series)',color="r", linewidth = 2);
    title(@sprintf("n=±%d",NN), fontsize=22, fontweight="bold");
    ylim([-0.1,1.1]);
    xlim([-L,L]);
    grid();
end