このような方法が一般的な同定の手法か存じ上げませんが、矩形波を最小二乗法でフィッティングした例をご紹介いたします。
目的関数が非線形ですので初期値によってはうまくフィッティングできないかと思います。
多峰性の目的関数の最適化にはGlobal Optimization Toolboxがご利用いただけます。
また、信号解析にはSignal Processing Toolboxなどもご利用いただけます。
ご参考になれば幸いです。
% 時系列ベクトル
Fs = 10;
t = 0:1/Fs:20;
% 測定値
measured = sin(t-2) + sin(3*(t-2))/3 + sin(5*(t-2))/5 + sin(7*(t-2))/7 + sin(9*(t-2))/9;
% 矩形波のモデル式(xは振幅,角周波数,位相,オフセット)
model = @(x) x(1)*sign(sin(x(2)*t-x(3)))+x(4);
% 適当なしきい値で2値化(highのインデックスを取得)
index = measured>0.5;
% highとlowのメディアン値を取る
hi_med = median(measured(index))
lo_med = median(measured(~index))
% 変化点抽出
edge_points = find(diff(index)>0);
period = mean(edge_points(2:2:end)-edge_points(1:2:end-1))/Fs
phi = edge_points(1)/Fs;
% 矩形波のパラメータの初期値
x0 = [(hi_med-lo_med)/2 2*pi/period phi (hi_med-lo_med)/2+lo_med];
% 目的関数を定義
% 例として測定値とモデル式の2乗誤差を目的関数とする
% 目的関数が多峰性(multimodal)のため局所解に陥る、初期値x0が重要
fun = @(x) norm(model(x)-measured);
% 最適化計算
options = optimset('Display','iter','PlotFcns',@optimplotfval);
[x,fval,exitflag,output] = fminsearch(fun,x0,options)
% より高度な最適化の方法として
% Optimization ToolboxやGlobal Optimization Toolboxを使う方法もある
% 結果の可視化
figure, plot(t,[measured; model(x0); model(x)]);
legend('測定値','初期値','最適値');
xlabel('時刻 t'); ylabel('値'); title('矩形波のフィッティング');
最適化ループでの目的関数の値:
フィッティング結果:
