[Eeglablist] Problem with batch PSD

[Makoto Miyakoshi]

Hi Hanna,

 > - peak frequencies are routinely outside the respective frequency
bands (for example theta peak frequencies are around 8-10 Hz)

That's not rare actually. The problem is rather in the assumption that peak
power is always present in the desired frequency range. Therefore people
mechanically calculate average values for their defined freq band. It's a
common practice.

>   - I'm not sure how to work around the events where various lengths of
data was rejected - since big "jumps" in the data can introduce artefacts
and mess up the spectra, I'm sure I should avoid these events, but I'm not
sure about the right way to do it

For the description of the problem, see slide 33 of this ppt pdf.
https://sccn.ucsd.edu/mediawiki/images/2/27/On_preprocessing_EEG_data.pdf

I don't know what EEGLAB does, but I'm sure it uses 'boundary' event marker
which indicates the presence of discontinuity to remove the window with
'boundary' when calculating PSD.
You have two ways to go
(1) You prepare data so that your data length is the integer multiple of
the window used in your welch(). Then you remove the segments that contains
the 'boundary' events.
(2) Basically the same, but in the opposite order: You run spectrogram() to
obtain time-frequency conversion first, then remove the time frames that
contain 'boundary' events.

Somehow I have believed that (2) is the way to go, but it's probably just
as tedious anyway, particularly if you want to use 50% overlap (which is
used in the original definition of PSD together with Hamming window; EEGLAB
uses 0% overlap though as default). Maybe you want to write code for (1).
Your code looks good. Just change the input to pwelch(). Here are some
ideas in pseudocode

(1) Use your welch's window length (use EEG.srate i.e. 1 s because the next
process becomes easier) to epoch your data into data = EEG.data(:,
1:windowLength, numResultantEpoch) You need to trim the last part of the
data for this.
(2) Run something like boundaryEventIdx = find(contains({EEG.event.type},
'boundary')); boundaryLatencies = EEG.event(boundaryEventIdx).latency;
boundaryLatencyS = floor(boundaryLatencies/EEG.srate) This gives you
integer second values for the epochs that contains 'boundary'. Note that by
using 1-s window, the epoch index number is the same as the latency in
second here.
(3) data(:,:,boundaryLatencyS) = []
(4) pwelch(data(:,:))

Makoto


On Wed, Jan 31, 2024 at 10:38 AM Hanna Szakács via eeglablist <
eeglablist at sccn.ucsd.edu> wrote:

> Hello Everyone,
>
> I'm trying to batch process EEG files that have been cleaned with Clean
> rawdata. I'm interested in peak frequencies and average power in 6
> different frequency bands. I have two problems:
>
>    - peak frequencies are routinely outside the respective frequency bands
>    (for example theta peak frequencies are around 8-10 Hz)
>    - I'm not sure how to work around the events where various lengths of
>    data was rejected - since big "jumps" in the data can introduce
> artefacts
>    and mess up the spectra, I'm sure I should avoid these events, but I'm
> not
>    sure about the right way to do it
>
> Here are the relevant segments of my code:
>
> % Loop through each EDF file
> for f = 1:numFiles
>     filename = edfFiles(f).name;
>     fprintf('Processing %s (%d of %d)\n', filename, f, numFiles);
>
>     % Load EEG data
>     EEG = pop_biosig(filename);
>
>     % Initialize variables to accumulate band data
>     bandPeakFreqs = zeros(EEG.nbchan, size(freqBands, 1));
>     bandAvgPowers = zeros(EEG.nbchan, size(freqBands, 1));
>
>     % Loop through each electrode/channel
>     for chan = 1:EEG.nbchan
>         % Perform spectral analysis for each channel
>         [psd, freq] = pwelch(EEG.data(chan,:), windowSize, noverlap, nfft,
> EEG.srate);
>
>         % Analyze each frequency band
>         for j = 1:size(freqBands, 1)
>             fmin = freqBands{j, 2};
>             fmax = freqBands{j, 3};
>
>             % Find the indices of the frequencies within the band
>             bandIndices = freq >= fmin & freq <= fmax;
>
>             % Extract the power spectrum for the band
>             bandPower = psd(bandIndices);
>             bandFreq = freq(bandIndices);
>
>             % Calculate peak frequency and average power
>             [peakPower, peakIdx] = max(bandPower);
>             peakFreq = bandFreq(peakIdx);
>             avgPower = mean(bandPower);
>
>             % Accumulate results for each band
>             bandPeakFreqs(chan, j) = peakFreq;
>             bandAvgPowers(chan, j) = avgPower;
>         end
>     end
>
>     % Average the peak frequencies and average powers across all electrodes
>     avgPeakFreqs = mean(bandPeakFreqs, 1);
>     avgAvgPowers = mean(bandAvgPowers, 1);
>
> And the parameters:
> window size 512
> noverlap 256
> nfft 1024
>
> All EEGs have varying length and varying number of electrodes left
> (originally 128), sampling rate is 1000 Hz.
>
> Thank you in advance for your help,
> Hanna
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to
> eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to
> eeglablist-request at sccn.ucsd.edu
>