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
Curve Fitting of large Data Measurement?
14 次查看（过去 30 天）
显示 更早的评论
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.
2 个评论
Hossam Amin
2022-1-23
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.
采纳的回答
Star Strider
2022-1-24
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.)
.
26 个评论
Hossam Amin
2022-1-24
编辑：Hossam Amin
2022-1-24
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
2022-1-24
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.
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
2022-1-24
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
2022-1-24
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
2022-1-25
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
2022-1-25
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
2022-1-26
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
2022-1-26
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
2022-1-26
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
2022-1-26
... because what I will try to do next is try to estimate the error between both decay curves.
Star Strider
2022-1-26
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
2022-1-26
编辑：Hossam Amin
2022-1-26
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
2022-1-27
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
2022-1-27
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
2022-1-27
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
2022-1-27
编辑：Hossam Amin
2022-1-27
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
2022-1-28
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
2022-1-28
编辑：Hossam Amin
2022-1-28
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
2022-1-29
‘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
2022-1-29
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
2022-1-29
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!
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
2022-1-30
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.
Sebastian
2023-8-21
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
更多回答（0 个）
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)