Damped harmonic motion curve fit

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!

Answer 1

[Star Strider]

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:

[Stashu Kozlowski]

Hey, is there anyway I can use an .xls file instead?

[Star Strider]

Sure! I attached it to this Comment.

To be complete, the code I used to create it:

D = load('Stashu Kozlowski DHM.mat');
x = D.x;
y = D.y;
T = table(x, y);
writetable(T, 'Stashu Kozlowski DHM.xlsx')

and one way to read it:

T1 = readtable('Stashu Kozlowski DHM.xlsx');

The .mat file is easiest and most efficient way to store it and to read it.

[Mengyan Zhu]

Thank you for posting this answer! I am a beginner in MATLAB and I have a similar question. Only your answer greatly helpful to me. If you don't mind, may I ask you what are the meaings of b1, b2, b3 mean? How should I modify to make it become a cosine function like this: (t) = *e*** (2*pi*t/T +) and what is the last line of code mean?

[Star Strider]

Mengyan Zhu —

My pleasure!

You can change ‘fit’ to be whatever you want it to be, providing you update the rest of the code to match the number of parameters (if you change them).

In that event, ‘fit’ would become something like this:

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

The last line:

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')

is a text call that prints out the function and the parameter values on the plot. It is not necessary for the code, however the plot looks better with it (at least in my opinion). Run the code with it and without it (make it a comment) to see the difference. Note that you would have to change it to match the function in ‘fit’ for the parameters to make sense with it.

[Daniel Marti]

Hi Strider,

Very cool script to fit a damped sinusoidal.

I have two questions that I hope you can help me with:

-How would vary the code to fit a function with a larger 1st peak?

-How would I go about making the damping happen earlier?

Thank you in advance.

[Daniel Marti]

This is how I was able fit my data so far.

[Star Strider]

Daniel Marti —

Thank you!

It would be necessary to know the system that created that signal. It appears that a second exponential function would be necessary, however that is not as simple as adding another exponential term. There should be some sort of relationship between the first and second exponentials so that they are not simply modeling the same curve with essentially the same parameter. Even though a second exponential might be a better fit, the best fit would be to model the system that created those data. That way, the estimated parameters would have some physical meaning, beyond dimply fitting the curve.

[Daniel Marti]

Hi Strider,

Thank you for your answer.

The system that generated this data is basically the discharging of a big capacitor bank.

So there's obviously resistance, inductace and capacitance involved, thus making it an RLC series circuit.

[Daniel Marti]

I meant to say: *inductance.

I truly appreciate your input. Thanks again.

[Star Strider]

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’.

[Hossam Amin]

Hello Star Strider,

Thank you for posting this script. I have a similar damping measurement data ( torque and speed) which I have uploded, I am trying to utilize your script in order to fit the both speed and torque measuremenr dataset. But without success. Can you please help me with that? What do I need to do or change in the code in order to work?

Thank you in advance for your help.

[Star Strider]

I was able to get a decent fit to these data.

Please post this as a new Question and I will post the code I use here as an Answer to it.

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

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

.

[Star Strider]

Please post a link to your Question here.

[Hossam Amin]

Thanks Star Strider,

Here is the link of the question

https://www.mathworks.com/matlabcentral/answers/1634565-curve-fitting-of-large-data-measurement

Answer 2

[Alex Sha]

The fellow results are little better.

Parameter Best Estimate

---------- -------------

b1 -1.36099782974822

b2 -599.110824641553

b3 0.000259106153388041

b4 1.22915310606227

b5 0.0138196119722517