My EEG code is not reading my edf file properly. How can I alter the code to fit for my edf file structure of my EEG Recordings?

Hi, I am trying to analyze EEG data in an .edf file format through a series of band analysis and breaking down the data into recording periods. My output is producing only one array of the recording period, even though my EEG recordings have 4 channels and should be producing 4 different arrays. Some of the .edf files I input don't even show an array at all, which leads me to believe that the file is only reading one channel. I incorporated edfread in the code, but seems nothing is working. Can anyone guide me as to how should I change the code in order to read my .edf structures? Thank you so much.

Here is the code:

function [] = BandAnalysisAV(fs, cage_letter, hourStart, outputName)

if cage_letter == 'A'

d=dir('*.edf');

disp(d);

end

for k=1:length(d)

%amp_data=command_read_Intan_RHD2000_file(d(k).name);

% Define the segment size (e.g., 1 minute)

[hdr, amp_data] = edfreadAV(d(k).name, 'assignToVariables', false);

%[hFile] = ns_OpenFile(filename);

%[amp_data] = ns_GetAnalogData(hFile, EntityID, StartIndex, IndexCount);

end

% section out 30 minute windows from each of the 4 channels

sample_window = fs*60*30;

% preallocate matrix for storage

v = zeros(length(d), sample_window);

% Determine # of 30 min. recording periods

recPeriods = floor((d(1).bytes/sample_window)/2);

dataStore.delta = zeros(4,length(recPeriods));

dataStore.theta = zeros(4,length(recPeriods));

dataStore.alpha = zeros(4,length(recPeriods));

dataStore.beta = zeros(4,length(recPeriods));

dataStore.gamma = zeros(4,length(recPeriods));

dataStore.highFreq = zeros(4,length(recPeriods));

dataStore.alphadelta = zeros(4,length(recPeriods));

for t = 1:recPeriods

% Update on computing progress

fprintf('Proccessing recording period %d of %d', t, recPeriods);

% 2nd order filter to remove line frequency from the signal,

% which is frequencies between 59 and 61 Hz

filt = designfilt('bandstopiir','FilterOrder',2, ...

'HalfPowerFrequency1',59,'HalfPowerFrequency2',61, ...

'DesignMethod','butter','SampleRate',fs);

% Remove line frequency from each channel using filter

for n = 1:length(d)

fid = fopen(d(n).name);

fseek(fid, (t*sample_window)-sample_window,'bof');

v(n,:) = filter(filt,(fread(fid, sample_window, 'int16')));

fclose(fid);

end

% Preallocate arrays for the skew and root mean square error for every

% second of the recording. For a 30 min. recording period at

% a sampling frequency of 2000 Hz that is 1800 seconds.

skw = zeros(4, 1802); % skw stores the average in 1801 and channel type (good or bad) in 1802

rtms = zeros(4, 1801); % rtms stores the average in 1801

zscore_rtms = zeros(4,1800);

max_zscore = zeros(4,359);

% Find the skew and root mean square error of the 30 minute

% recording period

for i = 1:length(skw(:,1))

for n = 1:length(skw(1,:)) - 2

skw(i,n) = skewness(v(i,(2000*(n-1) + 1:2000*n)));

rtms(i,n) = rms(v(i,((2000*(n-1) + 1:2000*n))));

end

skw(i,1801) = mean(skw(i, 1:end-2));

rtms(i,1801) = mean(rtms(i, 1:end-1));

zscore_rtms(i,:) = zscore(rtms(i,1:end-1));

% Identify bad channels

% Algorithm detects bad channels with 90% accuracy

if rtms(i,1801) > 200 || rtms(i,1801) < 30 || skw(i,1801) > 0.4

skw(i,1802) = -1;

else

skw(i,1802) = 1;

end

end

% Detect the max z-score for 5 second epochs to detect bad epochs

for i=1:length(zscore_rtms(:,1))

a = 1;

for j=1:5:length(zscore_rtms(1,:))-5

max_zscore(i,a)=max(zscore_rtms(i, j:j+5));

a = a + 1;

end

end

% Set bad epochs in the data to 0

for o = 1:length(v(:,1))

for c = 1:359

if abs(max_zscore(o,c)) >= 3

v(o, (c*10000 - 9999):c*10000) = 0;

end

end

end

summary = zeros(4,7);

for b = 1:length(v(:,1))

tempData = nonzeros(v(b,:)'); % remove the bad epochs in the data

tempPower = zeros(convergent(length(tempData)/10000), 7); % preallocate an array to store the power for tempData

% Note: each 30 min recording period uses the same tempData file

% which is written over for each recording period

if skw(b, 1802) == 1 % only perform band analysis on good channels

for c = 1:(length(tempData)/10000) % divide by 12500, equivalent of 5 seconds with fs = 2500 Hz

tempPower(c,:) = bandAnalysis((tempData((c*10000 - 9999):c*10000, 1)), fs);

% see bandAnalysis function below for details

end

summary(b,1) = mean(tempPower(:,1));

summary(b,2) = mean(tempPower(:,2));

summary(b,3) = mean(tempPower(:,3));

summary(b,4) = mean(tempPower(:,4));

summary(b,5) = mean(tempPower(:,5));

summary(b,6) = mean(tempPower(:,6));

summary(b,7) = mean(tempPower(:,7));

else

% Mark bad channels as NaN

summary(b,1) = 0/0;

summary(b,2) = 0/0;

summary(b,3) = 0/0;

summary(b,4) = 0/0;

summary(b,5) = 0/0;

summary(b,6) = 0/0;

summary(b,7) = 0/0;

end

end

dataStore.delta(:,t) = summary(:,1);

dataStore.theta(:,t) = summary(:,2);

dataStore.alpha(:,t) = summary(:,3);

dataStore.beta(:,t) = summary(:,4);

dataStore.gamma(:,t) = summary(:,5);

dataStore.highFreq(:,t) = summary(:,6);

dataStore.alphadelta(:,t)= summary(:,7);

clc;

end

% This is revised code from writeRHS_bands_CSV.m

timeStamps = zeros(1, length(recPeriods));

for a = 1:recPeriods

if hourStart > 24

timeStamps(a) = hourStart-24;

else

timeStamps(a) = hourStart;

end

hourStart = hourStart + 0.5;

end

dataStore.delta=[timeStamps; dataStore.delta];

dataStore.theta=[timeStamps; dataStore.theta];

dataStore.alpha=[timeStamps; dataStore.alpha];

dataStore.beta=[timeStamps; dataStore.beta];

dataStore.gamma=[timeStamps; dataStore.gamma];

dataStore.highFreq=[timeStamps; dataStore.highFreq];

dataStore.alphadelta=[timeStamps; dataStore.alphadelta];

allName=[outputName, '.xlsx'];

writematrix(dataStore.delta,allName, 'Sheet', 'delta');

writematrix(dataStore.theta,allName, 'Sheet', 'theta');

writematrix(dataStore.alpha,allName, 'Sheet', 'alpha');

writematrix(dataStore.beta,allName, 'Sheet', 'beta');

writematrix(dataStore.gamma,allName, 'Sheet', 'gamma');

writematrix(dataStore.highFreq,allName, 'Sheet', 'fastGamma');

writematrix(dataStore.alphadelta,allName, 'Sheet', 'alphadelta');

end

function [power] = bandAnalysis(data, fs)

[butterFilt, coButter] = butter(6, ((100/fs)*2)); % create a 6th order butterworth filter

dataF = filter(butterFilt, coButter, transpose(data)); % filter the data passed to the function

dataF = transpose(dataF);

N = length(dataF);

xdf = real(fft(dataF)); % All this is the fast-fourier transform and

xdfs = xdf(1:((N/2) + 1)); % the adjustments required to get the power

pxx = (1/(fs*N)) * abs(xdfs).^2; % spectral density

pxx(2:end-1) = 2*pxx(2:end-1);

pfreq = 0:fs/N:fs/2; % Create a frequency vector

normPxx = pxx/trapz(pxx); % Normalize the PSD to the passed in data

% Preallocate the 'power' vector

power = zeros(1,7);

% Separate the 6 EEG frequencies

power(1) = trapz(normPxx(pfreq >= 0.1 & pfreq <= 4)); % Delta bands

power(2) = trapz(normPxx(pfreq > 4 & pfreq <= 8)); % Theta bands

power(3) = trapz(normPxx(pfreq > 8 & pfreq <= 13)); % Alpha bands

power(4) = trapz(normPxx(pfreq > 13 & pfreq <= 25)); % Beta Bands

power(5) = trapz(normPxx(pfreq > 25 & pfreq <= 50)); % Gamma Bands

power(6) = trapz(normPxx(pfreq > 50 & pfreq <= 100)); % High-Freq Bands

power(7) = trapz(normPxx(pfreq > 8 & pfreq <= 13)./(trapz(normPxx(pfreq >= 0.1 & pfreq <= 4)))); % Alpha-Delta Ratio

end

Answer 1

[Karan Singh]

Hey Alvin, I had an overview of your code, and from the parts of you describing it to me I have a few suggestions.

1. Check the function that reads the .EDF files: It's important to ensure that the function you are using (“edfreadAV”) correctly reads the .EDF files and assigns the data to the “amp_data” variable. You can test this function separately to make sure it is working as expected.
2. Extract data for each channel: Currently, the code assumes that there are four channels, but it doesn't explicitly extract the data for each channel. To address this, you can modify the code inside the loop to extract the data for each channel using indexing.

For example, instead of:

v(n,:) = filter(filt, (fread(fid, sample_window, 'int16')));

You can use:

v(n,:) = filter(filt, amp_data(n, (t*sample_window)-sample_window+1 : t*sample_window));

• Close the file after use: It's a good practice to close the file after reading its contents. You can add the line “fclose”(fid) to close the file handle and free up system resources.
• Check the dimensions of the “v” matrix: It's a good idea to verify that the “v” matrix has the expected dimensions. In this case, it should have dimensions of 4 x “sample_window”. You can check this by adding a check, such as “size(v)”, to ensure it matches the expected dimensions.

I hope these suggestions help!

[Alvin]

Hello Karan! Thank you so much for the feedback. I incorporated your suggestions to the code, such as making sure the "v" matrix is in the right dimensions by changing the preallocated dimension to 4 and changing the filer for v. But I am now running into a matrix error within the filter code. The error is occuring when the code reaches t=5, where an array bound error occurs since t=5 exceeds the maximum.

ERROR: Index in position 2 exceeds array bounds. Index must not exceed 16292000.

For context, the dimensions of "amp_data" double is 5 x 16292000. I believe the error is occuring because of the (t*sample_window)-sample_window+1 : t*sample_window)) exceeding bounds at t=5, given that the sample_window is equal to 3,600,000. Could you possibly guide me through this? Thank you for all your help--this is my first time writing code for EEG time-frequency analysis, so any help/suggestions would be very useful.

Attached here is the revised (and simplified) version:

function [] = BandAnalysisAV(fs, cage_letter, hourStart, outputName)

if cage_letter == 'A'

d=dir('*.edf');

disp(d);

end

for k=1:length(d)

[hdr, amp_data] = edfreadAV(d(k).name, 'assignToVariables', false);

% section out 30 minute windows from each of the 4 channels

sample_window = fs*60*30;

% preallocate matrix for storage

v = zeros(4, sample_window);

size(v)

% Determine # of 30 min. recording periods

recPeriods = floor((d(1).bytes/sample_window)/2);

dataStore.delta = zeros(4,length(recPeriods));

dataStore.theta = zeros(4,length(recPeriods));

dataStore.alpha = zeros(4,length(recPeriods));

dataStore.beta = zeros(4,length(recPeriods));

dataStore.gamma = zeros(4,length(recPeriods));

dataStore.highFreq = zeros(4,length(recPeriods));

dataStore.alphadelta = zeros(4,length(recPeriods));

for t = 1:recPeriods

% Update on computing progress

fprintf('Proccessing recording period %d of %d', t, recPeriods);

% 2nd order filter to remove line frequency from the signal,

% which is frequencies between 59 and 61 Hz

filt = designfilt('bandstopiir','FilterOrder',2, ...

'HalfPowerFrequency1',59,'HalfPowerFrequency2',61, ...

'DesignMethod','butter','SampleRate',fs);

% Remove line frequency from each channel using filter

for n = 1:length(d)

fid = fopen(d(n).name);

fseek(fid, (t*sample_window)-sample_window,'bof');

v(n,:) = filter(filt, amp_data(n,(t*sample_window)-sample_window+1 : t*sample_window));

fclose(fid);

end