0
$\begingroup$

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

With data a channel x samples array (e.g. 64x250) and 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)[0], :]

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))[0]
delta_band = np.where(np.logical_and(fft_freq>=1, fft_freq<=4))[0]

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))

The array 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.

$\endgroup$
1
$\begingroup$

I can think of two things to consider: discontinuities and phase.

The first is that whenever you take a finite chunk from some conceptually infinite signal, you get a step-discontinuity (ie. click) in the envelope of the signal, which generally results in some broadband energy. This can be reduced by multiplying the chunk with a window function, at the cost of slight blurring of the desired spectrum). I don't see any windowing in your code, so you might want to at least think about this. Standard windows might not be ideal though, see below.

The other thing to consider is that if you take a larger chunk and remove a part of it form the middle, then there will again be a click. If the two chunks are moved next to each other, then the phases from the two chunks probably won't align very well and the measured magnitudes will be lower in a frequency dependent way (probably unwanted, since this is generally "unpredictable").

That said, the phase-problems at least can be avoided if rather than removing a chunk, you simply set it to all zero instead so the time-alignment is maintained (ie. temporarily mute the signal; I apologize if you are already doing this, my Python is rusty). The muted part obviously won't contribute, but otherwise the only remaining impact I can think of should be the extra envelope transients.

To reduce the impact of the envelope transients (assuming the broadband energy causes problems), you might want to do some fade-in/fade-out, for example by using the cosine-tapers of a Tukey-window wherever it starts or ends (ie. around muted parts + beginning and end of the whole thing), though I don't understand your specific application well enough to really give more specific advice than this.

Finally I want to note about the squaring: while $a^2/b^2 = (a/b)^2$ the average of squared magnitudes (or the root of average of squared magnitudes) is not the same as the average of absolute magnitudes; the relative weighting of values is different. If you want to compute the average power, then you should square before averaging.

$\endgroup$
7
  • $\begingroup$ Thank you for the valuable input. I would like to understand a bit more about what you described. First of all, for the squaring, as I was not averaging across epochs/trials, I did not pay close enough attention. As you said, I am averaging across frequencies within a band and thus need to first raise the FFT to the power of 2. Secondly, is it best practice to apply a window to a signal, especially to a short signal, before applying the Fourier transform? I guess the window applied goal is to attenuate the signal on the edges (and eventually on artifacts), correct? $\endgroup$ – Mathieu Dec 27 '20 at 9:21
  • $\begingroup$ And finally, instead of replacing a bad chunk with 0, could I instead consider each chunk individually, keep the good ones, apply a window and FFT; and then average the FFT (squared for power ;) ) across chunks? In my sample code above, the line data.reshape(64, 5, 500) creates 5 chunks of 500 points, and I could avoided having to reconstruct a 64x2500 elements array if I can work directly on the chunks individually. $\endgroup$ – Mathieu Dec 27 '20 at 9:23
  • $\begingroup$ If you average (over frequencies) the absolute magnitudes, then you are computing the mean of the amplitudes (Manhattan norm divided by the bandwidth). If you compute the average of squares, you get the Euclidean norm divided by the bandwidth. The latter is preserved by FFT, so theoretically (ignoring windowing considerations) equivalent to computing the RMS (up to the final square root) of a band-pass filtered chunk in time-domain. $\endgroup$ – myzz Dec 28 '20 at 3:43
  • $\begingroup$ As for windowing, in some sense there is always one (the rectangular window is nothing else). Any sharp discontinuities at the edges where the signal is chopped are broadband events, which causes broadband "spectral leakage" in the FFT. This can be reduced by using a smoother window, but the tradeoff is some blurring of the spectrum (ie. some spectral resolution is lost). In most cases some window is desired, though which one depends on the trade-off desired. $\endgroup$ – myzz Dec 28 '20 at 4:02
  • $\begingroup$ As for processing each chunk separately and taking a moving average over the results of multiple chunks, sure. In this case though, average the two bands separately and take the ratio of the averages (because again, averaging before and after division will give you different results). This is a valid approach (and might actually be the best choice) if you only need a moving average over longer timespans and don't care for any resolution advantages of a larger FFT and will also make it easier to apply standard windows if it seems necessary. Hope this helps. $\endgroup$ – myzz Dec 28 '20 at 4:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.