eeg svm

import numpy as np

#import data
allData = np.random.randint(0, 100, (64, 256, 913))
labels_all = np.random.randint(2, size=913)


from scipy.signal import butter, filtfilt

freqs = {'theta': [4, 8], 'alpha': [8, 15], 'beta': [15, 30], 'gamma': [30, 50]}
chan_occipital = [26, 63, 29] #list of electrodes to filter for each brain region
chan_leftTemporal = [16, 14]
chan_rightTemporal= [53, 51]
filter_order = 5
sampling_frequency= 256.0

#generate filter coefficients for each frequency band. Apply these within the loop below to each electrode of interest
b_gamma,a_gamma = butter(filter_order, np.array(freqs['gamma'])/sampling_frequency, btype='bandpass')
b_theta,a_theta = butter(filter_order, np.array(freqs['theta'])/sampling_frequency, btype='bandpass')
b_alpha,a_alpha = butter(filter_order, np.array(freqs['alpha'])/sampling_frequency, btype='bandpass')
b_beta,a_beta = butter(filter_order, np.array(freqs['beta'])/sampling_frequency, btype='bandpass')

#set initial size of features array
features_all = np.zeros(shape=(allData.shape[2], 32))

for trial in range(allData.shape[2]):
    #extract power of specific bands from each electrode within each brain region. add their sum to the feature array
    runningPower = np.array([])

    #Occipital cortex
    for channel in chan_occipital:
        #extract gamma
        gammaPower = (filtfilt(b_gamma,a_gamma,allData[channel,:,trial])) ** 2
        gamma_firstHalf = sum(gammaPower[:allData.shape[1]/2])
        gamma_secondHalf = sum(gammaPower[allData.shape[1]/2:])
        #extract theta
        thetaPower = (filtfilt(b_theta,a_theta,allData[channel,:,trial])) ** 2
        theta_firstHalf = sum(thetaPower[:allData.shape[1]/2])
        theta_secondHalf = sum(thetaPower[allData.shape[1]/2:])

        runningPower = np.hstack([runningPower, gamma_firstHalf, gamma_secondHalf, theta_firstHalf, theta_secondHalf])

    #Left temporal cortex
    for channel in chan_leftTemporal:
        #extract alpha
        alphaPower = (filtfilt(b_alpha,a_alpha,allData[channel,:,trial])) ** 2
        alpha_firstHalf = sum(alphaPower[:allData.shape[1]/2])
        alpha_secondHalf = sum(alphaPower[allData.shape[1]/2:])
        #extract theta
        thetaPower = (filtfilt(b_theta,a_theta,allData[channel,:,trial])) ** 2
        theta_firstHalf = sum(thetaPower[:allData.shape[1]/2])
        theta_secondHalf = sum(thetaPower[allData.shape[1]/2:])

        runningPower = np.hstack([runningPower, alpha_firstHalf, alpha_secondHalf, theta_firstHalf, theta_secondHalf])

    #Right temporal cortex
    for channel in chan_leftTemporal:
        #extract beta
        betaPower = (filtfilt(b_beta,a_beta,allData[channel,:,trial])) ** 2
        beta_firstHalf = sum(betaPower[:allData.shape[1]/2])
        beta_secondHalf = sum(betaPower[allData.shape[1]/2:])
        #extract theta
        thetaPower = (filtfilt(b_theta,a_theta,allData[channel,:,trial])) ** 2
        theta_firstHalf = sum(thetaPower[:allData.shape[1]/2])
        theta_secondHalf = sum(thetaPower[allData.shape[1]/2:])
        #extract gamma
        gammaPower = (filtfilt(b_gamma,a_gamma,allData[channel,:,trial])) ** 2
        gamma_firstHalf = sum(gammaPower[:allData.shape[1]/2])
        gamma_secondHalf = sum(gammaPower[allData.shape[1]/2:])

        runningPower = np.hstack([runningPower, beta_firstHalf, beta_secondHalf, theta_firstHalf, theta_secondHalf, gamma_firstHalf, gamma_secondHalf])

    features_all[trial] = runningPower


#normalize feature matrix
normalize = True

if normalize == False:
    features_normed = (features_all - features_all.min(axis=0)) / (features_all.max(axis=0) - features_all.min(axis=0))
    features_all = features_normed


#PCA
usePCA = False

if usePCA == True:
    from sklearn.decomposition import RandomizedPCA

    n_components = 15

    pca = RandomizedPCA(n_components=n_components, whiten=True).fit(features_all)

    features_all = pca.transform(features_all)


#train model
trainModel = True

if trainModel == True:
    from sklearn.cross_validation import train_test_split
    from sklearn.svm import SVC

    features_train, features_test, labels_train, labels_test = train_test_split(features_all, labels_all, test_size=0.2, random_state=5)

    estimator = SVC(kernel='rbf', C=1, gamma=1)
    estimator.fit(features_train, labels_train)
    pred = estimator.predict(features_test)
    score_train = estimator.score(features_train, labels_train)
    score_test = estimator.score(features_test, labels_test)
    print "Training score: " + str(score_train)
    print "Test score: " + str(score_test)