データに対して正弦波で近似を行いたい

添付しているデータを正弦波で近似したいのですが，

アドオンのcurve fitting toolをつかってもうまくいきません．

具体的には，正弦波の和の近似を行うと振幅，位相，周期の項はわかりますが，定数項はわかりません．

また，カスタム式

y = f(x) = a*sin(b*x+c)+d　　(aは振幅，bは周期，cは位相，dは定数項)

で近似を行いましたが，うまくいきませんでした．

助言等よろしくお願いいたします．

なお，添付ファイルの1列目は時間のデータ，2，3はその時間に対する正弦波のデータです．

Answer 1

[Hernia Baby]

ぱっとデータ見ました。

定数項はあらかじめ平均値をとって引くのはどうですか？

clc,clear,close all;
A = readmatrix('20230212-0.6mm.csv','NumHeaderLines',3);

ここで2列目に定数項を与えます

t = A(:,1);
x = A(:,2)+4;

定数項dを平均値として推定し、センタリングを行います。

x_mu    = mean(x);
x_1     = x -x_mu;
[fitresult, gof] = createFit(t, x_1);

以下はアプリで作成したものです。

function [fitresult, gof] = createFit(t, x_1)

%% 近似: 'MyPoly'。
[xData, yData] = prepareCurveData( t, x_1 );

% 近似タイプとオプションを設定します。
ft = fittype( 'sin1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [-Inf 0 -Inf];
opts.Normalize = 'on';
opts.StartPoint = [0.411512137125749 3.62760416984975 1.49630542841642];

% モデルをデータに近似します。
[fitresult, gof] = fit( xData, yData, ft, opts );

% データの近似をプロットします。
figure( 'Name', 'MyPoly' );
h = plot( fitresult, xData, yData);
legend( h, 'x_1 vs. t', 'MyPoly', 'Location', 'NorthEast', 'Interpreter', 'none' );
% ラベル Axes
xlabel( 't', 'Interpreter', 'none' );
ylabel( 'x_1', 'Interpreter', 'none' );
grid on
end

[Hernia Baby]

ちなみに定数項dも入れたいよって場合は以下をご参考ください

load("sample.mat")
t = A(:,1);
x = A(:,2) + 4;

定数項dを引いてフィッティングします

d = mean(x);
x_c = x - d;
[fitresult, norm, gof] = createFit(t, x_c);

式を見てみましょう

fitresult

fitresult = 
     General model Sin1:
     fitresult(x) =  a1*sin(b1*x+c1)
       where x is normalized by mean 0.3 and std 0.1732
     Coefficients (with 95% confidence bounds):
       a1 =     -0.5324  (-0.5328, -0.532)
       b1 =        2.74  (2.74, 2.741)
       c1 =       1.471  (1.47, 1.471)

係数を抜き出します

cv = coeffvalues(fitresult);

一般モデルを作りそこに係数を与えていきます

ここで正規化のパラメタを norm からとってきています

f = fittype('a*sin(b*(x-mu)/std+c)+d')

f = 
     General model:
     f(a,b,c,d,mu,std,x) = a*sin(b*(x-mu)/std+c)+d

c = cfit(f,cv(1),cv(2),cv(3),d,norm.mean,norm.std)

c = 
     General model:
     c(x) = a*sin(b*(x-mu)/std+c)+d
     Coefficients:
       a =     -0.5324
       b =        2.74
       c =       1.471
       d =       4.123
       mu =         0.3
       std =      0.1732

図示します

plot(c,t,x)

関数はこちら

function [fitresult, norm, gof] = createFit(t, x)

%% 近似: 'MyPoly'。
[xData, yData] = prepareCurveData(t, x);

% 加えた部分 ----
norm.mean = mean(xData);
norm.std  = std(xData);
% ---- 加えた部分

% 近似タイプとオプションを設定します。
ft = fittype( 'sin1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [-Inf 0 -Inf];
opts.Normalize = 'on';
opts.StartPoint = [0.411512137125749 3.62760416984975 1.49630542841642];

% モデルをデータに近似します。
[fitresult, gof] = fit( xData, yData, ft, opts );

% データの近似をプロットします。
figure( 'Name', 'MyPoly' );
h = plot( fitresult, xData, yData);
legend( h, 'x_1 vs. t', 'MyPoly', 'Location', 'NorthEast', 'Interpreter', 'none' );
% ラベル Axes
xlabel( 't', 'Interpreter', 'none' );
ylabel( 'x_1', 'Interpreter', 'none' );
grid on
end

[Akira Agata]

+1

fminsearch 関数をうまく使うと、以下のように計算することができます。

ちなみに下記で計算した結果、チャンネルAの定数項の推定値は cOpt(4) = 2.9939e-04 となります。

% データを読み込み
t = readtable('20230212-0.6mm.csv', ...
  VariableNamingRule = 'preserve');
x = t.("時間");
y = t.("チャンネルA");

% 近似式を y = f(x) = c(1)*sin(c(2)*x+c(3))+c(4) と想定
fcnChA = @(c) c(1)*sin(c(2)*x + c(3)) + c(4);

% 近似式との二乗誤差の総和
fcnErr = @(c) sum(abs(y - fcnChA(c)).^2);

% グラフより、おおよそ振幅0.6, 周期0.4であることから
% 初期値を以下のように設定
c0 = [0.6 2*pi/0.4 0 0];

% 二乗誤差の総和が最小となるよう c(1)~c(4)を最適化
cOpt = fminsearch(fcnErr, c0);

% 定数項 c(4) を表示
disp(cOpt(4))

% 最適化されたcによる近似値
yEst = fcnChA(cOpt);

% 結果を確認
figure
scatter(x, y)
hold on
plot(x, yEst, LineWidth=2)
legend({'Data', 'Estimated'})
grid on
box on

Answer 2

[Hiro Yoshino]

初期値設定に問題があるのかなと思います。

アプリで実行しているのは、基本的には @Akira Agata さんと同じ方法なので (最小二乗法)、アプリの初期値を見直してみましょう。

アプリ > 詳細オプション > 係数の制約

から開始点 (初期値) を設定できます。@Akira Agata さんの初期条件をそのままパクると、普通に収束しました。

バイアス値 = -0.003814 (-0.003834, -0.003795) << 信頼区間も出る

※ちなみに、データ B の方です。

追加でアドバイスをするとすれば ..... 残差のプロットを見るとまだ特徴が残っています。十分にデータを近似できてはいない状況だと思います。もう少しリッチなモデルを検討しても良いかと。

[優輔 井内]

Hiro Yoshino 様

ご助言いただきありがとうございます．

おっしゃる通り上限値下限値を設定してみると収束いたしました．

残差の周期性については波形に対して正弦波の和などの近似を行いアプローチしてみます．