Epilepsy 101 for Computational People

Disclaimer: This document cannot be used as any ground for medical and legal purposes. I write this up just because I am tired of explaining basic epilepsy terms to my colleagues (informatists, not neurologists) again and again.

1. Epilepsy is a neurological disorder while seizure is a characteristic symptom of epilepsy. Jerks are one kind of seizures.
2. In 1981, ILAE made a seizure classification that is widely used until today
   
   1. Focal Seizures: seizure limited in or originated from one hemisphere of the brain
      
      1. simple partial seizures, (simple: no interruption on consciousness)
      2. complex partial seizures, (complex: with interruption on consciousness) and
      3. partial seizure evolving into secondary generalized seizures, about 1/3 of partial seizures
   2. (Primary) Generalized Seizures: seizure all over the brain
      
      1. absence seizures (kinda like blackout), formerly called petit mal
      2. tonic-clonic seizures (kinda like whole-body jerk), formerly called grand mal
      3. many other generalized seizures
   3. Unclassified epileptic seizures
3. In 1989, ILAE made an epilepsy classification that is widely used until today
   
   1. Focal Epilepsies
      
      1. sympotomatic or cryptogenic (cause unknown): e.g., temporal lobe epilepsy (TLE)
      2. idopathic (genetic causes): e.g., benign childhood epilepsy
   2. Generalized Epilepsies
      
      1. sympotomatic or cryptogenic (cause unknown): e.g., West syndrome, Lennox-Gastaut syndrom
      2. idiopathic: e.g., childhood absence epilepsy, juvenile absence epilepsy,
   3. Epilepsies undertermined whether focal or generalized
4. When seizure is onset, we say the subject is in ictal state. Otherwise, interictal state. Because seizure duration (few seconds/minutes) is much shorter compared with interictal period, interictal EEG is more accessible.
5. EEG is a gold standard for epilepsy diagnosis. However, neurologists heavily rely on either ictal EEG or interictal epileptiform discharges (IEDs, such as sharp wave and spikes) in diagnosis. The IEDs are the distinctive EEG patterns for epilepsy when subjects are not in seizure. Because seizure and IEDs happen unpredictably and sporadically, the way to catch them is to use long-term EEG recordings which can last from hours to days. Such tedious procedure is costly and inconvenient.
6. Typical signal processing and machine learning problems related to epileptic EEG
   
   1. coarse epilepsy diagnosis: distinguishing epileptic interictal EEG (with or without IEDs) and non-epileptic EEG
   2. fine epilepsy/seizure diagnosis: identifying the type of seizure and/or epilepsy based on ictal or interictal EEG
   3. seizure detection: detecting seizure activities from epileptic subjects' EEG
   4. seizure prediction: predicting the onsets of seizure activities from epileptic subjects' EEG
   5. focus localization: locating the epileptogenic zone (for focal seizures only)

References:
1. For complete list of seizure types and epilepsy types, check Tables 9-1 and 9-2 on Pages 121-122 of this doc from British gov: http://www.nice.org.uk/guidance/cg137/resources/cg137-epilepsy-full-guideline3 which are reprinted from 1981 and 1989 ILAE classifications.

Signal Filtering using inverse FFT in Python

A straight forward way of doing signal filtering is zeroing out terms in inverse FFT result. I saw a good post online. However, it does not encapsulate into a function nor allow users to specify passing bands in terms of physical frequency. So I modified it into a function below.

def bandpass_ifft(X, Low_cutoff, High_cutoff, F_sample, M=None):
    """Bandpass filtering on a real signal using inverse FFT
    
    Inputs
    =======
    
    X: 1-D numpy array of floats, the real time domain signal (time series) to be filtered
    Low_cutoff: float, frequency components below this frequency will not pass the filter (physical frequency in unit of Hz)
    High_cutoff: float, frequency components above this frequency will not pass the filter (physical frequency in unit of Hz)
    F_sample: float, the sampling frequency of the signal (physical frequency in unit of Hz)    
    
    Notes
    =====
    1. The input signal must be real, not imaginary nor complex
    2. The Filtered_signal will have only half of original amplitude. Use abs() to restore. 
    3. In Numpy/Scipy, the frequencies goes from 0 to F_sample/2 and then from negative F_sample to 0. 
    
    """        
    
    import scipy, numpy
    if M == None: # if the number of points for FFT is not specified
        M = X.size # let M be the length of the time series
    Spectrum = scipy.fft(X, n=M) 
    [Low_cutoff, High_cutoff, F_sample] = map(float, [Low_cutoff, High_cutoff, F_sample])
    
    #Convert cutoff frequencies into points on spectrum
    [Low_point, High_point] = map(lambda F: F/F_sample * M /2, [Low_cutoff, High_cutoff])# the division by 2 is because the spectrum is symmetric 

    Filtered_spectrum = [Spectrum[i] if i >= Low_point and i <= High_point else 0.0 for i in xrange(M)] # Filtering
    Filtered_signal = scipy.ifft(Filtered_spectrum, n=M)  # Construct filtered signal 
    return Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point

Here is an example showing its usage.

First of all, let's generate a signal of different frequencies.

import numpy
N = 400 # signal length of number
x = numpy.arange(0, N, 1) # generate the time ticks
Sines = [numpy.sin(x*n)*(1-n) for n in [.9, .75, .5, .25, .12, .03, 0.025]]# different frequency components
y = numpy.sum(Sines, axis=0)   # add them by column, low frequencies have higher amplitudes

import matplotlib.pyplot as plt
plt.plot(x, y) # visualize the data

The raw signal looks like this:



Then, let's assume the sampling frequency is 500Hz and we only want signals between 5 to 30 Hz to go thru the filter.

Low_cutoff, High_cutoff, F_sample = 5, 30, 500
Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point = bandpass_ifft(y, Low_cutoff, High_cutoff, F_sample)

# Below is visualization
fig1 = plt.figure()
plt.stem(numpy.fft.fftfreq(N)*F_sample, abs(Spectrum))
plt.axis([0, F_sample/2, None, None]) # since the signal is real, the spectrum is symmetric|
plt.axvspan(Low_cutoff, High_cutoff,fc='g',alpha='.5') # The frequencies that we wanna get rid of

fig2 = plt.figure()
plt.plot(x, Filtered_signal)

In the first figure plotted, the green box shows frequencies that will be kept by the filter.



Here is the filtered signal. As you can see, high frequencies are gone. The signal looks very smooth.


Of course, this approach is very costly: FFT and iFFT are both $O(n\lg n)$ complexity algorithms. It is much more costly than convolving signal with filter terms (what people normally do in digital signal processing). Oh, and it cannot be done in real-time - you must wait until all data is sampled.

There is also a minor implementation issue here. Both Scipy and Numpy provides FFT/iFFT functions. But their speeds vary [1]. If you wanna go faster, try Python binding to some C/C++-based FFT libraries.

References:
1. http://stackoverflow.com/questions/6363154/what-is-the-difference-between-numpy-fft-and-scipy-fftpack