Damped Cosine Wave Fitting

Hey everyone, I am having trouble fitting a damped cosine wave function form experimental data that I took using a Cavendish Balance. I am given the time(x axis) and mrads(yaxis). Here is the code presented along with attached data. I need to fit this data with this function: "theta + A .* exp(-b*t) * cos(omega*t + delta)" I have done research on the 'lsqcurvefit' function in matlab only I am having trouble. Some help would be awesome thanks!

 clc;clear;
filename = 'Clockwise1.xls';
time1 = xlsread(filename,'A:A');
mrads1 = xlsread(filename,'B:B');
figure(1)
plot(time1,mrads1,'r.')
title('Cavendish Balance Decay (CW)')
xlabel('Time (sec)')
ylabel('Boom Position (mrads)')
xlim([0 1700]);

Answer 1

[Star Strider]

This works:

[d,si,r] = xlsread('Rob Mullins Clockwise1.xls');

x = d(:,1);
y = d(:,2);

yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,[1 0]) <= 0);                                         % Find zero-crossings
per = 2*mean(diff(zx));                                                         % Estimate period
ym = mean(y);                                                                   % Estimate offset

fit = @(b,x)  b(1).*(cos(2*pi*x./b(2) + 2*pi/b(3))) .* exp(b(4).*x) + b(5);     % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                              % Sum-squared-error cost function
s = fminsearch(fcn, [yr;  per;  -1;  -1;  ym])                                  % Minimise Least-Squares

xp = linspace(min(x),max(x));

figure(1)
plot(x,y,'b', 'LineWidth',1)
hold on
plot(xp,fit(s,xp), '--r', 'LineWidth',1.25)
hold off
grid
text(450, 7.6, sprintf('%.2f\\cdotcos(2\\cdot\\pi\\cdot%.4f\\cdotx %+.2f\\cdot2\\cdot\\pi)\\cdote^{%.4f\\cdotx} + %.2f', s(1),1/s(2),1/s(3),s(4),s(5)))
xlabel(char(si(1)))
ylabel(char(si(2)))
legend('Data', 'Regression')

[DENESH KUMAR M]

Works perfectly .Thanks!

[uzzi]

I also tried this code but it is showing some errors. Can you pls take a look, @Star Strider?

x = readtable('t.txt');
y = table2array(readtable('S.txt'));
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   %Range of 'y'
yz = y - yu+(yr/2);
zx = x(yz .* circshift(yz,[1 0]) <= 0);                                         % Find zero-crossings
per = 2*mean(diff(zx));                                                         % Estimate period
ym = mean(y);                                                                   % Estimate offset
fit = @(b,x) b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) .* exp(b(4).*x ) + b(5);     % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                              % Sum squared error cost function
s = fminsearch(fcn, [(yr); per; -1; -1; (ym)])                                   % minimize least squares
xp = linspace(min(x),max(x));
figure(1)
plot(x,y, 'b' , 'LineWidth' ,1)
hold on
plot(xp,fit(s,xp), '--r' , 'LineWidth' ,1.25)
hold off
grid
text(450, 7.6, sprintf( '%.2f\\cdotcos(2\\cdot\\pi\\cdot%.4f\\cdotx %+.2f\\cdot2\\cdot\\pi)\\cdote^ {%.4f\\cdotx} + %.2f' , s(1),1/s(2),1/s(3),s(4),s(5)))
xlabel('time')
ylabel('var')
legend( 'Data' , 'Regression' )

[Les Beckham]

I only had to make two changes to your code to get it to run. See below.

x = table2array(readtable('t.txt')); % create a datetime array instead of a table
x = x.Second - x(1).Second; % convert to an x vector that is type double and starts at zero
y = table2array(readtable('S.txt'));
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   %Range of 'y'
yz = y - yu+(yr/2);
zx = x(yz .* circshift(yz,[1 0]) <= 0);                                         % Find zero-crossings
per = 2*mean(diff(zx));                                                         % Estimate period
ym = mean(y);                                                                   % Estimate offset
fit = @(b,x) b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) .* exp(b(4).*x ) + b(5);     % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                              % Sum squared error cost function
s = fminsearch(fcn, [(yr); per; -1; -1; (ym)])                                  % minimize least squares

s = 5×1

0.1670 0.5416 -1.2986 -0.0690 49.9887

xp = linspace(min(x),max(x));
figure(1)
plot(x,y, 'b' , 'LineWidth' ,1)
hold on
plot(xp,fit(s,xp), '--r' , 'LineWidth' ,1.25)
hold off
grid
text(450, 7.6, sprintf( '%.2f\\cdotcos(2\\cdot\\pi\\cdot%.4f\\cdotx %+.2f\\cdot2\\cdot\\pi)\\cdote^ {%.4f\\cdotx} + %.2f' , s(1),1/s(2),1/s(3),s(4),s(5)))
xlabel('time')
ylabel('var')
legend( 'Data' , 'Regression' )

[Les Beckham]

You probably want to change the x, y location in the text() function call. I tried to edit my comment to show that, but when I re-ran the code with that change, it gave me errors in the first two lines of code that I didn't even change.

[Star Strider]

This code really was not designed for your signal.

I made some changes to it, so it sort of works —

x = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1284080/t.txt');
y = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1284085/S.txt');
% x
x = datetime(x{:,1}, 'InputFormat','yyyy-MM-dd HH:mm:ss.SSS');
x = second(x)-second(x(1))

x = 805×1

0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900

y

y = 805×1

49.9720 49.9740 49.9800 49.9720 49.9770 49.9840 49.9730 49.9610 49.9470 49.9400

figure
plot(x, y)
grid

ym = mean(y);                                                                   % Estimate offset
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   %Range of 'y'
yz = y - yu+(yr/2);
zx = find(diff(sign(y-ym)));                                                      % Find zero-crossings
per = 2*mean(diff(x(zx)))                                                         % Estimate period

per = 0.5393

fit = @(b,x) b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) + b(4).*x + b(5);     % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                              % Sum squared error cost function
s = fminsearch(fcn, [(yr); per; -0.1; -0.1; (ym)])                                   % minimize least squares

s = 5×1

-0.0364 0.6024 -0.0990 0.0043 49.9711

xp = linspace(min(x),max(x));
figure(1)
plot(x,y, 'b' , 'LineWidth' ,1)
hold on
plot(xp,fit(s,xp), '--r' , 'LineWidth' ,1.25)
hold off
grid
text(1, 49.8, sprintf( 'f(x) = %.2f\\cdotcos(2\\cdot\\pi\\cdot%.4f\\cdotx %+.2f\\cdot2\\cdot\\pi)\\cdote^ {%.4f\\cdotx} + %.2f' , s(1),1/s(2),1/s(3),s(4),s(5)))
xlabel('time')
ylabel('var')
legend( 'Data' , 'Regression' )

.

[uzzi]

Thank you for your kindness and answer, @Les Beckham and @Star Strider. Sorry for troubling you again. I want to fit this data with this function: f(t)=[sin(2*pi*f*t) exp^(-lambda*t)]

f = frequency

Is it possible to change?

[Star Strider]

Try this —

x = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1284080/t.txt');
y = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1284085/S.txt');
% x
x = datetime(x{:,1}, 'InputFormat','yyyy-MM-dd HH:mm:ss.SSS');
x = second(x)-second(x(1))

x = 805×1

0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900

y

y = 805×1

49.9720 49.9740 49.9800 49.9720 49.9770 49.9840 49.9730 49.9610 49.9470 49.9400

% figure
% plot(x, y)
% grid

ym = mean(y);                                                                   % Estimate offset
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   % Range of 'y'
yz = y - yu+(yr/2);
zx = find(diff(sign(y-ym)));                                                    % Find zero-crossings
per = 2*mean(diff(x(zx)))                                                       % Estimate period

per = 0.5393

fit = @(b,x) b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) + x.*exp(b(4).*x) + b(5);    % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                              % Sum squared error cost function
opts = optimset('MaxFunEvals',5000, 'MaxIter',5000);
s = fminsearch(fcn, [(yr); per; -0.1; -5; (ym)], opts)                                % minimize least squares

s = 5×1

1.0e+03 * 0.0001 0.0005 -0.0001 -2.5713 0.0500

xp = linspace(min(x),max(x));
figure(1)
plot(x,y, 'b' , 'LineWidth' ,1)
hold on
plot(xp,fit(s,xp), '--r' , 'LineWidth' ,1.25)
hold off
grid
text(1, 49.8, sprintf( 'f(x) = %.2f\\cdotcos(\\cdot%.2f\\cdot\\pi\\cdotx %+.2f\\cdot\\pi) + x\\cdot e^{%.2f\\cdotx} %+.2f' , s(1),2/s(2),2/s(3),s(4),s(5)))
xlabel('time')
ylabel('var')
legend( 'Data' , 'Regression' )

This should actually be a new question, so I will delete the last few comments in this thread in a few days. They actually have nothing to do with the original Question.

.

[uzzi]

Thanks a lot for your answer

Would you mind explain these 3 lines for me?

text(1, 49.8, sprintf( 'f(x) = %.2f\\cdotcos(\\cdot%.2f\\cdot\\pi\\cdotx %+.2f\\cdot\\pi) + x\ \cdot e^{%.2f\\cdotx} %+.2f' , s(1),2/s(2),2/s(3),s(4),s(5)))
fit = @(b,x) b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) + x.*exp(b(4).* x) + b(5);    % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);

That will be very helpful. Thanks a lot :)

[Star Strider]

The first line prints the function text in the plot.

The second one is the objective function (describing the data) that is used to fit the data.

The third one is the function that fminsearch uses to fit the data by comparing the objective function to the data.

Thjese most recent Comments have no direct relation to the original post, so I will delete them in a few days. They should actually be a new Question.

.

[uzzi]

Thank you very much,