Index exceeds array bounds despite a loop to prevent this?

function [peak_dat_avg] = FindMuscStrength8chan_Cfs(wavedata,channel,stim_freq,stim_time,lat1,lat2)

    artefact_dat = wavedata(:,9,:);
    emg_dat = wavedata(:,channel,:); 

    nframes = size(wavedata,3);
    npulse = single(stim_freq*stim_time);

    emgpeak_dat = zeros(npulse,1,nframes);
    peak_vals = zeros(npulse,1);

    for k = 1:nframes 

        [~, peak_locs] = findpeaks(artefact_dat(:,:,k),'NPeaks',npulse,'MinPeakProminence',0.025,'MaxPeakWidth',5,'MinPeakDistance',700);
        
        start_idx = round(peak_locs + lat1);
        end_idx = round(peak_locs + lat2);

        numb_peaks = numel(peak_locs);

        for i = 1:numb_peaks
            for n = 1:numb_peaks
                if (start_idx(n) > 6000)
                    start_idx(n) = 6000;
                    end_idx(n) = 6000;
                end
            end
            peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
        end

        emgpeak_dat(:,:,k) = peak_vals;

    end

    peak_dat_avg = mean(nonzeros(emgpeak_dat,1));

end

This function is designed to extract a small window of EMG data after locating a stimulation artefact on channel 9 of the data. The issue comes on line 28 where the error 'Index in position 1 exceeds array bounds; Index can't exceed 6000' pops up. I understand this as when trying to select the window of emg_dat it is attempting to start from a sample higher than 6000. However, I tried to implement the if loop above to locate any index values greater than the range of the data and set them to the maximum. I would really appreciate help on fixing this issue

[dpb]

W/o a sample of the data to work with, it's hard to tell what it is you're actually trying to operate on and return which isn't clearly defined here in your explanation.

However, in

    for k = 1:nframes 
      [~, peak_locs] = findpeaks(artefact_dat(:,:,k),'NPeaks',npulse,'MinPeakProminence',0.025,'MaxPeakWidth',5,'MinPeakDistance',700);
      start_idx = round(peak_locs + lat1);
      end_idx = round(peak_locs + lat2);
      numb_peaks = numel(peak_locs);
      for i = 1:numb_peaks
        for n = 1:numb_peaks
          if (start_idx(n) > 6000)
            start_idx(n) = 6000;
            end_idx(n) = 6000;
          end
        end
        peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));

you first set the entire value of the variables start_idx, end_ix then set an element of an array using index n, then refer to the i element of the array. None of those indexing expressions are consistent with each other and so you're reference to the indices in the selection statement is not using the bounded values you just tried to create.

It would appear that what you need to do instead is to ensure when you create start_idx, end_ix in each pass over the number of frames would be to ensure the end_idx array is no greater than the length of the signal you've passed into findpeaks.

I don't follow the purpose of having the doubly nested loop of i and n, it would seem only the i loop would be sufficient to handle each peak within each frame, but you would need to have the peak_vals array doubley dimensioned over nframesXmax(numbpeaks) to save each peak by frame or use a cell array to store the peak values by frame in an array for each frame.

    for k = 1:nframes
      trace=artefact_dat(:,:,k);
      N=numel(trace);
      [~, peak_locs] = findpeaks(trace,'NPeaks',npulse,'MinPeakProminence',0.025,'MaxPeakWidth',5,'MinPeakDistance',700);
      start_idx = round(peak_locs + lat1);
      end_idx = round(peak_locs + lat2);
      end_idx=min(end_idx,N);                   % ensure don't run over end
      ...

The above assumes you wouldn't set the value of lat1 such that start_idx would be past the end of the trace, but bounding it similarly would ensure you wouldn't exceed the data length. But, it would appear that if you can't pull a full peak out of the trace, you might want to simply just skip over the last one in a trace in which case testing for index >N and using continue over the loop over number of peaks would just ignore it leaving you with one less artifact. Of course, if lat2 were too large, you could have the case that the starting location for a subsequent artifact/peak could be before the end of the previous; hopefully your screening with findpeaks is eliminating that potential problem.

Attach a data file and somebody is bound to come along and look at it thoroughly...

[Alexander H]

@dpb Thank you very much for the help. The delay values lat1 and lat2 are fixed muscle latencies so they won't change but I've solved the issue by just ensuring both start_idx and end_idx can't exceed the data size as suggested; I've also removed my nested loops which were a silly effort to solve the problem previously.

Answer 1

[Voss]

The problem is that some element of end_idx is greater than 6000 but the corresponding element of start_idx is not greater than 6000. In that case, the if condition is not true (because start_idx(n) <= 6000), so the code setting start_idx(n) and end_idx(n) to 6000 doesn't get executed, which leaves end_idx(n) greater than 6000 and causes the error.

Also, there is no need for two nested for loops, checking all n for each i. That's redundant, but it doesn't cause any problems. Fixing that, your code could be written as:

for i = 1:numb_peaks
    if (start_idx(i) > 6000)
        start_idx(i) = 6000;
        end_idx(i) = 6000;
    end
    peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
end

Now, to fix the problem outlined above, I assume that in such a case, you'd want only end_idx(i) to be set to 6000, rather than changing both end_idx(i) and start_idx(i), since start_idx(i) is ok (not greater than 6000), in which case you can limit start_idx and end_idx independently:

for i = 1:numb_peaks
    if start_idx(i) > 6000
        start_idx(i) = 6000;
    end
    if end_idx(i) > 6000
        end_idx(i) = 6000;
    end
    peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
end

Another way to do that is using the min function:

for i = 1:numb_peaks
    start_idx(i) = min(start_idx(i),6000);
    end_idx(i) = min(end_idx(i),6000);
    peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
end

And either of those approaches can be done to all elements at once, before the loop:

start_idx(start_idx > 6000) = 6000;
end_idx(end_idx > 6000) = 6000;
for i = 1:numb_peaks
    peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
end

and:

start_idx = min(start_idx,6000);
end_idx = min(end_idx,6000);
for i = 1:numb_peaks
    peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
end

Finally, for robustness you should avoid hard-coding a value like 6000 in the code. (Consider what might happen if you run this function on an emg_dat that has fewer or more than 6000 rows.) Instead get the actual limit from the data itself:

max_idx = size(emg_dat,1);

and use max_idx in place of 6000 in the code, e.g.:

max_idx = size(emg_dat,1);
start_idx = min(start_idx,max_idx);
end_idx = min(end_idx,max_idx);
for i = 1:numb_peaks
    peak_vals(i) = peak2peak(emg_dat(start_idx(i):end_idx(i),:,k));
end

[Alexander H]

Thanks for all the advice it's really helpful :)

[Voss]

You're welcome!