MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/
はじめに
の、
MATLAB,Python,Scilab,Julia比較 第5章 その47【フーリエ係数⑪】
を書き直したもの。
フーリエ係数に至る道。
今回はフーリエ係数を求めるプログラムをJuliaで実現。
【再掲】フーリエ係数を求めるプログラムフロー
まずは、プログラムフローを再掲。
- csvファイル読み込み
- 各種変数初期化
- フーリエ係数算出
- n=10,50,200のパターンでフーリエ級数で波形を合成
- グラフにプロット
今回はフーリエ係数を求めるプログラムをJuliaで実現。
Juliaコード
まず、使用する波形を取り込んだcsvファイル
Scilabコードは以下になる。
using PyPlot
using CSV
using DataFrames
using Printf
N=1000; # 係数算出項数(同定元波形のplotよりも少なく)
csv_data=CSV.read("wave.csv", DataFrame); # 同定波形読み込み
wave = csv_data[:,1];
points=length(wave); # 波形のplot数取得
fx=wave'; # 波形を行ベクトルへ
dx=2*pi/points; # 1plotあたりのx軸幅
x=range(-pi,pi,length=points); # -π~+πの範囲で波形plot数分の等差数列
a = zeros(1,N); # a係数群格納用
b = zeros(1,N); # b係数群格納用
for n = 1:N
# 係数a_n算出
# a_n = (1/π)∫f(fx)cos(nx)dx
a[n] = fx*cos.(n*x)*dx/pi;
# 係数b_n算出
# a_n = (1/π)∫f(fx)cos(nx)dx
b[n] = fx*sin.(n*x)*dx/pi;
end
# 係数a_0算出
a0=sum(fx)*dx/pi;
Ns = [10,50,200];
for i = 1:length(Ns)
NN = Ns[i]; # 今回のa_n,b_n項数
# f(x)=a_0+Σ(a_n cos(nx)+ b_n sin(nx))
Fourier_series=ones(1,points)*a0/2;
for n = 1:NN
Fourier_series = Fourier_series+(a[n]*cos.(n*x)'+b[n]*sin.(n*x)');
end
# 元波形とフーリエ級数波形の表示
subplot(length(Ns),1,i);
plot(x, fx,linewidth = 3);
plot(x, Fourier_series',color="r", linewidth = 2);
title(@sprintf("n=%d",NN));
ylim([-0.1,1.1]);
xlim([-pi,pi]);
grid();
end
処理結果
処理結果は以下。
考察
CSVファイルの読み込みとDataFrames扱うパッケージが必要となる。
よって、事前に以下を実施しておく必要がある。
import Pkg;
Pkg.add("CSV");
Pkg.add("DataFrames");
コード自体はほぼMATLABと一緒だが、
ベクトルを扱うため、各所でドット演算子が必要になる点に注意。
まとめ
- フーリエ係数を求めるプログラムをJuliaで実現。
- おおよそ元の波形を再現できる係数が算出できている。
- 不連続点では流石に振動している。
MATLAB、Python、Scilab、Julia比較ページはこちら
マンガでわかるフーリエ解析
Amazon.co.jp
手を動かしてまなぶ フーリエ解析・ラプラス変換
Amazon.co.jp
物理数学 量子力学のためのフーリエ解析・特殊関数
物理数学 量子力学のためのフーリエ解析・特殊関数 | 柴田 尚和, 是常 隆 | 数学 | Kindleストア | Amazon
Amazonで柴田 尚和, 是常 隆の物理数学 量子力学のためのフーリエ解析・特殊関数。アマゾンならポイント還元本が多数。一度購入いただいた電子書籍は、KindleおよびFire端末、スマートフォンやタブレットなど、様々な端末でもお楽しみい...
単位が取れるフーリエ解析ノート
https://amzn.to/3V83fIl
今日から使えるフーリエ変換 普及版 式の意味を理解し、使いこなす
https://amzn.to/3ysbfvf
コメント