3-term exponential fit

49 次查看（过去 30 天）
William
William 2014-4-6
回答： Arturo Gonzalez 2020-9-8
Hi, I want to fit a 3-term exponential function i.e.
y = a*exp(-b*x) + c*exp(-d*x) + e*exp(-f*x)
to get coefficients b,d,f.
Matlab's curvefitting toolbox is great for 2 term fitting, but that is it's limit. Tried log fits, polyfit fit but had no luck. Any ideas?
Thanks, Will

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2014-4-6
I don’t have the Curve Fitting Toolbox (Optimization and Statistics instead). This uses fminsearch:
y = @(b,x) b(1).*exp(-b(2).*x) + b(3).*exp(-b(4).*x) + b(5).*exp(-b(6).*x);
p = [3; 5; 7; 11; 13; 17]*1E-1; % Create data
x = linspace(1, 10);
yx = y(p,x) + 0.1*(rand(size(x))-0.5);
OLS = @(b) sum((y(b,x) - yx).^2); % Ordinary Least Squares cost function
opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
B = fminsearch(OLS, rand(6,1), opts); % Use ‘fminsearch’ to minimise the ‘OLS’ function
fcnfit = y(B,x); % Calculate function with estimated parameters
figure(1)
plot(x, yx, '*b')
hold on
plot(x, fcnfit, '-r')
hold off
grid
  3 个评论
显示 1更早的评论
William
William 2014-4-13
编辑：William 2014-4-13
Think I'm getting the hang of it now. I had an extra noise component in my signal that was throwing off the fit. When I take this out of the signal, and do a two term fit, the method works great.
I've attached the noisy data and the output fit I get. The noise results in the oscillation seen in the surf image, and messes up the fit.
The three terms have
b(1)=1/120*10^-3, b(2)=1/80*10^-3,b(3)=1/30*10^-3
Below shows the spread of the fits. Apologies if I haven't explained it very well, I have a feeling the noise ( with b(2)) introduces an unavoidable error.
p.s. I've nulled fit values below zero and above 200
Star Strider
Star Strider 2014-4-13
The best model and fit will actually require you to go back to the differential equations that model the process you’re fitting, integrate them (analytically if possible), and fit that solution. Considering your data demonstrate an exponential-periodic behaviour, chances are that the parameters are not actually independent and in all likelihood are functionally related. I would not consider the oscillations ‘noise’ unless you know they are an artifact of your instrumentation or some such.

请先登录，再进行评论。

更多回答（2 个）

John D'Errico
John D'Errico 2014-4-13
Note that polyfit is meaningless in this context. You can't fit an exponential with a polynomial tool!
Next, fitting sums of exponentials is one form of an ill-posed problem. When you have too many terms in that sum (and 3 terms is starting to be a fair amount) then the problem starts to get nasty. Think of it like this, we are trying to form a linear combination of terms that all look pretty much alike! Estimation of the coefficients will generate singular (or nearly singular) matrices. If the rate parameters of several of the terms are too close, the problem will become difficult to solve using double precision arithmetic.
  1 个评论
John D'Errico
John D'Errico 2014-4-13
I'll try to add an example later when I have a chance. I'll also show how the fit process can be made a bit easier using a partitioned least squares.

请先登录，再进行评论。

Arturo Gonzalez
Arturo Gonzalez 2020-9-8
Per this answer, you can do it with the following matlab code
https://math.stackexchange.com/a/3808325/605948
clear all;
clc;
% get data
dx = 0.001;
x = (dx:dx:1.5)';
y = -1 + 5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);
% calculate n integrals of y and n-1 powers of x
n = 3;
iy = zeros(length(x), n);
xp = zeros(length(x), n+1);
iy(:,1) = cumtrapz(x, y);
xp(:,1) = x;
for ii=2:1:n
iy(:, ii) = cumtrapz(x, iy(:, ii-1));
xp(:, ii) = xp(:, ii-1) .* x;
end
xp(:, n+1) = ones(size(x));
% get exponentials lambdas
Y = [iy, xp];
A = pinv(Y)*y;
Ahat = [A(1:n)'; [eye(n-1), zeros(n-1, 1)]];
lambdas = eig(Ahat);
lambdas
% get exponentials multipliers
X = [ones(size(x)), exp(lambdas'.*x)];
P = pinv(X)*y;
P
% show estimate
y_est = X*P;
figure();
plot(x, y); hold on;
plot(x, y_est, 'r--');

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by