fitting peak points maxima

I am facing some problems trying to create a curve that will pass through maximum peaks. Will a solution be to export peaks into a new array, plot them and then create a fitting curve or there is another way?

Thanks in advance!

Answer 1

[Image Analyst]

Just try interp1() or spline(). Pass in your peak points and then the x values of the whole array, something like

[peakValues, indexes] = findpeaks(Sig);
yOut = spline(indexes, peakValues, 1:length(Sig));
plot(yOut, 'b-', 'LineWidth', 2);
grid on;

Let me know if it works. Attach your Sig data if it doesn't. Accept this answer if it does.

[Ancalagon8]

Have you seen it?

[Image Analyst]

Well that data looks totally different than the first set, but it shouldn't matter. But, no I didn't see your code - you forgot to post it so I can't fix it. All I saw was your mat file, not the code - where is the code that you tried attached? Was it something like this:

s = load('sig.mat')
y = s.Sig3b;
x = 1 : length(y);
plot(x, y, 'b.-', 'MarkerSize', 20);
grid on
[peakValues, peakIndexes] = findpeaks(y, 'MinPeakDistance', 15);
hold on;
plot(peakIndexes, peakValues, 'r*');
% Interpolate
yFit = spline(peakIndexes, peakValues, x);
plot(x, yFit, 'c-');

Do you want something that exactly hits the top of each peak (like the above code), or a regressed line that fits the peaks with an exponential decay, like the attached demo?

[Image Analyst]

OK. Try to adapt it and show your results. The demo is well commented so a smart engineer like you should be able to follow it.

[Image Analyst]

It helps to manually solve the equation using a single point identified on the fit to get the estimated coefficients. But it doesn't look like it's a pure exponential decay, or else we need to do a better job with findpeaks() to get the right peaks.

clc;    % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;

s = load('sig.mat')
y = s.Sig3b;
% Get rid of negative.
y(y<0) = [];
x = 1 : length(y);
plot(x, y, 'b', 'MarkerSize', 15);
grid on
[peakValues, peakIndexes] = findpeaks(y, 'MinPeakDistance', 15);
hold on;
plot(peakIndexes, peakValues, 'r.');
% Interpolate
yFit = spline(peakIndexes, peakValues, x);
plot(x, yFit, 'r-');

% Define coefficients and function that the X values obey.
a = 100 % Arbitrary sample values I picked.
b = 0.4
c = 5
tbl = table(x(:), y(:));
% Define the model as Y = a * exp(-b*x) + c
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * exp(-b(2)*x(:, 1)) + b(3);  

% Guess values to start with.  Just make your best guess.
aGuessed = 8000 % Arbitrary sample values I picked.
bGuessed = 0.005
cGuessed = 0
beta0 = [aGuessed, bGuessed, cGuessed]; % Guess values to start with.  Just make your best guess.

% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);

coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * exp(-coefficients(2)*x) + coefficients(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(x, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 30;

% Place formula text roughly in the middle of the plot.
formulaString = sprintf('Y = %.3f * exp(-%.3f * X) + %.3f', coefficients(1), coefficients(2), coefficients(3))
xl = xlim;
yl = ylim;
xt = xl(1) + abs(xl(2)-xl(1)) * 0.325;
yt = yl(1) + abs(yl(2)-yl(1)) * 0.59;
text(xt, yt, formulaString, 'FontSize', 25, 'FontWeight', 'bold');

set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 

[Image Analyst]

We need to make the table from ONLY the peak points, not all of them, so

% Define coefficients and function that the X values obey.
tbl = table(x(peakIndexes)', y(peakIndexes)');

NOT

% Define coefficients and function that the X values obey.
tbl = table(x(:), y(:));

New code:

clc;    % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;

s = load('sig.mat')
y = s.Sig3b;
% Get rid of negative.
y(y<0) = [];
x = 1 : length(y);
plot(x, y, 'b', 'MarkerSize', 15);
grid on
[peakValues, peakIndexes] = findpeaks(y, 'MinPeakDistance', 15);
hold on;
plot(peakIndexes, peakValues, 'r.', 'MarkerSize', 30);
% Interpolate
yFit = spline(peakIndexes, peakValues, x);
plot(x, yFit, 'r-');

% Define coefficients and function that the X values obey.
tbl = table(x(peakIndexes)', y(peakIndexes)');
% Define the model as Y = a * exp(-b*x) + c
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * exp(-b(2)*x(:, 1)) + b(3);  

% Guess values to start with.  Just make your best guess.
aGuessed = 6000 % Arbitrary sample values I picked.
bGuessed = 0.004
cGuessed = 0
beta0 = [aGuessed, bGuessed, cGuessed]; % Guess values to start with.  Just make your best guess.

% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);

coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * exp(-coefficients(2)*x) + coefficients(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(x, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 30;

% Place formula text roughly in the middle of the plot.
formulaString = sprintf('Y = %.3f * exp(-%.3f * X) + %.3f', coefficients(1), coefficients(2), coefficients(3))
xl = xlim;
yl = ylim;
xt = xl(1) + abs(xl(2)-xl(1)) * 0.325;
yt = yl(1) + abs(yl(2)-yl(1)) * 0.59;
text(xt, yt, formulaString, 'FontSize', 25, 'FontWeight', 'bold');

set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 

Note that there are some peaks early on below 3000, so you'll need to experiment around with findpeaks() if you want to eliminate those and get a better fit.

Note: If you want to specify a point for the y axis intercept, that is possible by adjusting the model. As it is, it tries to give its best guess.

[Image Analyst]

Yes. Did you try to change the model? It's not hard. See attached.