How to find and fit the objective function for a damped oscillation?

The excel file contains force data of a damped oscillation system. In order to find the natural frequency of the system I am trying to fit the graph with an objective function as follows.

I have defined the function of the osillation.

function F = forcesim(t,A,omega,phi,C,K)
F=A*cos(omega*t+phi)*exp(-t*C) + K;
end

I have tried with a code to find the estimated frequency as follows.

dataset = readtable ("F300.xlsx");
Force = dataset.Force2;
Time = dataset.Time2;

plot(Time,Force)
xlabel("Time (s)")
ylabel("Force (N)")
grid on

ind = Time>0.3 & Time<1.3;
[pks,locs] = findpeaks(Force(ind),Time(ind),'MinPeakProminence',1);

peaks = [pks locs]

% visualizse the selected peaks and label them.
hold on
plot(locs,pks,'v')
text(locs+.02,pks,num2str((1:numel(pks))'))
hold off

% calculate the approximate frequency
meanCycle = length(pks)/(locs(end)-locs(1))

Hence, that frequency is an estimated one, and therefore I am planning to find the frequency of the graph by matching the objective function mentioned above.

It would be really appreciated if you could tell me how to do that and find the frequency ω.

Answer 1

[Star Strider]

Using fminsearch —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1211598/F300.xlsx')

T1 = 3001×2 table
    Time2    Force2 
    _____    _______

        0    0.07471
    0.001    0.05859
    0.002    0.06006
    0.003    0.04541
    0.004    0.09375
    0.005    0.07178
    0.006    0.05713
    0.007    0.05859
    0.008    0.07031
    0.009    0.06006
     0.01    0.04541
    0.011     0.0293
    0.012    0.07031
    0.013    0.06299
    0.014    0.06299
    0.015    0.07031

x = T1.Time2;
y = T1.Force2;

start = find(y > 1, 1, 'first')                                                 % Find Offset

start = 160

y = y(start:end);
x = x(start:end);

y = detrend(y,1);                                                               % Remove Linear Trend
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                            % Returns Approximate Zero-Crossing Indices Of Argument Vector
zt = x(zci(y));
per = 2*mean(diff(zt));                                                         % Estimate period
ym = mean(y);                                                                   % Estimate offset

fit = @(b,x)  b(1) .* exp(b(2).*x) .* (sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5);   % Objective Function to fit
fcn = @(b) norm(fit(b,x) - y);                                                  % Least-Squares cost function
[s,nmrs] = fminsearch(fcn, [yr; -10;  per;  -1;  ym])                           % Minimise Least-Squares

s = 5×1

155.8223 -6.4293 0.0335 -1.4713 0.0000

nmrs = 61.2464

xp = linspace(min(x),max(x), numel(x)*10);

figure(3)
plot(x,y,'b', 'LineWidth',1.5)
hold on
plot(xp,fit(s,xp), '-r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Original Data',  'Fitted Curve')
xlim([min(x) 1])
text(0.3*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.1f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.0f%.3f)%+.1f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
format short eng
grid

mdl = fitnlm(x(:),y(:),fit,s)                                               % Get Statistics

mdl = 
Nonlinear regression model:
    y ~ b1*exp(b2*x)*(sin(2*pi*x/b3 + 2*pi/b4)) + b5

Estimated Coefficients:
           Estimate         SE          tStat       pValue
          __________    __________    __________    ______

    b1        155.82        1.6153        96.468      0   
    b2       -6.4293      0.041868       -153.56      0   
    b3      0.033533    7.6151e-06        4403.5      0   
    b4       -1.4713     0.0036676       -401.16      0   
    b5    2.2889e-16       0.03774    6.0647e-15      1   


Number of observations: 2842, Error degrees of freedom: 2837
Root Mean Squared Error: 1.15
R-Squared: 0.971,  Adjusted R-Squared 0.971
F-statistic vs. constant model: 2.37e+04, p-value = 0

.

[Piyumal Samarathunga]

Thank you very much for your time and kind consideration.

[Star Strider]

As always, my pleasure!

[Piyumal Samarathunga]

Hello, can you please let me know how I can use this method without using the detrend function? Thank you.

[Star Strider]

If your data have a mean of zero, and no trend in the baseline, the detrend function is not needed. If the data have a non-zero baseline, use the mean of the signal as the initial estimate for ‘b(5)’.

If there is a trend to the signal, it is possible to deal with it by adding a sixth parameter to the objective function to estimate the slope of the signal baseline as well:

fit = @(b,x) b(1) .* exp(b(2).*x) .* (sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5) + b(6).*x; % Objective Function to fit

and change the text call as well:

text(0.3*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.1f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.0f%.3f)%+.1f %+.1f\\cdot x$', [s(1:2); 1./s(3:4); s(5:6)]), 'Interpreter','latex')

It will be necessary to provide an estimate for the slope and add that to the end of the vector that is the second argument to the fminsearch call (the initial parameter estimate vector).

Adding an extra parameter to account for a sloping baseline could make the fit a bit more difficult, possibly requiring more than one attempt to get the fit correct. Adding a slope estimate to the ‘fit’ objective function should accommodate a non-zero and sloping baseline without needing to use the detrend function.

.

[Piyumal Samarathunga]

Thank you very much for your reply.

I have tryied with the code without using detrend function.

When I use the detrend function it gives a good match with the requency and the amplitute, but when I remove the detrend function, it does not give a good match for the amplitude and the exponential envilop.

Could you please let me know what is procedure to make the estimation better.

Thank you.

Data = readtable('2023_02_01_F300_NSH_d0_2.xlsx'); % Read data from the file

t = Data.T;
F = Data.F;

 region = t>=0.34 & t<=1; % 
 F=F(region);
 t=t(region);

start = find(F > 1, 1, 'first') % Find Offset

start = 3

Fu = max(F);
Fl = min(F);

Fr = (Fu-Fl); % Range of F
Fz = F-Fu+(Fr/2);

% Returns Approximate Zero-Crossing Indices Of Argument Vector
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); 
zt = t(zci(F));
per = 2*mean(diff(zt)); % Estimate period
Fm = mean(F); % Estimate offset

% Objective Function to fit
fit = @(b,t)  b(1) .* exp(b(2).*t) .* (sin(2*pi*t./b(3) + 2*pi/b(4))) + b(5) + b(6).*t;

% Least-Squares cost function
fcn = @(b) norm(fit(b,t) - F);

[s,nmrs] = fminsearch(fcn, [Fr; -10;  per;  -1;  Fm; 3.28]);

 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 25.530902 

tp = linspace(min(t),max(t), numel(t)*10);
figure(3)
plot(t,F,'b', 'LineWidth', 1)
hold on

plot(tp,fit(s,tp), '-r')
hold off
grid
xlabel('Time')
ylabel('Force')
legend('Original Data',  'Fitted Curve')
xlim([min(t) 1])
title ("d = 0 (Test 2)")

text(0.4*max(xlim),0.7*min(ylim), ...
    sprintf('$y = %.3f\\cdot e^{%.1f\\cdot t}\\cdot sin(2\\pi\\cdot t\\cdot %.0f%.3f)%+.1f %+.1f\\cdot t$', ...
    [s(1:2); 1./s(3:4); s(5:6)]), 'Interpreter','latex')

format short eng
grid on

mdl = fitnlm(t(:),F(:),fit,s) % Get Statistics

mdl = 
Nonlinear regression model:
    y ~ b1*exp(b2*t)*(sin(2*pi*t/b3 + 2*pi/b4)) + b5 + b6*t

Estimated Coefficients:
          Estimate         SE         tStat        pValue   
          _________    __________    ________    ___________

    b1       2173.8        175.27      12.403     6.9717e-32
    b2      -14.652       0.21338     -68.666    1.0188e-301
    b3      0.03038    2.9266e-05      1038.1              0
    b4     -0.08849    9.2642e-05     -955.18              0
    b5       3.3012      0.062133      53.132    7.6814e-240
    b6    -0.024256      0.089048    -0.27239        0.78541


Number of observations: 661, Error degrees of freedom: 656
Root Mean Squared Error: 0.434
R-Squared: 0.967,  Adjusted R-Squared 0.966
F-statistic vs. constant model: 4.75e+03, p-value = 0

[Star Strider]

Adding a second exponential term is an option —

Data = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1288995/2023_02_01_F300_NSH_d0_2.xlsx')

Data = 5000×2 table
      T        F   
    _____    ______

        0    0.1597
    0.001    0.1187
    0.002    0.1421
    0.003    0.1523
    0.004    0.1611
    0.005    0.1333
    0.006    0.1348
    0.007    0.1421
    0.008    0.1143
    0.009     0.123
     0.01    0.1187
    0.011     0.123
    0.012    0.1113
    0.013    0.1348
    0.014    0.1216
    0.015    0.1143

% Data = readtable('2023_02_01_F300_NSH_d0_2.xlsx'); % Read data from the file

t = Data.T;
F = Data.F;

 region = t>=0.34 & t<=1; % 
 F=F(region);
 t=t(region);

start = find(F > 1, 1, 'first') % Find Offset

start = 3

Fu = max(F);
Fl = min(F);

Fr = (Fu-Fl); % Range of F
Fz = F-Fu+(Fr/2);

% Returns Approximate Zero-Crossing Indices Of Argument Vector
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); 
zt = t(zci(F));
per = 2*mean(diff(zt)); % Estimate period
Fm = mean(F); % Estimate offset

% Objective Function to fit
fit = @(b,t)  b(1) .* exp(b(2).*t) .* (sin(2*pi*t./b(3) + 2*pi/b(4))) + b(5) .* exp(b(6).*t) .* (cos(2*pi*t./b(7) + 2*pi/b(8))) + b(9) + b(10).*t;

% Least-Squares cost function
fcn = @(b) norm(fit(b,t) - F);

[s,nmrs] = fminsearch(fcn, [Fr; -10;  per;  -1; randn(4,1);  Fm; 3.28]);

 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 22.348051 

tp = linspace(min(t),max(t), numel(t)*10);
figure(3)
plot(t,F,'b', 'LineWidth', 1)
hold on

plot(tp,fit(s,tp), '-r')
hold off
grid
xlabel('Time')
ylabel('Force')
legend('Original Data',  'Fitted Curve')
xlim([min(t) 1])
title ("d = 0 (Test 2)")

text(0.4*max(xlim),0.7*min(ylim), ...
    sprintf('$y = %.3f\\cdot e^{%.1f\\cdot t}\\cdot sin(2\\pi\\cdot t\\cdot %.0f%.3f)%+.1f %+.1f\\cdot t$', ...
    [s(1:2); 1./s(3:4); s(5:6)]), 'Interpreter','latex')

format short eng
grid on

mdl = fitnlm(t(:),F(:),fit,s) % Get Statistics

Warning: Some columns of the Jacobian are effectively zero at the solution, indicating that the model is insensitive to some of its parameters. That may be because those parameters are not present in the model, or otherwise do not affect the predicted values. It may also be due to numerical underflow in the model function, which can sometimes be avoided by choosing better initial parameter values, or by rescaling or recentering. Parameter estimates may be unreliable.

mdl = 
Nonlinear regression model:
    y ~ F(b,t)

Estimated Coefficients:
           Estimate         SE           tStat         pValue   
           _________    __________    ___________    ___________

    b1        2170.5    0.00091142     2.3814e+06              0
    b2       -14.648      0.019285        -759.58              0
    b3       0.03038    2.9362e-05         1034.6              0
    b4     -0.088487    9.3005e-05        -951.42              0
    b5       -1.3208      0.078404        -16.846     3.5996e-53
    b6      -0.51561         1.252       -0.41182         0.6806
    b7       -2312.8    5.0301e-07    -4.5979e+09              0
    b8       0.67552      0.095935         7.0414     4.8205e-12
    b9        2.0477      0.027837         73.563    8.4308e-319
    b10      0.45528       0.87088        0.52279         0.6013


Number of observations: 661, Error degrees of freedom: 655
Root Mean Squared Error: 0.434
R-Squared: 0.967,  Adjusted R-Squared 0.966
F-statistic vs. constant model: 3.8e+03, p-value = 0

However it does not appear to make a significant difference in the fit. (I did not change the text arguments to reflect the new ‘fit’ function, so that will not change.)

.

[Piyumal Samarathunga]

Thank you very much for your explanation.

By changing the initial parameter estimate vector, I could have good fittings, but sometimes it takes time to select proper values.

[Star Strider]

As always, my pleasure!

The Global Optimization Toolbox has functions such as ga that can make the parameter space search easier and more efficient.

Answer 2

[Matt J]

Using fminspleas from the File Exchange,

https://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

[t,F]= readvars ("https://www.mathworks.com/matlabcentral/answers/uploaded_files/1211598/F300.xlsx");

region=t>=0.2 &t<=0.8;

F=F(region);
t=t(region);


flist={@(p,x)cos(p(1)*t+p(2)).*exp(-t*p(3)),1  };

p0=[0.15,0,1];
[p,AK]=fminspleas(flist,p0,t,F);

Warning: Rank deficient, rank = 1, tol = 3.271538e-12.

Warning: Rank deficient, rank = 1, tol = 3.271538e-12.

Warning: Rank deficient, rank = 1, tol = 3.271538e-12.

Warning: Rank deficient, rank = 1, tol = 3.271538e-12.

f=@(t) AK(1)*cos(p(1)*t+p(2)).*exp(-t*p(3)) +AK(2);

plot(t,F,'o',t,f(t)); legend('Data','Fit')

[Piyumal Samarathunga]

Thank you very much for your reply. I really appreciate it.

[Piyumal Samarathunga]

Hello, can you plese let me know, how I can tune this method for different data set. Becasue I have tried it with other data set, but it did not give a good fit. Thank you.