fitting envelope peak points

I am facing some problems trying to create a fitting curve that will pass through maximum peaks. I tried findpeaks

findpeaks(Sig)
hold all
plot(Sig,'k')

as well as envelope

[up, lo] = envelope(Sig,'peak');
plot(Sig, 'o-k');
hold on
plot(up,'r');
plot(lo,'b');
legend ('bb','up','down')
grid

How can i create a fitting curve which will pass through max peak points? 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]

Image Analyst thanks for your answer but i does not work. See the data please. I want something like this:

[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?

[Ancalagon8]

Thanks for your answer again! I want a regressed line that fits the peaks with an exponential decay, like the demo you attached.

[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.

[Ancalagon8]

Image Analyst sorry again, I tried first to use the code:

s = load('sig.mat')
y = Sig3;
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-');
xlim ([0 5000])

and then made some changes to the demo, but i get only a straight line..

format long g;
format compact;
fontSize = 20;
% Define coefficients and function that the X values obey.
a = 10 % 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 = 11 % Arbitrary sample values I picked.
bGuessed = 0.5
cGuessed = 6
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') 

I think the solution is really close!

[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') 

[Ancalagon8]

Thank you for the immediate response and all this help! A result close to the yellow or the green curve would be great.. Maybe a good idea would be to find tha max peak and start the curve from this point, so as to be close enough to envelope..

[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.

[Ancalagon8]

@Image Analyst would it be possible to adjust the above code for power law and not exponential fitting?

[Image Analyst]

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

Answer 2

[Ancalagon8]

@Image Analyst thank you once more for your answer!I tried to change the code so as to receive something like this but i did not work for me:

[Image Analyst]

Wow, that's bizarre, but we've seen it before. See this link.

It seems to be related to the OpenGL renderer.

Or see this answer.