Extracting EEG Bands from raw EEG signal to match with the Bands calculated on the sensor level

Question

I have recently got my hands on a commercially available single channel EEG Sensor from Neurosky. This sensor acquire raw EEG signal at 512Hz from the frontal part of the human head and also process them on board level to various data such as -

1. EEG Bands: high alpha, low alpha, high beta, low beta....and so on.
2. Algorithms: Attention, Meditation
3. Raw Signal

My question is really simple but difficult to answer I suppose.

Problem:

I have captured all the signal and saved them into a csv file and loaded them in jupyter notebook. I am applying STFT on the raw signal. As I believe, STFT is nothing but FFT on window of the signal which in my case is 1 sec long window (512 data points). The window is overlapping with N data points.

Now I got my FFT plotting dynamically as if I am applying FFT on a real time signal. In the next step, I am trying to calculate band average of a given freq ranges -

freq = {
    'low_alpha':[7.5,9.25],
    'high_alpha':[10,11.75],
    'low_beta':[13,16.75],
    'high_beta':[18,29.75],
    'low_gamma':[31,39.75],
    'mid_gamma':[41,49.75],
    'delta':[0.5,2.75],
    'theta':[3.5,6.75]
} 

My problem is my band averages are not matching with the band avaerage calculate on the chip. Not even close. I am not expecting exact match rather something close which shows my understanding of the processing the signal is correct.

Here is how I am calculating the FFT -

def getFFT(signal,samplingFrequency):
    """Given some data and rate, returns FFTfreq and FFT (half)."""  
    signal_fft = np.fft.rfft(signal)
    fft_vals = np.absolute(signal_fft)
    psd = 20*np.log10(abs(signal_fft))
    return fft_vals,psd

And here is how I am plotting all my graphs with my calculations -

binSize = 512
x = np.arange(0,binSize,1)  
maxPlots = 3
fig,axs = plt.subplots(maxPlots,2) 

for i in range(maxPlots):  
    for j in range(2):
        plt.sca(axs[i][j]) 


dt = 5
t = 0
m = 0
calculated_band_Arr = [0] * binSize
actual_band_Arr = [0] * binSize

fft_freq = np.fft.rfftfreq(binSize, 1.0/binSize)

def plotFFT(i): 
    global t
    global m  
     
    # Sample a part of the signal 512
    sample = raw_signal[t:binSize+t] 
    fft,psd = getFFT(sample,binSize)  
     
    # plot sample of the signal
    axs[0][0].clear() 
    axs[0][0].plot(sample, color='c',linewidth=1,label="SIGNAL")  
    axs[0][0].legend(loc='upper left') 
    
    # plot fft of the signal
    axs[1][0].clear() 
    axs[1][0].plot(abs(fft[L]), color='c',linewidth=1,label="FFT")  
    axs[1][0].legend(loc='upper left') 
    
    # plot fft of the signal
    axs[1][1].clear() 
    axs[1][1].plot(abs(psd[L]), color='k',linewidth=1,label="PSD")  
    axs[1][1].legend(loc='upper left') 
    
    # Calculate band average according to my understanding and put it into an array
    # freq_ix = np.where((fft_freq >= freq[band_name][0]) & (fft_freq <= freq[band_name][1]))[0]
    # signal = butter_bandpass_filter(fft, freq[band_name][0], freq[band_name][1], binSize/2, 5)
    freq_res = fft[1] - fft[0]
    idx_delta = np.logical_and(fft >= freq[band_name][0], fft <= freq[band_name][1])
    band_power = simps(psd[idx_delta], dx=freq_res)
    calculated_band_Arr.pop(0)
    calculated_band_Arr.append(band_power)

    # Plot Selected Band calculated Manually
    # arr = np.array(calculated_band_Arr) 
    axs[2][0].clear()  
    axs[2][0].plot(calculated_band_Arr,color='r',linewidth=1,label="Calulated "+band_name)     
    axs[2][0].legend(loc='upper left')
    
    # Create Array with Selected Band from the sensor 
    actual_band_Arr.pop(0)
    actual_band_Arr.append(band_signal[m])
    
    # Plot Selected Band from the sensor
    axs[2][1].clear()
    axs[2][1].plot(actual_band_Arr,color='b',linewidth=1,label="Actual "+band_name)
    axs[2][1].legend(loc='upper left') 
    
    if m>= bandLen:
        m=0
    else:
        m=m+1
    
    if binSize+t>= bandLen:
        t = 0
    else:
        t = t+1
        
ani = FuncAnimation(plt.gcf(),plotFFT,interval=dt) 
plt.show() 

     

Answer 1

Proper STFT isn't simply putting a window on data and taking its FFT; I wouldn't recommend reinventing it unless knowing exactly what you're doing. There's open source implementations: librosa, ssqueezepy.

For matching against pre-computed values, it's important to account for any pre- or post-processing steps, such as baseline normalization or the log transform you took (relevant lecture).

I'd also recommend CWT over STFT for analyzing EEG - see this post.