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

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]

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.

Answer 1

[rtransform]

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

[INDS]

This code seems to work in terms of plotting the peaks on the graph. However, I get this error for the following line

Undefined function or variable 'peak'.

Error in findpeaks (line 24)

peak = peak(1:counter);

function [peak] = findpeaks(ECG,threshold)
 
    counter = 1;
    storageVec = zeros(length(ECG),2)   
    
 
f_s = 350;
N = length(ECG);
t = [0:N-1]/f_s;

    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
            storageVec(counter,1) = ECG(k);
            storageVec(counter,2) = t(k);
            counter = counter+1; 
            plot (t, ECG, t(k), ECG(k), '.')
            hold on
        end
    end
   peak = peak(1:counter); 
   
end

[rtransform]

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.

Answer 2

[Meg Noah]

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. 
mask = zeros(size(ECG)); mask(idx) = 1;
labeledHeartbeats = bwlabel(mask);
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:

% load the data
ECG = dlmread('ECG.csv');
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.