MATLAB Answers

Splitting Graph into Sections based on Peaks Found

조회 수: 1(최근 30일)
David Cowan
David Cowan 2021년 8월 24일
댓글: Mathieu NOE 2021년 9월 9일
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

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 8월 24일
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)
  댓글 수: 4
Mathieu NOE
Mathieu NOE 2021년 9월 9일
hello David
if my answer has fullfilled your expectations, do you mind accepting it ?
thanks

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by