MATLAB,Python,Scilab,Julia比較 第5章 その44【フーリエ係数⑧】

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

はじめに
登場人物
【再掲】フーリエ係数を求めるプログラムフロー
MATLABコード
処理結果
考察
まとめ

はじめに

フーリエ係数に至る道。
今回はフーリエ係数を求めるプログラムをMATLABで実現。

登場人物

博識フクロウのフクさん

イラスト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のパターンでフーリエ級数で波形を合成
グラフにプロット

フクさん

今回はフーリエ係数を求めるプログラムをMATLABで実現。

MATLABコード

フクさん

まず、使用する波形を取り込んだcsvファイル

wave.csv ダウンロード

フクさん

MATLABコードは以下になる。

N=1000;  % 係数算出項数(同定元波形のplotよりも少なく)
wave=csvread('wave.csv'); % 同定波形読み込み
points=length(wave); % 波形のplot数取得
fx=wave'; % 波形を行ベクトルへ
dx=2*pi/points; % 1plotあたりのx軸幅
x=linspace(-pi,pi,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);
    hold on
    plot(x, fx,'LineWidth',3);
    plot(x, Fourier_series,'-r','LineWidth',2);
    title(sprintf('n=%d',NN));
    ylim([-0.1,1.1]);
    xlim([-pi,pi]);
    grid();
end