nested loop or indexing in a 3D matrix - EEG/LFP data

Hello everyone,

i am struggling to find a solution to the following problem:

I have a 3d matrix, let's call it TD_32 (5000,32,86). It's a time series, where the first dimension represents the actual data, the second the channels and the third the trials (This data set is a EEG like data structure).

Threshold values have been found and calculated in another matrix (let's call it pre_rms_single_chan_4x) and it varies between all the trials and all the channels. This mean the threshold will end up being a 32x86 matrix. in the following lines of code, the matrix TD_32 have been inverted as the threshold is actually a negative peak.

[val_rms, idx_rms] = findpeaks(-TD_32_denoise_post_single_chan_single_tr, ...
    MinPeakHeight=pre_rms_single_chan_4x);

Now, I need to find all the values in the 3D matrix TD_32 which exceed the threshold, and store them into a new 3D matrix with boolean variables.

          idx_matrix(idx_rms) = true;

Unfortunately, there is a bug in my code which prevents it from working for all the dimensions. The code correctly compute only the values in the first dimension but stop working after the first trial for the first channel.

The error message i get from matlab:

Warning: Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.

> In findpeaks>removePeaksBelowMinPeakHeight (line 504)

In findpeaks (line 146)

I tried to manually compute the threshold and find the exceeding values from several cannels and trials and it works. This mean the issue(s) is in the code itself.

for ch=2:32
    for tr=2:86
          TD_32_denoise_post_single_chan = TD_32_denoise_post(:,ch,:);
          TD_32_denoise_post_single_chan = squeeze(TD_32_denoise_post_single_chan);
          pre_rms_single_chan_4x = pre_rms_all_chan_4x(ch,tr);
          pre_rms_single_chan_4x = squeeze(pre_rms_single_chan_4x);
          TD_32_denoise_post_single_chan_single_tr = TD_32_denoise_post(:,ch,tr);
          [val_rms, idx_rms] = findpeaks(-TD_32_denoise_post_single_chan_single_tr, ...
    MinPeakHeight=pre_rms_single_chan_4x);
          idx_matrix(idx_rms) = true;
    end
end

[Enzo]

@Star Strider @Stephen23 @dpb i summon the best in town. sorry for bothering you, but I am stuck and cannot move on in my work. thanks in advance

[Varun]

Hello! Can you share all the necessary variables in a .mat file? That will be helpful in investigating the issue better!

Answer 1

[Siraj]

Hii! It is my understanding that you have a matrix “TD_32” of dimensions 5000X32X86.

32 represents the number of channels and 86 represents the number of trials.

There also exists a threshold matrix “pre_rms_single_chan_4x” of dimensions 32X86.

The sign of this threshold matrix is inverted that is you now want to find values “smaller” than this threshold, which means that you want to find the values that are greater in magnitude than the threshold but have a negative sign, for each (channel, trial) pair.

“findpeaks(data)” returns a vector with the local maxima (peaks) of the input signal vector, data. A local peak is a data sample that is either larger than its two neighboring samples or is equal to Inf.

On running the attached code, we get the following warning.

Warning: Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.

> In findpeaks>removePeaksBelowMinPeakHeight (line 509)

In findpeaks (line 146)

In nested_loop_3d_matrix (line 16)

In my understanding the above warning is because of passing negative threshold and negative values to the “findpeaks()” function.

“findpeaks(-TD_32_denoise_post_single_chan_single_tr, MinPeakHeight=pre_rms_single_chan_4x);”

For example, let “TD_32_denoise_post_single_chan_single_tr” be [7,7,3,6,1] and the corresponding threshold be 2.

Here we can clearly see that we have some values present which are greater than 2(threshold).

Now “-TD_32_denoise_post_single_chan_single_tr” will be [-7, -7, -3, -6, -1], and the corresponding negative threshold will be -2.

Now we can clearly see that even though the values are larger in magnitude but eventually are smaller than -2.

Therefore, this results in the above warning.

Refer to the following code for better understanding.

rng(0);
TD_32_denoise_post = randi([0, 9], 5, 3, 2);
pre_rms_all_chan_4x = -randi([0 2], 3, 2);
idx_matrix = zeros(5,3,2);

%% Using -TD_32_denoise_post_single_chan_single_tr and threshold.
for ch =   1:3
    for tr = 1:2
        TD_32_denoise_post_single_chan_single_tr  = TD_32_denoise_post(:,ch,tr);
        disp("-TD_32_denoise_post_single_chan_single_tr");
        disp(-TD_32_denoise_post_single_chan_single_tr);

        threshold_single_chan_single_tr = pre_rms_all_chan_4x(ch, tr);
        disp("Threshold");
        disp(threshold_single_chan_single_tr);

        [val_rms, idx_rms] = findpeaks(-TD_32_denoise_post_single_chan_single_tr, "MinPeakHeight",threshold_single_chan_single_tr);
        idx_rms(idx_rms,ch,tr) = 1;


    end
end

-TD_32_denoise_post_single_chan_single_tr

    -8
    -9
    -1
    -9
    -6

Threshold

    -2

-TD_32_denoise_post_single_chan_single_tr

    -1
    -4
    -9
    -7
    -9

Threshold

     0

Warning: Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.

-TD_32_denoise_post_single_chan_single_tr

     0
    -2
    -5
    -9
    -9

Threshold

     0

-TD_32_denoise_post_single_chan_single_tr

    -6
     0
    -8
    -9
    -6

Threshold

     0

Warning: Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.

-TD_32_denoise_post_single_chan_single_tr

    -1
    -9
    -9
    -4
    -8

Threshold

     0

Warning: Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.

-TD_32_denoise_post_single_chan_single_tr

    -7
    -7
    -3
    -6
    -1

Threshold

    -2

Warning: Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.

%% Using TD_32_denoise_post_single_chan_single_tr and -threshold.
for ch =   1:3
    for tr = 1:2
        TD_32_denoise_post_single_chan_single_tr  = TD_32_denoise_post(:,ch,tr);
        disp("TD_32_denoise_post_single_chan_single_tr");
        disp(TD_32_denoise_post_single_chan_single_tr);

        threshold_single_chan_single_tr = pre_rms_all_chan_4x(ch, tr);
        disp("-Threshold");
        disp(-threshold_single_chan_single_tr);

        [val_rms, idx_rms] = findpeaks(TD_32_denoise_post_single_chan_single_tr, "MinPeakHeight",-threshold_single_chan_single_tr);
        idx_rms(idx_rms,ch,tr) = 1;


    end
end

TD_32_denoise_post_single_chan_single_tr

     8
     9
     1
     9
     6

-Threshold

     2

TD_32_denoise_post_single_chan_single_tr

     1
     4
     9
     7
     9

-Threshold

     0

TD_32_denoise_post_single_chan_single_tr

     0
     2
     5
     9
     9

-Threshold

     0

TD_32_denoise_post_single_chan_single_tr

     6
     0
     8
     9
     6

-Threshold

     0

TD_32_denoise_post_single_chan_single_tr

     1
     9
     9
     4
     8

-Threshold

     0

TD_32_denoise_post_single_chan_single_tr

     7
     7
     3
     6
     1

-Threshold

     2

Refer to the following documentation to see which data points qualify as “peaks” is “findpeaks()”.

https://in.mathworks.com/help/signal/ref/findpeaks.html

Hope this helps.