I am processing an EEG brain signal, which has up to 64 data channels sampled at 500 Hz. One of the analyses consists of extracting the ratio of alpha/delta power, where alpha represents the waveforms which frequencies are between 8 and 13 Hz, and delta the waveforms which frequencies are between 1 and 4 Hz.
The signal is first filtered with a 1-15 Hz BP filter (
scipy Butterworth 4th order, output 'sos'). Additionally, a common average projector and EOG SSP artifact correction projector are applied. Those are not relevant to this question.
My program compute the alpha/delta ratio at 2 different moments:
- A short window (500 ms or 1 second, i.e. 250 or 500 samples)
- A long window (4 to 10 seconds)
However, EEG signal is very susceptible to noise, i.e. if the participant moves his eyes, jaws, head, .. the brain signal will be completely masked by artifacts. Thus, a rejection criteria is also applied. It's a simple peak-to-peak rejection, i.e. if the
max - min is larger than the criteria on at least one of the channels, this window (epoch) is rejected.
For the short window, I can directly implement this as follows:
if any(np.ptp(data, axis=0) > reject['eeg']): continue # skip
channel x samples array (e.g.
reject['eeg'] the threshold value.
However, for the large window, I can not disregard the complete window. My idea was to cut down the e.g. 5 seconds signal into 1 second chunks, and apply the rejection criteria on the chunk. Something like:
data = raw.get_data(picks='eeg') # np.ndarray 64x2500 for 5 seconds data = data.reshape(64, 5, 500) # I hope this is the correct reshape, not tested. chunk2remove = np.any(np.ptp(data, axis=2) > reject['eeg'], axis=0) data = data[:, np.where(chunk2remove == False), :]
However, now I do not know if I can safely apply FFT on the signal. I do not know what the impact of this chunk rejection is, nor what the border effect might be. What is the correct way to apply an FFT to this chunked signal with discontinuities between the chunks?
Currently, the method I use to extract the alpha/delta power ratio is:
fs = 500. window = 5 # window length in seconds fft_freq = np.fft.rfftfreq(int(fs * window), 1.0/fs) alpha_band = np.where(np.logical_and(fft_freq>=8, fft_freq<=13)) delta_band = np.where(np.logical_and(fft_freq>=1, fft_freq<=4)) fftval = np.abs(np.fft.rfft(data, axis=1) / fs) alpha = np.average(np.multiply(np.abs(fftval[:, alpha_band]).T, weights)) delta = np.average(np.multiply(np.abs(fftval[:, delta_band]).T, weights))
weights is simply a
[0, 1] weight applied on each channel. As you can see, this is not exactly the power, but as I'm doing a ratio, I did not feel that the squaring was required. Feel free to correct me if I'm wrong.
If you have any comment, additional ideas for improvement, I'm also interested.