FFT of long signal by segments/chunks with discontinuities

Question

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.

Answer 1

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.