How to curve fitting a cosine function with exponential decay tau?

Hello,

I am new to MATLAB and I dont know how to fit a consine function. I have my data points and this is the function that I need to fit. I want to know the values of T, phi, and tau. How can I do that? If anyone could provide some insights, that would be greatly appreiciated!

This is my data for time and amplitude

x =

    0.1167
    0.2833
    0.3333
    0.3833
    0.4250
    0.4500
    0.4917
    0.5167
    0.5500
    0.6000
    0.6333
    0.6667
    0.7167
    0.7500
    0.8167
    0.9000
    1.0667
    1.1500
    1.2167
    1.2667
    1.3000
    1.3333
    1.3667
    1.4000
    1.4333
    1.4667
    1.5167
    1.5500
    1.5833
    1.6167
    1.7000
    1.8500
    1.9667
    2.0167
    2.0667
    2.1000
    2.1500
    2.1833
    2.2167
    2.2500
    2.2833
    2.3333
    2.3667
    2.4000
    2.4500
    2.5167
    2.6833
    2.7167
    2.7500
    2.8167
    2.8833
    2.9333
    2.9667
    3.0000
    3.0500
    3.0833
    3.1167
    3.1500
    3.2000
    3.2333
    3.2667
    3.3333
    3.4000
    3.5333
    3.6167
    3.6667
    3.7333
    3.7833
    3.8167
    3.8500
    3.8833
    3.9333
    3.9667
    4.0000
    4.0500
    4.0833
    4.1333
    4.2167
    4.3500
    4.6500
    4.5500
    4.6000
    4.6500
    4.6833
    4.7167
    4.7500
    4.8000
    4.8333
    4.8833
    4.9167
    4.9500
    5.0167
    5.0833
    5.2000
    5.2667
    5.3167
    5.4000
    5.4333
    5.4833
    5.5167
    5.5500
    5.5917
    5.6333
    5.6833
    5.7167
    5.7667
    5.8167
    5.9000
    6.0167
    6.1333
    6.2000
    6.2500
    6.3000
    6.3333
    6.3833
    6.4167
    6.4667
    6.5000
    6.5500
    6.5833
    6.6167
    6.6833
    6.8333
    6.9167
    6.9833
    7.0500
    7.0833
    7.1333
    7.1667
    7.2167
    7.2500
    7.2833
    7.3500
    7.3833
    7.4500
    7.5000
    7.5833
    7.6833
    7.7667
    7.8333
    7.8833
    7.9500
    7.9833
    7.0333

>> y

y =

    1.3963
    1.2217
    1.0472
    0.8727
    0.6981
    0.5236
    0.3491
    0.1745
         0
   -0.1745
   -0.3491
   -0.5236
   -0.6981
   -0.8727
   -1.0472
   -1.2217
   -1.2217
   -1.0472
   -0.8727
   -0.6981
   -0.5236
   -0.3491
   -0.1745
         0
    0.1745
    0.3491
    0.5236
    0.6981
    0.8727
    1.0472
    1.2217
    1.3963
    1.2217
    1.0472
    0.8727
    0.6981
    0.5236
    0.3491
    0.1745
         0
   -0.1745
   -0.3491
   -0.5236
   -0.6981
   -0.8727
   -1.0472
   -1.2217
   -1.3090
   -1.2217
   -1.0472
   -0.8727
   -0.6981
   -0.5236
   -0.3491
   -0.1745
         0
    0.1745
    0.3491
    0.5236
    0.6981
    0.8727
    1.0472
    1.2217
    1.3090
    1.2217
    1.0472
    0.8727
    0.6981
    0.5236
    0.3491
    0.1745
         0
   -0.1745
   -0.3491
   -0.5236
   -0.6981
   -0.8727
   -1.0472
   -1.1345
   -1.0472
   -0.8727
   -0.6981
   -0.5236
   -0.3491
   -0.1745
         0
    0.1745
    0.3491
    0.5236
    0.6981
    0.8727
    1.0472
    1.2217
    1.2392
    1.2217
    1.0472
    0.8727
    0.6981
    0.5236
    0.3491
    0.1745
         0
   -0.1745
   -0.3491
   -0.5236
   -0.6981
   -0.8727
   -1.0472
   -1.2043
   -1.0472
   -0.8727
   -0.6981
   -0.5236
   -0.3491
   -0.1745
         0
    0.1745
    0.3491
    0.5236
    0.6981
    0.8727
    1.0472
    1.2217
    1.2217
    1.0472
    0.8727
    0.6981
    0.5236
    0.3491
    0.1745
         0
   -0.1745
   -0.3491
   -0.5236
   -0.6981
   -0.8727
   -1.0472
   -1.1868
   -1.0472
   -0.8727
   -0.6981
   -0.5236
   -0.3491
   -0.0789

>> 

[Alex Sha]

Hi, are you sure two pairs of data "4.6500 -1.0472" and "7.0333 -0.0789" are correct?

Answer 1

[Star Strider]

Try this:

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(x./b(2)) .* (cos(2*pi*x./b(3) + 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
ylim(ylim*1.5)
xlabel('Time')
ylabel('Amplitude')
legend('Original Data',  'Fitted Curve')
text(0.2*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{\\frac{t}{%.3f}}\\cdot cos(2\\pi \\frac{t}{%.3f} %+.3f)$', [s(1:2); 1./s(3:4)]), 'Interpreter','latex', 'FontSize',14)

producing:

It is essentially the code in my previous Answer that you already discovered,simply changed to fit your equation.

[Mengyan Zhu]

Thank you very much! Do you know how can I add x and y error bars for the points?

[Star Strider]

My pleasure!

If you have error bars associated with them, use the errorbar function to plot the error bars. Recent releases can plot both x and y error bars, earlier releases can only plot y error bars. (I do not remember the release the x error bars were introduced.)