Fit a standard exponential fit to approximate data

1. For every exponential shaped signal, I want to fit a standard exponential decay fit (I've drawn in black) by joining the possible points in the green line. Also I need the count of such fits in the plot along with the duration of each.

2. I want to find out the slope (I've drawn in violet colour with a question mark) by drawing a line (I've drawn in dark cyan) that connects the starting and ending point of the standard fit. How to do?

Herewith I've attached the Excel sheet and the code that I've tried so far.

clc; close all; clear all;
a=readtable('Data.xlsx');
time=a{2:40000,1};amplitude=a{2:40000,3};
max_peak_values = [];
max_peak_locs = [];
window_size = 1000; 
for i = 1:window_size:length(amplitude)
    end_index = min(i + window_size - 1, length(amplitude));
    segment = amplitude(i:end_index); 
    [max_peak, loc] = max(segment);
    max_peak_values = [max_peak_values; max_peak]; 
    max_peak_locs = [max_peak_locs; loc + i - 1]; 
end
joined_signal = zeros(size(amplitude)); 
joined_signal(max_peak_locs) = max_peak_values; 
for j = 1:length(max_peak_locs)-1
    x = [max_peak_locs(j), max_peak_locs(j+1)];
    y = [max_peak_values(j), max_peak_values(j+1)];
    interp_x = x(1):x(2);
    interp_y = linspace(y(1), y(2), length(interp_x));
    joined_signal(interp_x) = interp_y; 
end
figure;
plot(time, amplitude, 'b'); 
hold on; grid on
plot(time(max_peak_locs), max_peak_values, 'ro'); 
plot(time, joined_signal, 'g'); 

[Image Analyst]

Can you explain how the green line is determined? It looks like you sort of want to fit only the peaks of your data rather than all of the data, however that assumption breaks down in a few places, like around x=0.065 where the green line does not hit a peak, and around 0.04 where it's skipping some peaks.

Can we assume that if there is a run of 0.025 where the values are less than 2, that defines a region that divides total signal into separate exponentially shaped signals? Finding those zones is trivial with thresholding and bwareaopen (to extract only zones of a certain length of longer) and then bwlabel (to give a unique ID number to each flat region.

mask = bwareaopen(signal < 2, minAllowableLength); % minAllowableLength is the minimum number of indexes you'll allow for a seeparation zone.
labeledMask = bwlabel(mask);

What do you need the fits (equation parameters) for? If all you need is the count of signal regions, and the slopes from highest peak to last/final peak, then you don't need to compute a fit to determine just those two things.

And why not fit the exponential curve through the data instead of just the peaks? And if you did that, would you want to exclude those points where the data falls to less than 2, so that you just include the data in the 4-peak regions, not the baseline data?

[Reji G]

The green line is obtained by joining the maximum points for every 1000 samples.

I want to fit a standard exponential fit and find the deviation between the standard fit and the green line. Slope is for further classification.

If I consider all the points then it'd take a longer time. Hence I took maximum points on 1000 samples.

Answer 1

[Sam Chak]

Hi @Reji G

You can use the approach described by @William Rose. Alternatively, you can use the envelope() function. The following example is a code snippet in which I employed the bi-exponential exp2 model instead of the monoexponential model.

data= importdata('Data.xlsx');
t   = data.data(:,1); 
y   = data.data(:,2:3);
yy  = y(:,2);
[mv, idx] = max(yy);

t   = t(idx:20001);
yy  = yy(idx:20001);

[up, ~] = envelope(yy, 2000, 'peak');
hold on
% plot(t, up, 'linewidth', 1.5)

f   = fit(t, up, 'exp2')

f = 
     General model Exp2:
     f(x) = a*exp(b*x) + c*exp(d*x)
     Coefficients (with 95% confidence bounds):
       a =        1544  (1526, 1562)
       b =      -105.7  (-106, -105.3)
       c =       2.438  (2.381, 2.494)
       d =       4.713  (4.451, 4.974)

plot(f, t, yy), grid on
title('Y2'); xlabel('t')

[William Rose]

@Sam Chak, nice!

[William Rose]

@Reji G, @Sam Chak, it is interesting that the fitted line is going up slightly at the end. The fitted values for c and d are positive, so there's a growing exponential in the fit, which is small in this time range.

[Reji G]

Thanks. But I could't understand the code. my data has 2 exponential like signals. What should I modify to fit exponential fits if my data contains multiple number of exponential like signals?

[Sam Chak]

Hi @William Rose,

Your observation is sharp. If the OP expects a monotonically decreasing envelope, then some data processing will be required before fitting the curve to the data.

Hi @Reji G

You should split the signal into the left half (first exponential pattern) and the right half (second exponential pattern). You may manually select the desired local maximum points that describe the envelope of the signal. I learned the 'envelope()' approach from @Star Strider's solutions in other threads. You can refer to past questions to find similar problems.

[William Rose]

@Reji G,

@Sam Chak fits the envelope() of the signal, above. Which is a smart idea, in my opinion. But the resulting fit includes a small rising exponential. If you find the maximum points with findpeaks(), and fit them with 'exp2', you get two decaying exponentials, which is more like what you want, I think. I don't know exactly why you get different results with envelope() versus the maxes from findpeaks(). envelope() seems more elegant but more opaque. With findpeaks(), I always have to experiment a bit with the options, to find the peaks I want. See my comment to my answer, below.

Answer 2

[prabhat kumar sharma]

Hello Reji,

To fit a standard exponential decay to your data and calculate the slope of the line connecting the start and end points of the fit, you can follow these steps in MATLAB. Since I don't have access to your actual data file (Data.xlsx), I'll provide you with a general approach using your existing code as a base.

1. Identify Regions for Exponential Fit: Use the joined_signal to identify segments that resemble exponential decays.
2. Fit Exponential Decay: For each identified segment, fit an exponential decay model.
3. Calculate Slope: Compute the slope of the line connecting the start and end points of the fitted exponential curve.

clc; close all; clear all;
a = readtable('Data.xlsx');
time = a{2:40000, 1};
amplitude = a{2:40000, 3};

% Parameters
window_size = 1000; 
decay_fits = [];

% Find peaks
max_peak_values = [];
max_peak_locs = [];
for i = 1:window_size:length(amplitude)
    end_index = min(i + window_size - 1, length(amplitude));
    segment = amplitude(i:end_index); 
    [max_peak, loc] = max(segment);
    max_peak_values = [max_peak_values; max_peak]; 
    max_peak_locs = [max_peak_locs; loc + i - 1]; 
end

% Create joined signal
joined_signal = zeros(size(amplitude)); 
joined_signal(max_peak_locs) = max_peak_values; 
for j = 1:length(max_peak_locs)-1
    x = [max_peak_locs(j), max_peak_locs(j+1)];
    y = [max_peak_values(j), max_peak_values(j+1)];
    interp_x = x(1):x(2);
    interp_y = linspace(y(1), y(2), length(interp_x));
    joined_signal(interp_x) = interp_y; 
end

% Plot original and joined signals
figure;
plot(time, amplitude, 'b'); 
hold on; grid on;
plot(time(max_peak_locs), max_peak_values, 'ro'); 
plot(time, joined_signal, 'g');

% Fit exponential decay to each segment
for k = 1:length(max_peak_locs)-1
    start_idx = max_peak_locs(k);
    end_idx = max_peak_locs(k+1);

    % Extract segment
    segment_time = time(start_idx:end_idx);
    segment_amplitude = joined_signal(start_idx:end_idx);

    % Fit exponential decay: y = a * exp(b * x)
    f = fit(segment_time, segment_amplitude, 'exp1');
    decay_fits = [decay_fits; f];

    % Plot fitted curve
    plot(f, segment_time, segment_amplitude);

    % Calculate and plot slope
    start_point = [segment_time(1), f.a * exp(f.b * segment_time(1))];
    end_point = [segment_time(end), f.a * exp(f.b * segment_time(end))];
    plot([start_point(1), end_point(1)], [start_point(2), end_point(2)], 'c', 'LineWidth', 1.5);

    % Calculate slope
    slope = (end_point(2) - start_point(2)) / (end_point(1) - start_point(1));
    disp(['Slope for segment ', num2str(k), ': ', num2str(slope)]);
end

hold off;

I hope it helps!

[Reji G]

I've attached Data.xlsx as an excel sheet. Its downloadable. Also the above code gives error

Answer 3

[William Rose]

@Reji G,

[Fix spelling error: fit1(x,y,'exp1') should be fit(x,y,'exp1'...).]

Thank you for providing Data.xlsx.

data=importdata('Data.xlsx');
x=data.data(:,1); y=data.data(:,2:3);
figure; 
subplot(211),plot(t,y(:,1),'-r'); 
grid on; title('Y1')
subplot(212),plot(t,y(:,2),'-b'); 
grid on; title('Y2'); xlabel('X')

The data y2 is the data that matches your drawing.Your code indicates that you are trying to find peaks and interpolate between the peaks. Please review findpeaks() and interp1(). They have many nice options to help you find the peaks you want, and interpolate. Then use fit(). to fit the data. You can use fittype 'exp1' to fit a monoexponential. fit(x,y,'exp1') fits the equation

Y = a*exp(b*x)

Note that the equation above decays to the asymptote x=0. If you want to include a non-zero asymptote in your fit, then you'll need to make a custom fit function, which you can do, see the help for fit().

You may want to offset the x-values: define x1=x-x0, where x0 is the time of the big peak. Then fit y versus x1. You don't have to do this. If you don't do this, you may be surprised when you get a large value for the fit estimate of a, in the equation above. That will happen because a is the fitted y-value at x=0, and if you extend the fitted exponential back to the original x=0, it will be rather large. A lot larger than the peak at x=0.04. Your choice.

For the slope (violet line in your drawing), you can compute that by the standard equation once you find the start and end points you want with findpeaks().

By the way, your data has quantization issues, as the plot below shows. It is because the data in the Excel file is rounded to one digit after the decimal point. I would invesitgate that, and try to get the data with more significant digits.

figure; plot(t,x(:,2),'-b'); 
xlabel('Time'); grid on; title('X2')
xlim([0.094,0.099]); ylim([2,4])

Good luck with your analysis.

[William Rose]

@Reji G,

data=importdata('Data.xlsx');
x=data.data(:,1); y=data.data(:,3);
%[pks,locs]=findpeaks(y,x);  % finds too many peaks
[pks1,locs1]=findpeaks(y,x,MinPeakProminence=1);
% pks1 includes peaks on the rising shoulder which I don;t want.
[pks2,locs2]=findpeaks(y,x,MinPeakProminence=1,MinPeakDistance=0.001);
% Adding MinPeakDistance gets rid of the unwanted peaks on the rising
% shoulders.See plot in top panel, below.
% Plot data and peaks
figure, plot(x,y,'-b',locs1,pks1,'gx',locs2,pks2,'ko'); 
grid on; title('Y'); xlabel('X')
legend('data','pks1','pks2')

% Extract the peaks and locs with locs<0.1, and move the origin 
% for locs to the first peak, as discussed in my answer above.
pks2=pks2(locs2<0.1);               % do this before adjusting locs2
locs2=locs2(locs2<0.1)-locs2(1);    % adjust locs2
% Fit the locs2,pks2 points with a monoexponential.
f1=fit(locs2,pks2,'exp1')

f1 = 
     General model Exp1:
     f1(x) = a*exp(b*x)
     Coefficients (with 95% confidence bounds):
       a =       18.05  (11.86, 24.25)
       b =      -52.45  (-81.96, -22.95)

% Plot the data and fit
figure 
plot(locs2(locs2<0.1),pks2(locs2<0.1),'ko')
hold on
plot(f1,'r')
grid on; xlabel('X'); ylabel('Y'); title('Fit to Peaks')

fit1 (see below) does not looks great. This must be why @Sam Chak suggested a biexponential. Let's try it.

f2=fit(locs2,pks2,'exp2')

f2 = 
     General model Exp2:
     f2(x) = a*exp(b*x) + c*exp(d*x)
     Coefficients (with 95% confidence bounds):
       a =       23.11  (21.61, 24.61)
       b =       -9698  (-1.478e+12, 1.478e+12)
       c =       8.488  (7.603, 9.373)
       d =       -17.8  (-21.38, -14.23)

You said originally that you also wanted slope from first to last peak (violet).

slopeEst=(pks2(end)-pks2(1))/(locs2(end)-locs2(1));

Plot fit2.

plot(f2,'g')
plot(locs2([1,end]),pks2([1,end]),Color=[.5,0,1],LineWidth=2);
legend('peaks','exp1','exp2',['slope=',num2str(slopeEst)]); 
hold off

The fitted values a=23.11, b=-9698 indicate that first exponential in fit2 is esentially a delta function with amplitude 23.11. The time constant is so large and negative that you can only say it is a "really fast decay". The 95% CI for b, for fit2, is non-sensical; a more reasonable CI for b is [-Inf,-5000].