0
$\begingroup$

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

     
$\endgroup$
0
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ thank you so much for replying. Can you please tell me why it STFT ismn't simple putting a window on data? I know that there could be edge artifacts in the signal but how badly it will distort the signal? Or is it something else? $\endgroup$
    – Greyfrog
    Jul 11 at 3:26
  • $\begingroup$ @Greyfrog Edge effects can be severe with large windows; the solution is padding (I recommend reflect). Besides that, if you only seek power or magnitude and don't care for strict time-frequency resolutions, then yes it's as simple as fft(window * x[slice]) - else it involves properly designing the window relative to length of FFT, and for some applications, demodulating. $\endgroup$ Jul 11 at 3:50
  • $\begingroup$ I just wanted to clarify something. How can I apply FFT on real time data coming from sensor? $\endgroup$
    – Greyfrog
    Jul 12 at 13:16
  • $\begingroup$ @Greyfrog If by "real-time" you mean online/streaming, there's overlap-add. $\endgroup$ Jul 12 at 16:13
  • $\begingroup$ Yes, I mean real time data coming directly from the sensor over bluetooth and applying FFT on it. $\endgroup$
    – Greyfrog
    Jul 12 at 16:33

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.