Help fixing code to do loop aswell as how to show start/end times of spike as seen on plot.

4 views (last 30 days)
This is what I'm getting at for my script. I am trying to add a loop or redo my script and add a loop for it to show when it shows a spike at a certain pH level, and says the Start time and End time as well as the duration of the spike. I have included the script file and Csv file so if you wish to run it to see what I can do to accomplish what I'm trying to get at. Thanks Please comment if you need more info. Thanks for stopping by!
Please also check the images!

Answers (1)

Mathieu NOE
Mathieu NOE on 8 Dec 2021
hello Erick
I have implemented a new portion of code to define the two crossing points with maximal time accuracy (by linear interpolation)
this is the result
I simply commented the portion of your code that I didn't need for my little demo , but you can easily bring that back alive.
code :
%25 November Erick Delgado EGR 122, Project 2
clearvars
clc
%get the data from table
T = readtable('pH_Data.csv'); %Reads Table or CSV File
T.Properties.VariableNames = {'Date_Time','Time_Run','Volt_Run','pH_Run'}; %Give Labels New names
% Date_and_Time_Run Time_(s)_Run_1 Voltage_(mV)_Run_l_ pH_Run_1
%set min pH limit for Carp fish
MinpH = 7.6;
% plot(T.Date_Time,T.pH_Run) %Give us a plot of pH and Time
% title('pH Data from November 2021');
% xlabel('Time');
% ylabel('pH');
% yline(MinpH)
%% my code starts here
x = T.Date_Time;
y = T.pH_Run;
threshold = MinpH; % your value here
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(y,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
figure(1)
plot(x,y,x,threshold*ones(size(x)),'k--',t0_pos,s0_pos,'dr',t0_neg,s0_neg,'dg','linewidth',2,'markersize',12);grid on
legend('signal','threshold','positive slope crossing points','negative slope crossing points');
title('pH Data from November 2021');
xlabel('Time');
ylabel('pH');
time_delta = t0_pos - t0_neg
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% %get pH data
% pH_data = T.pH_Run;
%
% %determine Whether water in our watershed exceeded any levels of pH
% flag1 = false;
%
% for index = 1:length(pH_data)
% if pH_data(index) < MinpH
% flag1 = true; %Identifies the event occurred
% break
% end
% end
%
% if flag1 == true
% disp('Water in our watershed exceeded min levels of pH')
%
% %determine how many times the limit(s) is exceeded (number of events)
% eventInfo = int16.empty;
% thisEventTime = 0;
% flag2= false;
% indexNow = 0;
% indexAfter = 0;
%
% eventTimes = double.empty(0,2);
% for index = 1:length(pH_data)
% if pH_data(index) < MinpH && flag2 == false %Found that the event started
% indexNow = index; %Get the time when the event started
% flag2 = true; %Mark event start
% elseif pH_data(index) >= MinpH && flag2 == true %Found that the event is over
% indexAfter = index; %Get the time when the event is over
% thisEventTime = round((indexAfter - indexNow) * 3 / 60); %Calculate the duration of this event
% eventTimes = [eventTimes; indexNow indexAfter];
% eventInfo(end + 1) = thisEventTime; %Add the duration of the event
% flag2 = false; %Mark event is over
% end
% end
% T.Date_Time(eventTimes(1):eventTimes(2))
% %display number of events and the total duration of time (in hours) the limit(s) was exceeded for each event
% disp('====================================================');
% disp('The number of elements shows the number of events');
% disp('Each number shows the total duration of time for each event in hours');
% disp(eventInfo);
% else
% disp('Water in our watershed did not exceeded min levels of pH')
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if ~isempty(ind)
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
else
% empty output
ind_pos = [];
t0_pos = [];
s0_pos = [];
% extract the negative slope crossing points
ind_neg = [];
t0_neg = [];
s0_neg = [];
end
end
%
  6 Comments
Mathieu NOE
Mathieu NOE on 13 Dec 2021
hello Erick
can you please clarify what do meaun by "multiple" spikes ? the data provided seems to show only one peak (now identified by the balck diamond)
the code helps you identify three points in your data
so which points are to be identified ? and do have have data with multiple spikes ?
%25 November Erick Delgado EGR 122, Project 2
clearvars
clc
%get the data from table
T = readtable('pH_Data.csv'); %Reads Table or CSV File
T.Properties.VariableNames = {'Date_Time','Time_Run','Volt_Run','pH_Run'}; %Give Labels New names
% Date_and_Time_Run Time_(s)_Run_1 Voltage_(mV)_Run_l_ pH_Run_1
%set min pH limit for Carp fish
MinpH = 7.6;
% plot(T.Date_Time,T.pH_Run) %Give us a plot of pH and Time
% title('pH Data from November 2021');
% xlabel('Time');
% ylabel('pH');
% yline(MinpH)
%% my code here
x=T.Date_Time;
y = T.pH_Run;
threshold = MinpH; % your value here
%% find crossing points
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(y,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
%% find peak point between crossing points
ind = find(x>t0_neg & x<t0_pos);
signal_extract = y(ind);
[val,loc] = min(signal_extract);
figure(1)
plot(x,y,x,threshold*ones(size(x)),'k--',x(ind(loc)),val,'dk',t0_pos,s0_pos,'dr',t0_neg,s0_neg,'dg','linewidth',2,'markersize',12);grid on
legend('signal','threshold','peak','positive slope crossing points','negative slope crossing points');
title('pH Data from November 2021');
xlabel('Time');
ylabel('pH');
time_delta = t0_pos - t0_neg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% %get pH data
% pH_data = T.pH_Run;
%
% %determine Whether water in our watershed exceeded any levels of pH
% flag1 = false;
%
% for index = 1:length(pH_data)
% if pH_data(index) < MinpH
% flag1 = true; %Identifies the event occurred
% break
% end
% end
%
% if flag1 == true
% disp('Water in our watershed exceeded min levels of pH')
%
% %determine how many times the limit(s) is exceeded (number of events)
% eventInfo = int16.empty;
% thisEventTime = 0;
% flag2= false;
% indexNow = 0;
% indexAfter = 0;
%
% eventTimes = double.empty(0,2);
% for index = 1:length(pH_data)
% if pH_data(index) < MinpH && flag2 == false %Found that the event started
% indexNow = index; %Get the time when the event started
% flag2 = true; %Mark event start
% elseif pH_data(index) >= MinpH && flag2 == true %Found that the event is over
% indexAfter = index; %Get the time when the event is over
% thisEventTime = round((indexAfter - indexNow) * 3 / 60); %Calculate the duration of this event
% eventTimes = [eventTimes; indexNow indexAfter];
% eventInfo(end + 1) = thisEventTime; %Add the duration of the event
% flag2 = false; %Mark event is over
% end
% end
% T.Date_Time(eventTimes(1):eventTimes(2))
% %display number of events and the total duration of time (in hours) the limit(s) was exceeded for each event
% disp('====================================================');
% disp('The number of elements shows the number of events');
% disp('Each number shows the total duration of time for each event in hours');
% disp(eventInfo);
% else
% disp('Water in our watershed did not exceeded min levels of pH')
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if ~isempty(ind)
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
else
% empty output
ind_pos = [];
t0_pos = [];
s0_pos = [];
% extract the negative slope crossing points
ind_neg = [];
t0_neg = [];
s0_neg = [];
end
end
%

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!