# Damped harmonic motion curve fit

45 views (last 30 days)
Stashu Kozlowski on 9 Oct 2019
Commented: Star Strider on 21 Oct 2020
Hey,
I have a data set in matlab, when plotted it looks like this:
My goal is to determen a damped sinusoidal equation that would fit this data set, I honestly dont even know how to start. I have included my code, but it isn't much.Any help is much appreciated. Thank you!
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Accepted Answer

Star Strider on 10 Oct 2019
Try this:
D = load('Stashu Kozlowski DHM.mat'); % File Attached
x = D.x;
y = D.y;
y = detrend(y); % 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
xp = linspace(min(x),max(x), 500);
figure
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')
text(0.3*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.0f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.0f%.3f)$', [s(1:2); 1./s(3:4)]), 'Interpreter','latex')
The estimated parameters are:
s =
-1.398211481931498e+00
-6.142349926864338e+02
2.591368008158479e-04
-5.442228857001487e+00
-3.075267405594925e-15
and the fit is nearly perfect:
##### 10 CommentsShowHide 9 older comments
Star Strider on 21 Oct 2020
My pleasure!
It shouldn’t be difficult to generate a time-domain model of your circuit. Use that to define ‘fit’ in my code, instead of the funciton I posted. It may be necessary to change the number or order of initial parameter estimates (the ‘[yr; -10; per; -1; ym]’ in the fminsearch call) to matrch those in your model.
An updated (and slightly more robust) version of ‘zci’ is:
zci = @(v) find(diff(sign(v)));
The rest of the code is unchanged, except as necessary to work with your expression of ‘fit’.

Sign in to comment.

### More Answers (1)

Alex Sha on 10 Oct 2019
The fellow results are little better.
Parameter Best Estimate
---------- -------------
b1 -1.36099782974822
b2 -599.110824641553
b3 0.000259106153388041
b4 1.22915310606227
b5 0.0138196119722517
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Community Treasure Hunt

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

Start Hunting!

Translated by