charter(JJ1PID)

MATLABで非線形曲線近似する

非線形曲線近似する最低限のソースコードを載せました

Case.1 減衰振動のデータの近似

 y = 12.5e^{-0.2t}\sin(t)で減衰振動する物体の時刻tでの位置yの計測値、tdataとydataを考えます


MATLAB 非線形近似 減衰振動

%  tdata: 0から20まで100等分した値
tdata = linspace(0, 20, 100);

% ydata: 計測して得られた物体の位置(を想定)
% ここでは理想のモデル式に乱数のノイズを足しています
ydata = 12.5*exp(-0.2*tdata) .* sin(tdata) + 0.3*randn(size(tdata));


% 係数ベクトルtと、tdataを入力とする関数を作成
% 考えているモデル式の係数の部分を係数で置き換えます
fun = @(t, tdata) t(1)*exp(t(2)*tdata) .* sin(t(3)*tdata);

% 係数ベクトルxの各要素の初期値(なんとなく予想できる値)
% 近ければOKです
t0 = [12, -0.1, 1.1];

% 各係数が出力されます
t = lsqcurvefit(fun, t0, tdata, ydata);
% t=  [12.4747   -0.1979    0.9979]だったとします

% 計測値を黒の丸で、近似曲線を青の実践でプロットします
% plot(tdata, ydata, 'ko' ,tdata, fun(t, tdata), 'b-');
% なら計測値と近似曲線を一気にプロット出来ますが、近似曲線がtdataの値ごとの点を結んでいるためカクカクになります
% そこで、近似曲線の式を別に新しく定義し、重ねます
plot(tdata, ydata, 'ko');
hold on;
f = @(t)  12.4747* exp(-0.1979* t) .* sin(0.9979* t);
fplot(f, [0 20], 'b-');

% 以下、包絡線のプロットを付け加えています
f1 =@(t) 12.4747 * exp(-0.1979*t);
f2 =@(t) -12.4747 * exp(-0.1979*t);
fplot(f1, [0 20]);
fplot(f2, [0 20]);
hold off;


あとはプロパティインスペクターでグリッド入れたり、ラベルやタイトルを入れます


Case.2 LPFの利得の周波数特性

利得の周波数特性が G = \frac{1000}{j\omega+1000}と表されるLPFについて、その計測値の近似曲線を求めます
dB単位では G = 20\log_{10}{\frac{1000}{j\omega+1000}}\,\rm[dB]となります


LPF 利得 周波数特性

% 10から10^6までログスケールで等間隔な100個の点ω [rad/s]
omega = logspace(1,6,100)

% Gdata: 計測して得られた利得(を想定)
% ここでは理想のモデル式に乱数のノイズを足しています
Gdata = 20*log10(abs(1000./(i*omega+1000))) + randn(100,1)


G = @(k,omega) 20*log10(abs(k(1)./(i*omega + k(2))))


% 実際にフィッティングするときは大雑把でいいです
k0 = [990,1020]
k = lsqcurvefit(G,k0,omega,Gdata)


% 片対数グラフを作成
% 今回は近似曲線が緩やかで直線部分が長いのでカクカクしません
semilogx(omega,Gdata,'ko',omega,G(k,omega),'b-')

おわりに

上の例では標本値を作るところから始めていますが、実際にはデータをテーブル形式でインポートし(EPLとする)、次のようにすればいいです

xData = EPL.x
yData = EPL.y


もっと詳しく知りたい方はこちらへ
非線形曲線近似 (データ適合) 問題を最小二乗近似的に解く - MATLAB lsqcurvefit - MathWorks 日本


線形・多項式モデルの近似はこちらへ
www.charter-blog.com