# How can I store data from an if loop into a 1D vector

5 views (last 30 days)
INDS on 10 Jan 2020
Answered: Meg Noah on 12 Jan 2020
function [Peak] = findpeaks(t,ECG)
% threshold = "user to insert threshold"
for k = 2 : length(ECG) - 1
if (ECG(k) > ECG(k-1) && ECG(k) > ECG(k+1) && ECG(k) > threshold)
disp('prominant peak found')
ECG(k)
k
plot (t, ECG, t(k), ECG(k), '.')
hold on
end
end
I currently have this as my code to find peaks of data given by an ECG. The data given is a 1D vector and I would like to save the output of this function as a 1D vector.
I was also wondering whether it could be done through the use of the diff function.
Last thing, I'm not sure how to get the user to insert a threshold. Thank you!
Guillaume on 10 Jan 2020
Yes, it can be done with diff. With and without diff you certainly don't need a loop.
You haven't said what the output of the function should be. A vector of indices (i.e. the k at which the peaks are found), maybe?
I'm not sure how to get the user to insert a threshold
From the point of view of the function, it should be an input just as t and ECG. Don't do user interaction inside a processing function. From the point of the view of the caller of the function, there are many ways to get that threshold from very basic (input) to slightly more fancy (inputdlg for example) to full-fledged GUI. Again, it shouldn't be the job of findpeaks to communicate with the user.

rtransform on 10 Jan 2020
Edited: rtransform on 10 Jan 2020
Dear INDS,
One way would be to use a counter-variable within the if-condition which increases every time the condition is met.
Then this counter can be used to index a storage vector.
e.g
function [peak] = findpeaks(t,ECG)
% threshold = "user to insert threshold"
counter = 1;
storageVec = zeros(length(ECG),2) % max possible length, you can extract only the non-zero values later on
% you would not even need to define this vector here, then you can also omit "peak = peak(1:counter);"
% but in general i would not recommend to have a variable change size in every iteration of a loop
prompt = 'Threshold: ';
threshold = input(prompt); % user needs to enter nuber and hit enter
% this loop is unnecessary, though. use logicals instead!
for k = 2 : length(ECG) - 1
if (ECG(k) > ECG(k-1) && ECG(k) > ECG(k+1) && ECG(k) > threshold)
disp('prominant peak found')
ECG(k)
k
peak(counter,1) = ECG(k);
peak(counter,2) = t(k);
counter = counter+1;
plot (t, ECG, t(k), ECG(k), '.')
hold on
end
end
peak = peak(1:counter); % removes the zero values at the end
end
I would highly reccomend, though, to check out MATLAB's built-in findpeaks function, which has a lot of options like minpeakdistance, minpeakprominence, threshold,...
This makes it really easy to look for R-peaks, T-waves,...
Considering using the diff-function, you could use it on the ecg signal to find maxima, but without additional steps such as thresholding, this gives you all maxima, not just the R-peaks.
For user input, see code above.
Edit:
Just like Guillaume pointed out, for something like a threshold query you don't need a loop, though.
try something like
function [peak, times] = findpeaks(t, ecg, threshold)
peak = ecg(ecg > threshold); % note these are not the indices, but the values of the R-peaks
times = t(ecg > threshold);
end
Kind regards
##### 2 CommentsShowHide 1 older comment
rtransform on 10 Jan 2020
my bad, of course storageVec and peak should be the same. (I updated the code above)
But I can only stress again, if you want your code to be fast, don't use a loop in this case.

Meg Noah on 12 Jan 2020
This function returns the indices where the peaks for ECG values in time-segments defined by a user-specified threshold
function [idxPeakVal] = findpeaksInECG(ECG,threshold)
% Step 1
% First, determine points of the ECG scan exceed a given threshold.
idx = find(ECG > threshold);
% Step 2
% Next, determine which of the data-points found in the previous step
% correspond to the same heartbeat. (Hint: data-points corresponding
% to one heartbeat will have consecutive indices, whereas transitions
% between heartbeats correspond to discontinuities in the index. diff
% function might help.
idxStartNewHeartBeat = [idx(1) idx(find(diff(idx)>1)+1)']';
% Step 3
% For each of the groups of data-points (corresponding to different
% heartbeats), find the location of the element with the maximum amplitude,
% and store this location in a vector.
stats = regionprops(labeledHeartbeats,ECG,'MaxIntensity'); %#ok<MRPBW>
maxAmplitudeValues = [stats.MaxIntensity];
idxPeakVal = zeros(size(maxAmplitudeValues));
for iheartbeat = 1:length(maxAmplitudeValues)
idxPeakVal(iheartbeat) = find(labeledHeartbeats == iheartbeat & ...
ECG == maxAmplitudeValues(iheartbeat));
end
end
To run the program:
nx = numel(ECG); x = 1:nx;
Fs = 350; % [samples/s] frequency
T = 1/Fs; % [s] sampling period
N = 3000; % [samples] Length of signal
t = (0:N-1)*T; % [s] Time vector
% visualization
figure('color','white','position',[70 100 900 500]);
plot(1e3*t,ECG,'k','DisplayName','ECG Amplitude')
title({'Heartbeat data'})
xlabel('t (milliseconds)')
ylabel('ECG(t)')
hold on
threshold = 100;
[idxPeakVal] = findpeaksInECG(ECG,threshold);
% plot the peak values;
plot(1e3*t(idxPeakVal),ECG(idxPeakVal),'or','MarkerFaceColor','r');
This must be some sort of classroom activity - there are a number of questions that are all sounding about the same.