Curve Fitting of large Data Measurement?

I have a measurement data which I have uploded, and I am trying to utilize a function to fit the measurement with a damping sinusoidal order of the measurement dataset. But without success. Can anyone please point me in the right direction? The curve fitting app in matlab is unable to make it. I tried a damping function in it, it didn't work. Also the lsqcurvefit nothing is working. I tried a sine function with not success too.

Thanks in advance.

[Ahmed raafat]

I have download your data , and plot the second and third col

the repeating sequence leads to failing the curve fitting like this (Untitled.png)

you need to add another factor

[Hossam Amin]

Maybe I forgot to mention. The data is composed of speed and torque measurements.

You would need to add the x-axis at time using linspace command.

Answer 1

[Star Strider]

Try this —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/871405/Measurement%20Data.xlsx', 'VariableNamingRule','preserve');
format longE
x = T1.Time;
yc(:,1) = T1.('Torque[Nm]');
yc(:,1) = smoothdata(yc(:,1), 'sgolay', 250);                                       % Remove Some Noise To Make The Fit Easier
yc(:,2) = T1.('RPM[1/min]');

for k = 1:size(yc,2)
    yk = yc(:,k);
    y = yk;
    % y = detrend(y);
    ym = mean(y);                                                                   % Estimate offset
    y = y - ym;
    yu = max(y);
    yl = min(y);
    yr = (yu-yl);                                                                   % Range of ‘y’
    yz = y-yu+(yr/2);
    zci = @(v) find(diff(sign(v)));                                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector
    zt = x(zci(y));
    per = 2*mean(diff(zt));                                                         % Estimate period

    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) - yk);                                                 % 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,yk,':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.1*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.3E\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3E%+.3f) %+.4f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
end

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

s = 5×1

1.0e+00 * 2.054070381848720e-01 -1.211851657300356e-04 3.381604644061252e+03 -1.273801702648572e+00 -6.517492492161443e-03

nmrs = 
     7.253749849224849e+00

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

s = 5×1

1.0e+00 * 5.565560752001775e+02 -3.434517388063775e-04 3.862185655766382e+03 -6.747881574969156e-01 3.384702768264112e+02

nmrs = 
     6.954008038834750e+03

Because it uses the least-square approach, the fit is dominated by the highest peaks, and fits them preferentially. Experiment with tweaking ‘fit’ to get different results. (This uses a slightly-modified version of my original code to process and fit the data.)

.

[Hossam Amin]

Thanks for the code.

Yeah I noticed that you used a smoothing function to filter the signal, and then a for loop that goes through the torque and speed measurement one after the other.

You didn't change the 'fit' function though, which begs to the question, what if I had a different measurement like the the one I newly uploaded, do I need to change the cost function (fminsearch) or the filter or the 'fit', it's also damping but the behavior is slighly different.

I also tried to chane the x-axis to represent it in seconds using the linspace instead of milisecond that exist now, it the it failed to plot correctly.

Thanks for the help. I appreciate it.

[Star Strider]

As always, my pleasure!

I changed the code slightly and created a separate variable ‘s0’ as the initial parameter estimates. Changing the second element to allowed a decent fit.

The smoothing was only to do the initial approximation, since the period (frequency) calculations are based on the signal zero-crossings, and excessive noise creates problems in that step. The actual fitting is done on the original signal, and the original signal is plotted with the fitted function.

Changing the ‘Time’ vector from milliseconds to seconds is straightforward, and produces a decent fit with the same initial parameter estimates and correspondingly different fitted parameters. I also retreved the variable names to use in various places in the code, including titling the plots.

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/871810/Measurement%20Data%20_2.xlsx', 'VariableNamingRule','preserve');

vn = T1.Properties.VariableNames;                                                   % Variable Names
ttlstr = vn([2 3]);                                                                 % Variable Names Used In The Code
x = T1.Time  * 1E-3;                                                                % Change 'Time' From Milliseconds To Seconds
yc(:,1) = T1.(ttlstr{1});
yc(:,1) = smoothdata(yc(:,1), 'sgolay', 250);                                       % Remove Some Noise To Make The Fit Easier
yc(:,2) = T1.(ttlstr{2});

for k = 1:size(yc,2)
    yk = yc(:,k);
    y = yk;
    % y = detrend(y);
    ym = mean(y);                                                                   % Estimate offset
    y = y - ym;
    yu = max(y);
    yl = min(y);
    yr = (yu-yl);                                                                   % Range of ‘y’
    yz = y-yu+(yr/2);
    zci = @(v) find(diff(sign(v)));                                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector
    zt = x(zci(y));
    per = 2*mean(diff(zt));                                                         % Estimate period

    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) - yk);                                                 % Least-Squares cost function
    s0 = [yr; -1E-2;  per;  -1;  ym];                                               % Initial Parameter Estimates
    [s,nmrs] = fminsearch(fcn, s0)                                                  % Minimise Least-Squares

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

    figure
    plot(x,yk,':b', 'LineWidth',1.5)
    hold on
    plot(xp,fit(s,xp), '--r')
    hold off
    grid
    xlabel('Time (s)')
    ylabel('Amplitude')
    title(ttlstr{k})
    legend('Original Data',  'Fitted Curve')
    text(0.1*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.3f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3f%+.3f) %+.4f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
end

s = 5×1

0.1789 -0.1048 3.1084 -1.2119 -0.0062

nmrs = 6.6842

s = 5×1

436.6429 -0.3349 3.5279 -2.0160 274.6755

nmrs = 6.5964e+03

This version should be reasonably robust to any other data provided to it.

.

[Hossam Amin]

That's great yes.

One last inquiry please, if I needed to determine the height "ampliude" of each oscillation peak and the distance (x-displacement) between each two peaks, in which they are all visible on the plot, how can I do so using findpeaks or any other function. I tried using 'findpeaks' on the torque measurement but it doesn't do the job perfectly.

On the speed measurement, I was able to get the peaks of the first 4. But not the distances between them.

Thanks in advance.

The code I used. Currently it is set to read the torque measurement only.

z=importdata('Measurement Data_Torque.xlsx');
ind = 1:15000;
Time = z.data(ind,1);  
y_meas = z.data(ind,2);  
[Ypk,Xpk,Wpk,Ppk] = findpeaks(y_meas(ind));
plot(Time,y_meas,Time(Xpk),Ypk,'dr');
n_peaks=numel(Xpk);  
%Vector length of interest 
t_new = Time(Xpk);  
y_new = Ypk;  
%---Calculate Logarithmic Decrement,  undamped natural frequency, damped natural frequency, damping ratio
Log_Dec = zeros(length(n_peaks));  
for nn = 1:n_peaks-1 %
    Log_Dec(nn)= log(y_new(nn)/y_new(nn+1));  % computed with n = 1 periods apart
end  
Mean_dec = mean(Log_Dec); %Calculate Average Logarithmic Decrement  
%Damping  
damp_ratio_logdec = 1/sqrt(1+((2*pi/(Mean_dec))^2)); %Assesses Damping Constant 

[Star Strider]

Try this —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/871810/Measurement%20Data%20_2.xlsx', 'VariableNamingRule','preserve');

vn = T1.Properties.VariableNames;                                                   % Variable Names
ttlstr = vn([2 3]);                                                                 % Variable Names Used In The Code
x = T1.Time  * 1E-3;                                                                % Change 'Time' From Milliseconds To Seconds
yc(:,1) = T1.(ttlstr{1});
yc(:,1) = smoothdata(yc(:,1), 'sgolay', 250);                                       % Remove Some Noise To Make The Fit Easier
yc(:,2) = T1.(ttlstr{2});

for k = 1:size(yc,2)
    yk = yc(:,k);
    y = yk;
    % y = detrend(y);
    ym = mean(y);                                                                   % Estimate offset
    y = y - ym;
    yu = max(y);
    yl = min(y);
    yr = (yu-yl);                                                                   % Range of ‘y’
    yz = y-yu+(yr/2);
    zci = @(v) find(diff(sign(v)));                                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector
    zt = x(zci(y));
    per = 2*mean(diff(zt));                                                         % Estimate period

    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) - yk);                                                 % Least-Squares cost function
    s0 = [yr; -1E-2;  per;  -1;  ym];                                               % Initial Parameter Estimates
    [s,nmrs] = fminsearch(fcn, s0)                                                  % Minimise Least-Squares

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

    figure
    plot(x,yk,':b', 'LineWidth',1.5)
    hold on
    plot(xp,fit(s,xp), '--r')
    hold off
    grid
    xlabel('Time (s)')
    ylabel('Amplitude')
    title(ttlstr{k})
    legend('Original Data',  'Fitted Curve')
    text(0.1*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.3f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3f%+.3f) %+.4f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
end

s = 5×1

0.1789 -0.1048 3.1084 -1.2119 -0.0062

nmrs = 6.6842

s = 5×1

436.6429 -0.3349 3.5279 -2.0160 274.6755

nmrs = 6.5964e+03

dt = mean(diff(x));

for k = 1:size(yc,2)
    yk = yc(:,k);
    [pks{k},locs{k}] = findpeaks(yk, 'MinPeakProminence',0.1*max(yk-mean(yk)), 'MinPeakDistance',750);
    ipkdist{k} = diff(locs{k});
    ipktime{k} = ipkdist{k} * dt;
    dipktime = mean(ipktime{k})
    figure
    plot(x,yk,':b', 'LineWidth',1.5)
    hold on
    plot(x(locs{k}),yk(locs{k}), '^r', 'MarkerFaceColor','r')
    hold off
    grid
    xlabel('Time (s)')
    ylabel('Amplitude')
    title(ttlstr{k})
    ipktc = compose('\\Delta t = %.3f s\n\\downarrow', ipktime{k});
    text(x(locs{k}(1:numel(locs{k})-1))+dipktime/2, pks{k}(1:numel(locs{k})-1), ipktc, 'Horiz','left','Vert','middle')
end

dipktime = 3.0478

dipktime = 3.0300

The time difference calculations may be robust to other data sets, hwoever since I have only used it with this one, I cannot be certain of that. The code might require some tweaking to work with the others.

.

[Hossam Amin]

Thank you. I appreciate it.

[Star Strider]

As always, my pleasure!

[Hossam Amin]

I was wondering if the fit could be better though? Even for the very first one, because I had a closer analysis onto it and there is phase shift near the end and the amplitude is not nearly fit.

The problem with my data set is that the decay rate and the logarithmic decrement is not constant. (as shown in my code I put earlier). The damping fit function you have used tries to have an approximate by keeping the decay rate and lograithmic decrement constant.

To fix that, is it possible to use two sinus function with the exponential ? I haven't found any similar example on that though. I would be glad if you have the time to help me on this regard. Appreciate it.

[Star Strider]

There are enough data so that the number of parameters estimated from it would not be a problem, however getting a nonlinear regression optimiisation routine to converge on it could be. (For that, one of the Global Optimization Toolbox functions. specifically ga, would likely be required.)

Extending it could go something like this —

 fit = @(b,x)  b(1) .* exp(b(2).*x) .* (sin(2*pi*x./b(3) + 2*pi/b(4))) + exp(b(5).*x) .* (sin(2*pi*x./b(6) + 2*pi/b(7))) + b(8);
 fcn = @(b) norm(fit(b,x) - yk);

with ‘fcn’ (the fitness function or objective function) being the same with the Global Optimization Toolbox as it is here.

Note — There are now 8 parameters, so this is likely beyond the ability of fminsearch (that can only estimate about 7 parameters reliably), and even lsqcurvefit could find it challenging.

.

[Hossam Amin]

OK nice I will try working around with all that information to get the best fit I want.

Lastly, sorry to keep asking so many questions, it's a thing I have to work on since it's not my expertise (data analysis).. :D

How can come up with a code that shows the fit decay rate like this image and also prints the equation of the fit on the plot.

Thanks in advance.

[Star Strider]

My original code prints the equation of the fit on the plot, as does this version.

Adding a plot of the exponential decay envelope simply requires a few changes to the original code —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/871810/Measurement%20Data%20_2.xlsx', 'VariableNamingRule','preserve');
vn = T1.Properties.VariableNames;                                                   % Variable Names
ttlstr = vn([2 3]);                                                                 % Variable Names Used In The Code
x = T1.Time  * 1E-3;                                                                % Change 'Time' From Milliseconds To Seconds
yc(:,1) = T1.(ttlstr{1});
yc(:,1) = smoothdata(yc(:,1), 'sgolay', 250);                                       % Remove Some Noise To Make The Fit Easier
yc(:,2) = T1.(ttlstr{2});
for k = 1:size(yc,2)
    yk = yc(:,k);
    y = yk;
    % y = detrend(y);
    ym = mean(y);                                                                   % Estimate offset
    y = y - ym;
    yu = max(y);
    yl = min(y);
    yr = (yu-yl);                                                                   % Range of ‘y’
    yz = y-yu+(yr/2);
    zci = @(v) find(diff(sign(v)));                                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector
    zt = x(zci(y));
    per = 2*mean(diff(zt));                                                         % Estimate period
    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
    expenv = @(b,x) [-1 1]*b(1) .* exp(b(2).*x) + b(5);                             % Exponential Envelope Decay Function
    fcn = @(b) norm(fit(b,x) - yk);                                                 % Least-Squares cost function
    s0 = [yr; -1E-2;  per;  -1;  ym];                                               % Initial Parameter Estimates
    [s,nmrs] = fminsearch(fcn, s0)                                                  % Minimise Least-Squares
    xp = linspace(min(x),max(x), 500);
    figure
    plot(x,yk,':b', 'LineWidth',1)
    hold on
    plot(xp,fit(s,xp), '--r')
    plot(x, expenv(s,x), '-', 'Color',[1 0.5 0], 'LineWidth',2)
    hold off
    grid
    xlabel('Time (s)')
    ylabel('Amplitude')
    title(ttlstr{k})
    legend('Original Data',  'Fitted Curve', 'Exponential Envelope', 'Location','best')
    text(0.1*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.3f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3f%+.3f) %+.4f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
end

s = 5×1

0.1789 -0.1048 3.1084 -1.2119 -0.0062

nmrs = 6.6842

s = 5×1

436.6429 -0.3349 3.5279 -2.0160 274.6755

nmrs = 6.5964e+03

Make appropriate changes to get the desired result.

.

[Hossam Amin]

And if I want the exponential decay fit on the real measurement in which the curve tries to touch the peaks of the oscillations of the original data, instead of the touching the fitted damping curve itself.?

[Hossam Amin]

... because what I will try to do next is try to estimate the error between both decay curves.

[Star Strider]

Use the envelope function to determine it (findpeaks is also an option, however the noise would make that more difficult), and then use an appropriate nonlinear parameter estimation method to estimate its parameters.

Use a slightly edited version of ‘expenv’ for that:

expenv = @(b,x) b(1) .* exp(b(2).*x) + b(3);                                    % Exponential Envelope Decay Function

Then use it similarly to the way I originally used ‘fit’ to estimate its parameters.

.

[Hossam Amin]

All right. Sounds good.

Is there a rule of thumb in chossing the "Initial Parameter Estimates" s0 = [yr; -1E-2; per; -1; ym]?

Because it requires a lot of trial and error for the different data measurements I have, what's the idea behind trying to choose correct values intuitively to get the fit faster?

[Star Strider]

The ‘s0’ vector elements were chosen to attempt to provide a bit of accuracy with respect to the previous calculations.

Beginning with that many random numbers will rarely result in a decent fit to the data, since there are simply too many values they could take, and too many local minima the function could fall into. For the extended version of ‘fit’ the ‘s0’ vector should be similar for the additional elements, however not the same. I would consider multiplying the initial frequency estimate by 2 (equivalent to dividing the period by 2) and perhaps using random values for the rest of the additional estimates, since the frequency of the second sin function will drive the fit of the associated parameters.

The preferred method however would be to use one of the Global Optimization Toolbox functions since they will intelligently attempt to fit the function with a large number of different parameters, and (given enough time and perhaps a few re-starts) will almost always produce an acceptable solution in a well-posed problem.

[Hossam Amin]

I see, well I don't have license for the Global Optimization Toolbox and the trial is only 30 Days. I guess this is sufficient for my need for now.

I will now subtract the fit curve from the signal in order to calculate the error.

And in order to do that, the best way is name the fit(s0,xp) with certain variable let's say Y_fit. I changed the xp to have the same length as the original data, so xp now = linspace(min(x),max(x), length(yk)). And when I did so, it didn't affect the fit, so we are good. And so now:

Sig_diff = yk - Y_fit and I believe it checks out.

Any insights or maybe a better method you propose ?

[Star Strider]

A simple genetic algorithm is straightforward (although not trivial) to code, and although it will lack the sophistication of the ga function, it should be good enough to get the correct result.

If your research involves a significant amount of optimisations, especially difficult optimisations, the Global Optimization Toolbox is worth getting. If the objective is to create control system models using the input-output relationships of the system, the System Identification Toolbox is definitely worth getting.

There are several ways to calculate the goodness-of-fit. I am not certain where we are just now, however if the design is to compare those two vectors, I would use the fitlm function, with one being the independent and one being the dependent variable. That will produce all the statistics necessary to evaluate how well they match, and additional calls to the ‘mdl’ result will produce additional information, including parameter confidence intervals, and confidence intervals on the regression itself. An alternative is the regress function that will also provide statistics on the fit, however fitlm will likely be better for this problem. (The fitnlm function is also an option, however that requires developing a nonlinear model for the regression.)

These functions are better than simply getting the mean of the differences or even the mean of the squared differences, and doing only simple statistics on the result.

[Hossam Amin]

Great. Thank you for all your help and information. :D

I think I am gonna go with the fitlm or fitnlm for the time being, as it's not the focus of my research to dive in complex optimization algorithms. But it's worth the time and knowledge in the future.

[Star Strider]

As always, my pleasure!

Just to be clear, I was considering fitnlm both with respect to the model regression (use the fminsearch results as the initial parameter estimates in order to get the statistics on the parameters and the fit) as well as comparing the envelope vectors, noting that it will then be necessary to provide a model of some sort to model the envelope functions against each other (most likely a simple exponential or something similar).

I would only consider fitlm to compare the two envelope vectors. They are likely similar enough that a linear hypothesis would work for testing them.

[Hossam Amin]

Aha, I see. You mean I can use fitnlm instead of fminsearch?.. I mean I looked into it and nlinfit and fitnlm are very similar to one another, I used [beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X,Y,modelfun,beta0) instead of the fminsearch ... I set the following:

• initial conditions of s0 as beta0,
• modelfun is the 'fit' which the exponential sinus damping function
• x and y are time and measurement signal
• the output beta is the variables of damping; b(1), b(2) ,,etc

The curve fit slightly changed but it is still a very good fit. Plus now I have the different other outputs in regards to the Covariance and MSE, etc.

Isn't that what you meant by considering fitnlm/nlinfit with respect to the model regression?

And that fitlm can be used for the envelope using a Model which has an exponential only function. However, I don't think I need to compare the two envelopes. By intuition I believe it is just the mirror of one another.

[Star Strider]

‘Isn't that what you meant by considering fitnlm/nlinfit with respect to the model regression?’

Yes! My original code uses fminsearch because everyone has it, and it works for this problem. The ‘extended’ version of the ‘fit’ function would not be able to use fminsearch because it has more parameters than fminsearch can comfortably estimate.

‘I need to compare the two envelopes. By intuition I believe it is just the mirror of one another.’

I thought the objective of that part of the research was to compare one part of the envelope (for example the positive one) as estimated by the regression with the (positive) signal envelope provided by the envelope function. That’s the reason I suggested a regression function that can produce a number of relevant statistics.

.

[Hossam Amin]

Right! .. so just to be clear when you say extended version of 'fit' function you mean the fitnlm/nlinfit and not the ModelFun which in your code is called 'fit' but has no relevance to the in-build matlab 'fit' function or the extended ones which requires the Statistics and Machine Learning Toolbox.

Secondly, on the other part, that might also be interesting to look at. But can you please show me an extended code example for it in regards to the envelope. Thank you.

[Star Strider]

I was in the middle of providing a response to this when Windows 11 crashed. Again. Taking everything I was working on with it. This also occurred with Windows 10, and I was hoping that it would be solve in Windows 11, but of course not. The instability remains. I’m an Instrument-Rated Private Piliot and am truly grateful that Micro$oft doesn’t produce airplanes!

The ‘extended version’ is what I referred to in a previous Comment.

Comparing the envelopes —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/871810/Measurement%20Data%20_2.xlsx', 'VariableNamingRule','preserve');
Nrows = size(T1,1)

Nrows = 31000

vn = T1.Properties.VariableNames;                                                   % Variable Names
ttlstr = vn([2 3]);                                                                 % Variable Names Used In The Code
x = T1.Time  * 1E-3;                                                                % Change 'Time' From Milliseconds To Seconds
yc(:,1) = T1.(ttlstr{1});
yc(:,1) = smoothdata(yc(:,1), 'sgolay', 250);                                       % Remove Some Noise To Make The Fit Easier
yc(:,2) = T1.(ttlstr{2});
for k = 1:size(yc,2)
    yk = yc(:,k);
    [ehi, elo] = envelope(yk, Nrows*0.2, 'peak');
    y = yk;
    % y = detrend(y);
    ym = mean(y);                                                                   % Estimate offset
    y = y - ym;
    yu = max(y);
    yl = min(y);
    yr = (yu-yl);                                                                   % Range of ‘y’
    yz = y-yu+(yr/2);
    zci = @(v) find(diff(sign(v)));                                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector
    zt = x(zci(y));
    per = 2*mean(diff(zt));                                                         % Estimate period
    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
    expenv = @(b,x) [-1 1]*b(1) .* exp(b(2).*x) + b(5);                             % Exponential Envelope Decay Function
    fcn = @(b) norm(fit(b,x) - yk);                                                 % Least-Squares cost function
    s0 = [yr; -1E-2;  per;  -1;  ym];                                               % Initial Parameter Estimates
    [s,nmrs] = fminsearch(fcn, s0)                                                  % Minimise Least-Squares
    xp = linspace(min(x),max(x), 500);
    figure
    hp{1} = plot(x,yk,':b', 'LineWidth',1, 'DisplayName','Original Data');
    hold on
    hp{2} = plot(xp,fit(s,xp), '--r', 'Linewidth',1.5, 'DisplayName','Fitted Curve');
    hp{3} = plot(x, expenv(s,x), '-', 'Color',[1 0.5 0], 'LineWidth',2, 'DisplayName','Exponential Regression Envelope');
    hp{4} = plot(x, [ehi elo], '-g', 'LineWidth',1.5, 'DisplayName','Data Envelope');
    hold off
    grid
    xlabel('Time (s)')
    ylabel('Amplitude')
    title(ttlstr{k})
    legend([hp{1};hp{2};hp{3}(1);hp{4}(1)], 'Location','best')
%     legend('Original Data',  'Fitted Curve', 'Exponential Envelope', 'Location','best')
    text(0.1*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.3f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3f%+.3f) %+.4f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
    
    expenvm = expenv(s,x);
    envr = @(b,x) b(1).*exp(b(2).*x) + b(3);
    B0 = [100; 0.01; mean(expenvm(:,2))];
    envmdl{k} = fitnlm(ehi, expenvm(:,2), envr, B0);
    [ernew,erci] = predict(envmdl{k}, ehi);
    
    envmdl{k}
    
    figure
    plot(ehi, expenvm(:,2), ':b', 'LineWidth',2)
    hold on
    plot(ehi, [ernew erci], '-r', 'LineWidth',1.25)
    hold off
    grid
    xlabel('Data Envelope')
    ylabel('Exponential Regression Envelope')
    title(ttlstr{k})
    legend('Data', 'Regression', '95% Regression Confidence Intervals', 'Location','best')




end

s = 5×1

0.1789 -0.1048 3.1084 -1.2119 -0.0062

nmrs = 6.6842

ans = 
Nonlinear regression model:
    y ~ b1*exp(b2*x) + b3

Estimated Coefficients:
          Estimate        SE        tStat      pValue 
          _________    ________    ________    _______

    b1       50.022      110.92     0.45099      0.652
    b2    0.0085013    0.018827     0.45155    0.65159
    b3      -50.011      110.92    -0.45088    0.65208


Number of observations: 31000, Error degrees of freedom: 30997
Root Mean Squared Error: 0.00673
R-Squared: 0.979,  Adjusted R-Squared 0.979
F-statistic vs. constant model: 7.09e+05, p-value = 0

s = 5×1

436.6429 -0.3349 3.5279 -2.0160 274.6755

nmrs = 6.5964e+03

ans = 
Nonlinear regression model:
    y ~ b1*exp(b2*x) + b3

Estimated Coefficients:
          Estimate        SE        tStat     pValue
          ________    __________    ______    ______

    b1      1.2361     0.0039735    311.09      0   
    b2    0.010609    5.8863e-06    1802.3      0   
    b3      249.22       0.05228    4766.9      0   


Number of observations: 31000, Error degrees of freedom: 30997
Root Mean Squared Error: 3.09
R-Squared: 0.999,  Adjusted R-Squared 0.999
F-statistic vs. constant model: 1.2e+07, p-value = 0

This illustrates one approach. Use whatever works in the context of your project.

At least this time Windows 11 managed to complete this without crashing!

.

[Hossam Amin]

Oh! Sorry to hear that. Yeah I guess they should have no buisness with airplanes LOL.

OK so know I got you all right. Because it is true that the extended ModelFun or 'fit' which include two sinus functions till variable B(8) didn't work with the fminsearch, so I would need to use nlinfit or fitnlm with it "I still need to try it though".

And thanks for sharing the code and for your time.

[Star Strider]

As always, my pleasure!

[Sebastian]

Hello Star Strider,

i´ve also tried to fit my messuered Data with your code. Unfortunately i do not receive any fitted curve. You write that the fit function should be varied. How do you proceed? Is it experience or can you give me an advice? I´ve uploaded my data as well.

Thank you in advance