MATLAB Answers

fitting scatter data into multiple cosine functions

3 ビュー (過去 30 日間)
an gao
an gao 2019 年 10 月 15 日
回答済み: Alex Sha 2019 年 10 月 15 日
hello everyone,
i am trying to fit a function:
y = 2*a^2*cos(w*t)+2*b^2*cos(w*t+2*phi)+2*a*b*cos(delta)*cos(w*t+phi);
with y and t known, and try to get the unknown parameters of a, b,w,phi and delta. the most important parameter is phi, the rest are less important.
the scatter data looks like as follows and is also attached in the matlab.mat file.
My question is that is it really possible to a get an unique/accurate result of the fitting?
thanks.
untitled.bmp

  0 件のコメント

サインインしてコメントする。

回答 (3 件)

Alex Sha
Alex Sha 2019 年 10 月 15 日
編集済み: Alex Sha 2019 年 10 月 15 日
Hi, how about the results below, the fitting is perfect, the reason for you cann't get good result is more possible the bad initial start valus. One more thing is the results are not unique.:
1:
Root of Mean Square Error (RMSE): 5.86503133415183E-10
Sum of Squared Residual: 6.91411710266713E-17
Correlation Coef. (R): 1
R-Square: 1
Adjusted R-Square: 1
Determination Coef. (DC): 1
Chi-Square: 1.19308902605578E-18
Parameter Best Estimate
---------- -------------
a 1.24243011796252
w -2000.00000005105
b -0.306494082231627
phi 0.448966482806276
delta 0.369303537459062
2:
Root of Mean Square Error (RMSE): 3.17092533626903E-9
Sum of Squared Residual: 2.02100826512676E-15
Correlation Coef. (R): 1
R-Square: 1
Adjusted R-Square: 1
Determination Coef. (DC): 1
Chi-Square: -1.07854914106232E-15
F-Statistic: 1.58977775365207E19
Parameter Best Estimate
---------- -------------
a 1.08573520679119
w 2000.0000002083
b -0.131006405724609
phi 3.73119497513834
delta 30.8055670286459
3:
Root of Mean Square Error (RMSE): 8.08540777797663E-10
Sum of Squared Residual: 1.31401376061692E-16
Correlation Coef. (R): 1
R-Square: 1
Adjusted R-Square: 1
Determination Coef. (DC): 1
Chi-Square: 2.45286536720573E-17
F-Statistic: 2.44514484842302E20
Parameter Best Estimate
---------- -------------
a 0.0596400967434023
w 2000.00000007458
b -1.14950810470384
phi -12.5353535070843
delta -17.9045320865089

  1 件のコメント

an gao
an gao 2019 年 10 月 15 日
thanks Alex, unfortunately the results are not what i am looking for. Is there a way to get a unique result?

サインインしてコメントする。


David Wilson
David Wilson 2019 年 10 月 15 日
編集済み: David Wilson 2019 年 10 月 15 日
Yes, but it is always tricky to fit such sinusoids without good starting estimates. I'm going to use lsqcurvefit from the optimisation toolbox.
%% fit a sinusoid
load matlab
plot(t,y)
It may well be prudent to start with a simplified version first. The trick is to get the dominant amplitude and frequency starting estimates to be reasonable. One way is simply to measure the time between peaks and use that for the period.
[PKS,LOCS]= findpeaks(y)
kPeriod = mean(diff(LOCS))
Ts = mean(diff(t));
Period = kPeriod*Ts;
w = 2*pi/Period
That gives us a frequency of about 2000. (I'm assuming t is regularly sampled, and you have findpeaks from the signal processing toolbox. If not, you just need to measure it yourself.)
The amplitude is about 3, so 3 = 2*a^2.
% Simplified version
f = @(p,t) 2*p(1)^2*cos(p(2)*t);
a = sqrt(3/2); w = 2000;
p0 = [a, w];
p = lsqcurvefit(f, p0, t,y)
ti = linspace(0,0.02)';
yi = f(p,ti);
y0 = f(p0, ti);v% initial guess for comparison
plot(t,y,'o', ti, yi,'-', ti, y0, '--')
Now we will try your (horribly overdone and over-parameterised version)
%% function of interest
f = @(p,t) 2*p(1)^2*cos(p(3)*t)+2*p(2)^2*cos(p(3)*t+2*p(4))+2*p(1)*p(2)*cos(p(5))*cos(p(3)*t+p(4));
a = sqrt(3/2); b = 0; w = 2000; phi = 0; delta = 0;
p0 = [a, b,w, phi, delta];
p = lsqcurvefit(f, p0, t,y)
ti = linspace(0,0.02)';
yi = f(p,ti);
y0 = f(p0, ti)
plot(t,y,'o', ti, yi,'-', ti, y0, '--')
results look OK.

  0 件のコメント

サインインしてコメントする。


Alex Sha
Alex Sha 2019 年 10 月 15 日
So, an gao, what is the result you want? Unless you add range constraint for all parameters (except w), because there are trigonometric functions in the fit formula.

  0 件のコメント

サインインしてコメントする。

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by