Damped harmonic motion curve fit

45 views (last 30 days)
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!

Accepted Answer

Star Strider
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:
Damped harmonic motion curve fit - 2019 10 09.png
  10 Comments
Star Strider
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
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

Community Treasure Hunt

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

Start Hunting!

Translated by