Getting decrease in a sine wave as a function of time

I am trying to fit an exponential function to the peaks in a decreasing sinusoid:

I tried importing the values for time and position like this because the decimal places are marked with commas in my data that should be dots in MatLab:

% Open data file
fid = fopen('Data.csv');

% Read data file
readData = textscan(fid,'%f %f','HeaderLines',1,'Delimiter',',');

% Extract data from readData
time = readData{1,1}(:,1);
position = readData{1,2}(:,1);

% Plot data
plot (time,position,'bo-')

Is this correct? The plot looks odd. Attached Data.csv if one wanted to check. The file consists of two long columns of values.

Next I'm thinking of using the peakfinder (I am not sure about its input and output at all):

C = vertcat(time, position);
peakfinder(C)

And finally fit the values using polyfit? It should be in the form position = a * exp(-b * time) where a and b are constants. I think I should start something like this:

fit = polyfit(x, log(y), 1)

Where x and y for me are time and position from peakfinder, not sure how to write those at all.

[per isakson]

• Concerning the decimal separator see http://uk.mathworks.com/matlabcentral/answers/57276-import-data-file-with-comma-decimal-point#answer_69283. Your code doesn't read the file correctly with R2016a.
• The Import Data tool interprets comma as decimal separator correctly.

&nbsp

• "to the peaks in a decreasing sinusoid" &nbsp see env_secant(x_data, y_data, view, side) by Andreas Martin

[Star Strider]

Thank you, Per. I didn’t realise there is a comma decimal separator. I solved that problem with textscan, then strrep, and str2double.

[Daniel Holmberg]

Thank you Per! That information was exactly what I was looking for, although now I'm going with Star Strider's method down below since it contains exactly what I need (the exponential term) and a bit more. I'm using matlab for the first time and it seems very useful. Thanks guys for great answers!

Answer 1

[Star Strider]

See: How to filter noise from time-frequency data and find natural frequency of a cantilever?

The data you attached are not the data in the plot. Please attach the correct data.

[per isakson]

"attached are not the data in the plot" &nbsp I think it is.

[Star Strider]

@ Per — Correct. Fixed in my code here.

The Code:

fidi = fopen('Daniel Holmberg Data.csv','rt');
yy = textscan(fidi, '%s%s', 'Delimiter','\t', 'HeaderLines',1, 'CollectOutput',1);
yy = strrep(yy{:}, ',', '.');
y = str2double(yy(:,2));
t = str2double(yy(:,1));

yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zt = t(yz(:) .* circshift(yz(:),[1 0]) <= 0);                                   % Find zero-crossings
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,t) - y);                                                 % Least-Squares cost function
s = fminsearch(fcn, [yr; -1E-1;  per;  -1;  ym])                            % Minimise Least-Squares

xp = linspace(min(t),max(t), 2000);

figure(3)
plot(t,y,'b')
hold on
% plot(t,yf,'y', 'LineWidth',2)
plot(xp,fit(s,xp), 'r', 'LineWidth',1.5)
hold off
grid
title('Sample Data 2')
xlabel('Time')
ylabel('Amplitude')
legend('Original Data', 'Fitted Curve')
txtstr = sprintf('y(t) = %.3f\\cdot e^{%.3f\\cdott}\\cdotsin(%.3f\\cdotx%+.3f)%+.3f', [s(1:2); 2*pi./s(3:4); s(5)]);
text(11, 0.72, txtstr)

The Plot:

[Daniel Holmberg]

The "original data plot" and my attached plot looks identical to me so it should really be the corresponding data. But hey, thanks for the great answer! I'll have a thourough look at the code.

[Star Strider]

My pleasure!

They are the same. I originally had problems reading your file because of the comma decimal separator. It is a ‘.csv’ file, so that created problems for Excel here in the U.S. After Per clued me into that, I was easily able to read it with textscan and an extra line of code. (Saving it as a ‘.txt’ file would eliminate that problem, since it would not default (on my computer) to Excel.)

[Daniel Holmberg]

I see. Do you know by the way how to evaluate the goodness of the fit? I read this but don't even know how to define what they have as "curvefit" in your code.

[Star Strider]

I do not have the Curve Fitting Toolbox (I don’t need it with the other Toolboxes I have), so I cannot help you with its functions.

There are a number of ways to ‘measure’ goodness-of-fit, the most commonly used being the Coefficient of Determination. Parameter confidence intervals tell you if a particular parameter is necessary in the model. If the confidence limits of a parameter include zero (are of opposite signs), that parameter is not necessary in the model.

The Statistics and Machine Learning Toolbox function fitnlm and nlinfit do everything I need. To use nlinfit with my code, add these lines:

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(t,y,fit,s);
betaci = nlparci(beta,R,'covar',CovB);
beta_mtx = [betaci(:,1) beta betaci(:,2)]

and to use the fitnlm function:

mdl = fitnlm(t,y,fit,s)

mdl = 

Nonlinear regression model:
    y ~ b1*exp(b2*x)*(sin(2*pi*x/b3 + 2*pi/b4)) + b5

Estimated Coefficients:
          Estimate         SE         tStat       pValue   
          _________    __________    _______    ___________

      b1      0.16606     0.0019625     84.616              0
      b2    -0.026461    0.00055049    -48.069    1.5324e-281
      b3       1.2335    0.00013222     9329.1              0
      b4     -0.90296     0.0015152    -595.94              0
      b5      0.55719    0.00049126     1134.2              0

Number of observations: 1200, Error degrees of freedom: 1195
Root Mean Squared Error: 0.017
R-Squared: 0.935,  Adjusted R-Squared 0.935
F-statistic vs. constant model: 4.28e+03, p-value = 0

That indicates an extremely good fit.

Note that the parameters estimated as ‘s’ will differ slightly from ‘beta’ returned by the Statistics and Machine Learning Toolbox functions because fminsearch and nlinfit use different optimisation algorithms, and come to different conclusions as to what constitutes the best fit.

[Daniel Holmberg]

Okay thanks a lot it looks fairly simple. I couldn't try for myself because nlinfit has not yet been implemented to the statistics package from octave forge.

[Star Strider]

My pleasure.

Answer 2

[Roger Stafford]

The data you show should probably be fitted with something like:

a*exp(-b*time)*sin(c*time+d)

with the four parameters, a, b, c, and d, to be adjusted. It’s known as an exponentially decaying sinusoidal signal.

[Daniel Holmberg]

Yes, that is correct. The thing is though I am interested only in the decrease/exponential decay hence I only wrote the exponential term. But the answers in this thread were even better than expected since I now know the approximation of the whole function. Both the for me the important decrease and how the wave looks.