データに対して正弦波で近似を行いたい
显示 更早的评论
添付しているデータを正弦波で近似したいのですが，
アドオンのcurve fitting toolをつかってもうまくいきません．
具体的には，正弦波の和の近似を行うと振幅，位相，周期の項はわかりますが，定数項はわかりません．
また，カスタム式
y = f(x) = a*sin(b*x+c)+d (aは振幅，bは周期，cは位相，dは定数項)
で近似を行いましたが，うまくいきませんでした．
助言等よろしくお願いいたします．
なお，添付ファイルの1列目は時間のデータ，2，3はその時間に対する正弦波のデータです．
采纳的回答
Hernia Baby
2023-2-21
编辑：Hernia Baby
2023-2-21
ぱっとデータ見ました。
定数項はあらかじめ平均値をとって引くのはどうですか？
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 个评论
ちなみに定数項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
係数を抜き出します
cv = coeffvalues(fitresult);
一般モデルを作りそこに係数を与えていきます
ここで正規化のパラメタを norm からとってきています
f = fittype('a*sin(b*(x-mu)/std+c)+d')
c = cfit(f,cv(1),cv(2),cv(3),d,norm.mean,norm.std)
図示します
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
2023-2-22
编辑：Akira Agata
2023-2-22
+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
更多回答（1 个）
Hiro Yoshino
2023-2-22
3 个投票
初期値設定に問題があるのかなと思います。
アプリ > 詳細オプション > 係数の制約
バイアス値 = -0.003814 (-0.003834, -0.003795) << 信頼区間も出る
※ちなみに、データ B の方です。
追加でアドバイスをするとすれば ..... 残差のプロットを見るとまだ特徴が残っています。十分にデータを近似できてはいない状況だと思います。もう少しリッチなモデルを検討しても良いかと。
类别
在 帮助中心 和 File Exchange 中查找有关 線形回帰と非線形回帰 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
