[Eeglablist] Problem with batch PSD

Hanna Szakács szakacsmiriamhanna at gmail.com
Wed Jan 31 00:13:12 PST 2024


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


More information about the eeglablist mailing list