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

15 views (last 30 days)
Piyumal Samarathunga on 28 Nov 2022
Commented: Star Strider on 3 Feb 2023 at 12:02
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.
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 ω.

Star Strider on 29 Nov 2022
Using fminsearch
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
.
Star Strider on 3 Feb 2023 at 12:02
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.
.

Matt J on 28 Nov 2022
Edited: Matt J on 28 Nov 2022
Using fminspleas from the File Exchange,
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')
##### 2 CommentsShowHide 1 older comment
Piyumal Samarathunga on 3 Feb 2023 at 10:15
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.

### Categories

Find more on Signal Operations in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!