Splitting Graph into Sections based on Peaks Found

Hi I have the following looking data and I am trying to do is to first extract the peaks and troughs which I have done as in the images. Once this is done i need to be able to create plots from the peak to the trough to create both expodental growth and expodential decay curves. I have been trying to automate this so that it can be done with many different plots and each plot is stored and then curve fitting can be done on each curve. I cant seem to get it to plot from one peak to the next trough without checking it every tiem and hard coding it in. If you could help that would be great.

Figure 1: Raw Data

Figure 2: Peaks Found

Figure 3: Troughs Found

Here is the code I am working on below. The part I am trying to automate is the (% Extract data from cell up till index section) and the time section after that. The time section will be identical to the extract data section it just comes from the time vector. The data file is attached as well so you can run the code.

% Extract Time and Cell 2
data = xlsread('Data.xlsx','Sheet1','A1:o327');

% Extract Time 
t = data(:,1);

% Extract Cell 2
c2 = data(:,3);

%Plot Cell 2 Data
figure(1)
plot(t,c2);

% Find Peaks and locations   
[peakValues, indexesOfPeaks] = findpeaks(c2,t, 'MinPeakHeight',2e4 , 'MinPeakDistance', 10);

% Find Troughs and Locations
[troughValues, indexesOfTrough] = findpeaks(-c2,t, 'MinPeakDistance', 10);

% Plot the peaks
figure(3)
findpeaks(c2,t, 'MinPeakHeight',2e4 , 'MinPeakDistance', 10);

% Plot the troughs
figure(4)
findpeaks(-c2,t, 'MinPeakDistance', 40);

% Make Trough Values positive
postiveTrough = abs(troughValues);

% Find and store all instances of index of Peak in data of cell
iPM = zeros(length(peakValues),1);
for i = 1:length(peakValues)
   
    indexPeakMatrix = find(c2 == peakValues(i));
    iPM(i+1) = indexPeakMatrix;
    
end

% Find and store all instances of index of Trough in data of cell
iTM = zeros(length(postiveTrough),1);
for i = 1:length(postiveTrough)
   
    indexTroughMatrix = find(c2 == postiveTrough(i));
    iTM(i+1) = indexTroughMatrix;
    
end

% Extract data from cell up till index 
p1 = c2(1:iPM(2));
dataArray{1} = p1;
% fprintf('Section 1: iPM(2) = %d\n',iPM(2)); 
p2 = c2(iPM(2):iTM(3));
dataArray{2} = p2;
% fprintf('Section 2: iPM(2) = %d and iTM(3) = %d\n',iPM(2), iTM(3)); 
p3 = c2(iTM(3):iPM(3));
dataArray{3} = p3;
% fprintf('Section 3: iTM(3) = %d and iPM(3) = %d\n',iTM(3), iPM(3)); 
p4 = c2(iPM(3):iTM(5));
dataArray{4} = p4;
% fprintf('Section 4: iPM(3) = %d and iTM(5) = %d\n',iPM(3), iTM(5)); 
p5 = c2(iTM(5):iPM(4));
dataArray{5} = p5;
% fprintf('Section 5: iTM(5) = %d and iPM(4) = %d\n',iTM(5), iPM(4));
p6 = c2(iPM(4):iTM(8));
dataArray{6} = p6;
% fprintf('Section 6: iPM(4) = %d and iTM(8) = %d\n',iPM(4), iTM(8));
p7 = c2(iTM(8):iPM(5));
dataArray{7} = p7;
% fprintf('Section 7: iTM(8) = %d and iPM(5) = %d\n',iTM(8), iPM(5));
p8 = c2(iPM(5):iTM(12));
dataArray{8} = p8;
% fprintf('Section 8: iPM(5) = %d and iTM(11) = %d\n',iPM(5), iTM(12));
p9 = c2(iTM(12):iPM(6));
dataArray{9} = p9;
% fprintf('Section 9: iTM(12) = %d and iPM(6) = %d\n',iTM(12), iPM(6));
p10 = c2(iPM(6):iTM(15));
dataArray{10} = p10;
% fprintf('Section 10: iPM(6) = %d and iTM(15) = %d\n',iPM(6), iTM(15));
p11 = c2(iTM(15):iPM(7));
dataArray{11} = p11;
% fprintf('Section 11: iTM(15) = %d and iPM(7) = %d\n',iTM(15), iPM(7));
p12 = c2(iPM(7):iTM(18));
dataArray{12} = p12;
% fprintf('Section 12: iPM(7) = %d and iTM(18) = %d\n',iPM(7), iTM(18));
p13 = c2(iTM(18):iPM(8));
dataArray{13} = p13;
% fprintf('Section 13: iTM(18) = %d and iPM(8) = %d\n',iTM(18), iPM(8));
p14 = c2(iPM(8):iTM(20));
dataArray{14} = p14;
% fprintf('Section 14: iPM(8) = %d and iTM(20) = %d\n',iPM(8), iTM(20));
p15 = c2(iTM(20):iPM(9));
dataArray{15} = p15;
% fprintf('Section 15: iTM(20) = %d and iPM(9) = %d\n',iTM(20), iPM(9));
p16 = c2(iPM(9):iTM(23));
dataArray{16} = p16;
% fprintf('Section 16: iPM(9) = %d and iTM(23) = %d\n',iPM(9), iTM(23));
p17 = c2(iTM(23):iPM(10));
dataArray{17} = p17;
% fprintf('Section 17: iTM(23) = %d and iPM(10) = %d\n',iTM(23), iPM(10));
p18 = c2(iPM(10):iTM(25));
dataArray{18} = p18;
% fprintf('Section 18: iPM(10) = %d and iTM(25) = %d\n',iPM(10), iTM(25));
p19 = c2(iTM(25):iPM(11));
dataArray{19} = p19;
% fprintf('Section 19: iTM(25) = %d and iPM(11) = %d\n',iTM(25), iPM(11));
p20 = c2(iPM(11):iTM(28));
dataArray{20} = p20;
% fprintf('Section 20: iPM(11) = %d and iTM(28) = %d\n',iPM(11), iTM(28));
p21 = c2(iTM(28):iPM(12));
dataArray{21} = p21;
% fprintf('Section 21: iTM(27) = %d and iPM(12) = %d\n',iTM(27), iPM(12));
p22 = c2(iPM(12):iTM(30));
dataArray{22} = p22;
% fprintf('Section 22: iPM(12) = %d and iTM(30) = %d\n',iPM(12), iTM(30));

% Extract data from time up till index
t1 = t(1:iPM(2));
timeArray{1} = t1;
t2 = t(iPM(2):iTM(3));
timeArray{2} = t2;
t3 = t(iTM(3):iPM(3));
timeArray{3} = t3;
t4 = t(iPM(3):iTM(5)); 
timeArray{4} = t4;
t5 = t(iTM(5):iPM(4));
timeArray{5} = t5;
t6 = t(iPM(4):iTM(8));
timeArray{6} = t6;
t7 = t(iTM(8):iPM(5));
timeArray{7} = t7;
t8 = t(iPM(5):iTM(12));
timeArray{8} = t8;
t9 = t(iTM(12):iPM(6));
timeArray{9} = t9;
t10 = t(iPM(6):iTM(15));
timeArray{10} = t10;
t11 = t(iTM(15):iPM(7));
timeArray{11} = t11;
t12 = t(iPM(7):iTM(18));
timeArray{12} = t12;
t13 = t(iTM(18):iPM(8));
timeArray{13} = t13;
t14 = t(iPM(8):iTM(20));
timeArray{14} = t14;
t15 = t(iTM(20):iPM(9));
timeArray{15} = t15;
t16 = t(iPM(9):iTM(23));
timeArray{16} = t16;
t17 = t(iTM(23):iPM(10));
timeArray{17} = t17;
t18 = t(iPM(10):iTM(25));
timeArray{18} = t18;
t19 = t(iTM(25):iPM(11));
timeArray{19} = t19;
t20 = t(iPM(11):iTM(28));
timeArray{20} = t20;
t21 = t(iTM(28):iPM(12));
timeArray{21} = t21;
t22 = t(iPM(12):iTM(30));
timeArray{22} = t22;


% Curve Fit all data Sections
% for i = 1:22
% 
%     % Curve fot each section
%     createFit(timeArray{i},dataArray{i});
% 
% end

Answer 1

[Mathieu NOE]

hello David

quick and dirty trial to make your code work automatically and as good as it can in my limited time availability

there is for sure room for further improvements

the "b" coefficients of the exponential solutions are stored in b_up and b_down

see code below :

clc
clearvars

% Extract Time and Cell 2
data = xlsread('Data.xlsx','Sheet1','A1:o327');
% Extract Time 
t = data(:,1);
dt = mean(diff(t)); % avearge time interval between peaks

% Extract Cell 2
c2 = data(:,3);
%Plot Cell 2 Data
figure(1)
plot(t,c2);
% Find Peaks and locations   
[peakValues, TimeOfPeaks] = findpeaks(c2,t, 'MinPeakHeight',2e4 , 'MinPeakDistance', 10);
% Plot the peaks
figure(2)
findpeaks(c2,t, 'MinPeakHeight',2e4 , 'MinPeakDistance', 10);

period = mean(diff(TimeOfPeaks)); % average time interval between peaks
buffer_before_peak = round(0.25*period); % seconds
buffer_after_peak = round(0.7*period); % seconds

for ci = 1:length(TimeOfPeaks)
    
    %% before peak exp fit 
    t_start = TimeOfPeaks(ci) - buffer_before_peak;
    t_stop = TimeOfPeaks(ci) - 0*dt; % stop time is equal to peak time
    % my test code 
    ind = find(t>=t_start & t<=t_stop);
    x_up = t(ind);
    x_up = x_up - x_up(1); % time starts at zero always to simplify the exp fit 
    y2fit_up = c2(ind);

    [sol_up,y_fit_up] =createFit_up(x_up,y2fit_up);
    
    % store b coefficient 
    b_up(ci) = sol_up(2);

    
     %% after peak exp fit (decay)
    t_start = TimeOfPeaks(ci) +2*dt ; % start time is 2 samples after peak
    t_stop = TimeOfPeaks(ci) + buffer_after_peak;
    % my test code 
    ind = find(t>=t_start & t<=t_stop);
    x_down = t(ind);
    x_down = x_down - x_down(1); % time starts at zero always to simplify the exp fit  
    y2fit_down = c2(ind);

    [sol_down,y_fit_down] =createFit_decay(x_down,y2fit_down);
    % store b coefficient 
    b_down(ci) = sol_down(2);

    figure(ci+2),
    subplot(211),plot(x_up,y2fit_up,'r',x_up,y_fit_up,'-.k');
    title(['Positive slope fit / peak # ' num2str(ci) ]);
    xlabel('seconds');
    legend('data','exp fit');   
    subplot(212),plot(x_down,y2fit_down,'r',x_down,y_fit_down,'-.k');
    title(['Negative slope fit / peak # ' num2str(ci) ]);
    xlabel('seconds');
    legend('data','exp fit');   
    
end
    

function [sol,y_fit] =createFit_up(x,y2fit)
    % exponential decay fit
    % code is giving good results with template equation : % y = a.*exp(b.*(x-c)) + d;

    f = @(a,b,c,d,x) a.*exp(b.*(x-c)) + d;
    obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y2fit);
    sol = fminsearch(obj_fun, [y2fit(end)-y2fit(1),0,0,y2fit(1)]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    d_sol = sol(4);
    y_fit = f(a_sol, b_sol,c_sol ,d_sol, x);

end

function [sol,y_fit] =createFit_decay(x,y2fit)
    % exponential decay fit
    % code is giving good results with template equation : % y = a.*(1-exp(b.*(x-c))) + d;

    f = @(a,b,c,d,x) a.*(1-exp(b.*(x-c))) + d;
    obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y2fit);
    sol = fminsearch(obj_fun, [y2fit(end)-y2fit(1),0,0,y2fit(1)]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    d_sol = sol(4);
    y_fit = f(a_sol, b_sol,c_sol ,d_sol, x);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

the 11 last plots show in the upper subplot the "rising" portion of the signal before each peak, the lower subplot is for the decay side

to make my life easier , each time a portion of signal is considered , its corresponding time vector always start at zero - you can change that yourself if you prefer

here two examples (first and last peak)

[David Cowan]

Hi Mathieu,

This works really well and is what I had been struggling with thank you very much. The next part I have done on the old code that needs to be moved to this new setup is the ablity to get the overall r-squared value of the growth curve and the decay curves. I have been able to do this with the following setup using the fit function and saving the goodness to a cell array and extracting the r-squared value at each iteration then calcauting the mean value. Can this be implemented in the new code? Also I would aslo like the standard devation but i have not been able to find a solution to that.

Thanks

David

% Curve Fit all growth sections
 for i = 1:2:22
     
     % Curve Fitting
     [curve, goodness_growth{i}, output_growth{i}] = fit(timeArray{i},dataArray{i},'exp1');
     
     % Plot each Curve
     figure(4)
     plot(curve,timeArray{i},dataArray{i});
     
     hold on;
     
     % Print out each R-Square Fit Value
     fprintf('R-Squared Value at Curve %d = %d\n', i, goodness_growth{1, i}.rsquare);
     
     % Obtain the mean and std of r-sqaure fit values
     mean_growth = mean(goodness_growth{1, i}.rsquare);

 end
 
 % print overall R-Square Value fit
 fprintf('Mean R-Squared Value Overall Growth = %d\n\n', mean_growth);

[Mathieu NOE]

hello David

I added the R2 computation in my code , it also appears in the legend box FYI

I created a small subfunction for that as I don't have the curve fitting toolbox

you can probably work around how you want to save or use your fprintf lines of code to display it.

hope it helps

clc
clearvars
% close all

% Extract Time and Cell 2
data = xlsread('Data.xlsx','Sheet1','A1:o327');
% Extract Time 
t = data(:,1);
dt = mean(diff(t)); % avearge time interval between peaks

% Extract Cell 2
c2 = data(:,3);
%Plot Cell 2 Data
figure(1)
plot(t,c2);
% Find Peaks and locations   
[peakValues, TimeOfPeaks] = findpeaks(c2,t, 'MinPeakHeight',2e4 , 'MinPeakDistance', 10);
% Plot the peaks
figure(2)
findpeaks(c2,t, 'MinPeakHeight',2e4 , 'MinPeakDistance', 10);

period = mean(diff(TimeOfPeaks)); % averarge time interval between peaks
buffer_before_peak = round(0.25*period); % seconds
buffer_after_peak = round(0.7*period); % seconds

for ci = 1:length(TimeOfPeaks)
    
    %% before peak exp fit 
    t_start = TimeOfPeaks(ci) - buffer_before_peak;
    t_stop = TimeOfPeaks(ci) - 0*dt; % stop time is equal to peak time
    % my test code 
    ind = find(t>=t_start & t<=t_stop);
    x_up = t(ind);
    x_up = x_up - x_up(1); % time starts at zero always to simplify the exp fit 
    y2fit_up = c2(ind);

    [sol_up,y_fit_up] =createFit_up(x_up,y2fit_up);
    
    % store b coefficient 
    b_up(ci) = sol_up(2);
    
    % Rsquared coefficient
    Rsquared_up(ci) = my_Rsquared_coeff(y2fit_up,y_fit_up);

    
     %% after peak exp fit (decay)
    t_start = TimeOfPeaks(ci) +2*dt ; % start time is 2 samples after peak
    t_stop = TimeOfPeaks(ci) + buffer_after_peak;
    % my test code 
    ind = find(t>=t_start & t<=t_stop);
    x_down = t(ind);
    x_down = x_down - x_down(1); % time starts at zero always to simplify the exp fit  
    y2fit_down = c2(ind);

    [sol_down,y_fit_down] =createFit_decay(x_down,y2fit_down);
    
    % store b coefficient 
    b_down(ci) = sol_down(2);
    
    % Rsquared coefficient
    Rsquared_down(ci) = my_Rsquared_coeff(y2fit_down,y_fit_down);

    figure(ci+2),
    subplot(211),plot(x_up,y2fit_up,'r',x_up,y_fit_up,'-.k');
    title(['Positive slope fit / peak # ' num2str(ci) ]);
    xlabel('seconds');
    legend('data',['exp fit / R2 = ' num2str(Rsquared_up(ci)) ]);   
    subplot(212),plot(x_down,y2fit_down,'r',x_down,y_fit_down,'-.k');
    title(['Negative slope fit / peak # ' num2str(ci) ]);
    xlabel('seconds');
    legend('data',['exp fit / R2 = ' num2str(Rsquared_down(ci)) ]);   
    
end
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [sol,y_fit] =createFit_up(x,y2fit)
    % exponential decay fit
    % code is giving good results with template equation : % y = a.*exp(b.*(x-c)) + d;

    f = @(a,b,c,d,x) a.*exp(b.*(x-c)) + d;
    obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y2fit);
    sol = fminsearch(obj_fun, [y2fit(end)-y2fit(1),0,0,y2fit(1)]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    d_sol = sol(4);
    y_fit = f(a_sol, b_sol,c_sol ,d_sol, x);

end

function [sol,y_fit] =createFit_decay(x,y2fit)
    % exponential decay fit
    % code is giving good results with template equation : % y = a.*(1-exp(b.*(x-c))) + d;

    f = @(a,b,c,d,x) a.*(1-exp(b.*(x-c))) + d;
    obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y2fit);
    sol = fminsearch(obj_fun, [y2fit(end)-y2fit(1),0,0,y2fit(1)]);
    a_sol = sol(1);
    b_sol = sol(2);
    c_sol = sol(3);
    d_sol = sol(4);
    y_fit = f(a_sol, b_sol,c_sol ,d_sol, x);

end

function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Mathieu NOE]

hello David

if my answer has fullfilled your expectations, do you mind accepting it ?

thanks

[Mathieu NOE]

hello David

if my answer has fullfilled your expectations, do you mind accepting it ?

thanks