Clear Filters
Clear Filters

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

127 views (last 30 days)
優輔 井内
優輔 井内 on 21 Feb 2023
Commented: 優輔 井内 on 24 Feb 2023
添付しているデータを正弦波で近似したいのですが,
アドオンのcurve fitting toolをつかってもうまくいきません.
具体的には,正弦波の和の近似を行うと振幅,位相,周期の項はわかりますが,定数項はわかりません.
また,カスタム式
y = f(x) = a*sin(b*x+c)+d  (aは振幅,bは周期,cは位相,dは定数項)
で近似を行いましたが,うまくいきませんでした.
助言等よろしくお願いいたします.
なお,添付ファイルの1列目は時間のデータ,2,3はその時間に対する正弦波のデータです.

Accepted Answer

Hernia Baby
Hernia Baby on 21 Feb 2023
Edited: Hernia Baby on 21 Feb 2023
ぱっとデータ見ました。
定数項はあらかじめ平均値をとって引くのはどうですか?
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
  2 Comments
Hernia Baby
Hernia Baby on 21 Feb 2023
ちなみに定数項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
Akira Agata on 22 Feb 2023
Edited: Akira Agata on 22 Feb 2023
+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

Sign in to comment.

More Answers (1)

Hiro Yoshino
Hiro Yoshino on 22 Feb 2023
初期値設定に問題があるのかなと思います。
アプリで実行しているのは、基本的には @Akira Agata さんと同じ方法なので (最小二乗法)、アプリの初期値を見直してみましょう。
アプリ > 詳細オプション > 係数の制約
から開始点 (初期値) を設定できます。@Akira Agata さんの初期条件をそのままパクると、普通に収束しました。
バイアス値 = -0.003814 (-0.003834, -0.003795) << 信頼区間も出る
※ちなみに、データ B の方です。
追加でアドバイスをするとすれば ..... 残差のプロットを見るとまだ特徴が残っています。十分にデータを近似できてはいない状況だと思います。もう少しリッチなモデルを検討しても良いかと。
  1 Comment
優輔 井内
優輔 井内 on 24 Feb 2023
Hiro Yoshino
ご助言いただきありがとうございます.
おっしゃる通り上限値下限値を設定してみると収束いたしました.
残差の周期性については波形に対して正弦波の和などの近似を行いアプローチしてみます.

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!