Brain-Computer Interface


Charles Moyes (cwm55) and Mengxiang Jiang (mj294)


We built a robust Brain-Computer Interface (BCI) using single-channel electroencephalography (EEG) with an AVR microcontroller, and we were able to play Pong using our brain waves (and monitor/record our sleep).


Watch a demo video »

Introduction

Our goal was to build a brain-computer interface using an AVR microcontroller. We decided that the least invasive way of measuring brain waves would be using electroencephalography (EEG) to record microvolt-range potential differences across locations on the user's scalp. In order to accomplish this, we constructed a two-stage amplification and filtering circuit. Moreover, we used the built-in ADC functionality of the microcontroller to digitize the signal. Passive silver-plated electrodes soaked in a saline solution are placed on the user's head and connected to the amplifier board. The opto-isolated UART sends the ADC digital values over USB to a PC connected to the microcontroller. The PC runs software written in MATLAB and C to perform FFT and run machine learning algorithms (SVM) on the resultant signal. From there, we were able to control our own OpenGL implementation of the classic PC game Pong using our mind's brain waves. We also wrote software to record our sleep and store the EEG signal inside a data file.



Figure: Recording a User's Brain Waves Using EEG
Source: http://www.enotes.com/electroencephalogram-reference/eeg-machine


High-Level Design

Rational and Inspiration for Project Idea

Our project idea was inspired by Charles's severe obstructive sleep apnea (OSA) disorder. In order to diagnose sleep apnea, a clinical sleep study is performed where the patient is attached to EEG electrodes, along with SpO2, EMG, and respiration sensors. The patient's sleeping patterns are recorded overnight, and apneas (periods of sleep without breathing) can be identified within the collected data. This process is costly and requires an overnight stay at a hospital or sleep lab. Moreover, the patient often is denied access to their own data since a licensed sleep specialist interprets it for them. Our goal was to build a low-cost alternative that would allow users to take their health in their own hands by diagnosing and attempting to treat their own sleep disorders. Moreover, our project has diverse applications in the areas of neurofeedback (aiding meditation and treatment of ADHD disorder), along with brain-computer interfaces (allowing the disabled to control wheelchairs and spell words on a computer screen using their thoughts).


Background Math


Support Vector Machines

The machine learning algorithm we used was a support vector machine (SVM), which is a classifier that operates in a higher dimensional space and attempts to label the given vectors using a dividing hyperplane. The supervised learning method takes a set of training data and constructs a model that is able to label unknown test data. A brief explanation of the mathematics behind SVMs follows. During training, the SVM is given a set of instance-label pairs of the form $\{(\vec{x_i}, y_i) : i=1, \ldots, l\}$ where the instances are $n$-dimensional vectors such that $\vec{x_i} \in \mathbb{R}^n$. The $n$ dimensions represent $n$ separate "features." In addition, the labels are in the form $y \in \{1, -1\}$, where 1 and -1 designate target and non-target instances respectively. To "train" the support vector machine to recognize unknown input vectors, the following minimization problem is solved:

$$\min_{\mathbf{w},b,\mathbf{\xi}} \frac{1}{2}\mathbf{w}^T\mathbf{w}+C\sum_{i=1}^{l}\xi_i$$subject to:$$y_i(\mathbf{w}^T\phi(\mathbf{x}_i)+b)\geq 1 - \xi_i$$ $$\xi_i \geq 0$$Source: http://www.csie.ntu.edu.tw/~cjlin/papers/guide/guide.pdf


Note that $\phi$ is a function that maps the training vectors $\vec{x}_i$ into a higher-dimensional space, while $C>0$ and $\xi_i$ act as error terms (so-called "slack variables"). Moreover, $K$ is the kernel function which is defined as $K(\phi(\mathbf{x_i})^T\phi(\mathbf{x_j}))$. For our purposes, we used a radial basis function (RBF) kernel which has a $K$ function of $K(\mathbf{x_i}, \mathbf{x_j})=\exp(-\gamma||\mathbf{x_i}-\mathbf{x_j}||^2)$ where $\gamma > 0$ represents a user-tunable parameter.


DFT

The discrete Fourier transform (DFT) transforms a sequence of $N$ complex numbers (an $N$-point signal) in the time-domain into another $N$-sequence in frequency domain via the following formula:

$$X_k=\sum_{n=0}^{N-1}x_n \cdot e^{-2\pi\frac{k}{N}n}$$

The Fourier transform is denoted by $\mathcal{F}$, where $\mathbf{X}=\mathcal{F}(\mathbf{x})$. Source: http://en.wikipedia.org/wiki/Discrete_Fourier_transform

An algorithm, the Fast Fourier Transform (FFT) by Cooley and Tukey, exists to perform DFT in $O(n\log n)$ computational complexity as opposed to $O(n^2)$. We take advantage of this speed-up to perform DFTs in real-time on the input signals.


Filters

In order to filter the brain wave data in MATLAB, we use a finite impulse response (FIR) filter which operates on the last $N+1$ samples received from the ADC. In signal processing, the output $y$ of a linear time-invariant (LTI) system is obtained through convolution of the input signal $x$ with its impulse response $h$. This $h$ function "characterizes" the LTI system. The filter equation in terms of the output sequence $y[n]$ and the input sequence $x[n]$ is:

$$y[n]=h_0x[n]+h_1x[n-1]+\ldots+h_Nx[n-N]$$

Note that only $N$ coefficients are used for this filter (hence, "finite" impulse response). If we let $N \rightarrow \infty$, then the filter becomes an infinite impulse response (IIR) filter. Source: http://en.wikipedia.org/wiki/FIR_filter


EEG Signal Analysis

The EEG signal itself has several components separated by frequency. Delta waves are characteristic of deep sleep and are high amplitude waves in the frequency range $0 \leq f \leq 4$ Hz. Theta waves occur within the 4-8 Hz frequency band during meditation, idling, or drowsiness. Alpha waves have frequency range 8-14 Hz and take place while relaxing or reflecting. Another way to boost alpha waves is to close the eyes. Beta waves reside in the 13-30 Hz frequency band and are characteristic of the user being alert or active. They become present while the user is concentrating. Gamma waves in the 30-100 Hz range occur during sensory processing of sound and sight. Lastly, mu waves occur in the 8-13 Hz frequency range while motor neurons are at rest. Mu suppression takes place when the user imagines moving or actually moves parts of their body. An example diagram of the EEG signal types follows:



Figure: EEG Wave Frequency Ranges
Source: http://neurodevelopmentcenter.com/

EEG signals also contain event-related potentials (ERPs). An example is the P300 signal, which occurs when the user recognizes an item in a sequence of randomly presented events occurring with a Bernoulli distribution. It is emitted with a latency of around 300-600 ms and shows up as a deflection in the EEG signal:



Figure: P300 Event-Related Potential (ERP)
Source: http://en.wikipedia.org/wiki/Event-related_potential

Other artifacts present themselves in the EEG signal as well such as eye blinking and eye movement. An illustration of an example signal corrupted by eye blinking follows:



Figure: Eye Blinking within an EEG Signal
Source: http://psy.hull.ac.uk/Staff/m.large/Infonew/


Logical Structure

The overall structure of the project consists of an amplifier pipeline consisting of a differential instrumentation amplifier (where common-mode noise is measured using a right-leg driver attached to the patient's mastoid or ear lobe), along with an operational amplifier and some filters (to remove DC offsets, 60 Hz power-line noise, and other artifacts). From there, the signal passes to the microcontroller, where it is digitized via an ADC. Next, it is send over an isolated USB UART connection to a PC via an FTDI chip. The PC then performs signal processing and is able to output the results to the user, creating a neurofeedback loop which allows the user to control the PC using their brain waves. A functional block diagram of the overall structure follows:




Figure: High-Level Block Diagram of Project

Hardware/Software Trade-offs

Performing FFT in hardware using a floating-point unit (FPU) or a field programmable gate array (FPGA) would have allowed us to realize a considerable speed-up; however, our budget was only limited to $75, so this was not an option. Another trade-off we encountered was the use of MATLAB versus C. MATLAB is an interpretted language whose strength lies in performing vectorized matrix and linear algebra operations. It is very fast when performing these operations, but it is an interpreted language that does not run as native code. This speed penalty affected us when we attempted to collect data in real-time from the serial port at 57600 baud. To combat this speed penalty, I wrote a much faster OpenGL serial plotting application in C that runs at 200-400 frames per second on my machine (well above the ADC sample rate of 200 Hz) and is able to perform FFTs in real-time as the data comes in. Furthermore, yet another trade-off was the decision to use the PC to output the EEG waveforms rather than a built-in graphical LCD in hardware. Once again, budget constraints limited us, along with power usage since for safety reasons, our device uses four AA batteries instead of a mains AC power supply.


Relationship of Your Design to Available Standards

There exists a Modular EEG serial packet data format that is typically used to transmit EEG data over serial; however, we used ASCII terminal output (16-bit integer values in ASCII separated by line breaks) for simplicity, ease of debugging, and compatibility with MATLAB. Moreover, serial communications followed the RS232/USB standards. Another consideration was the IEC601 standard. IEC601 is a medical safety standard for medical devices that ensures that they are safe for patient use. Unfortunately, testing for IEC601 compliance was very much out-of-budget. Nevertheless, we discuss the many safety considerations that we absolutely adhered by in the Safety subsection (under Results) of this report.


Hardware Design

Amplifier Board Design

We built an analog amplification circuit with a total gain of 1,500 based on the design by Chip Epstein at https://sites.google.com/site/chipstein/home-page/eeg-with-an-arduino with modified gain and filter stages. The first stage uses an AD620 instrumentation amplifier for differential common mode signal rejection to reduce noise. The gain of the AD620 is approximately 23. A voltage divider and a 3140 opamp buffer provide a 2.5 V virtual ground for the instrumentation amplifier. After passing through the instrumentation amplifier, the signal is filtered using an RC high pass filter with $f_c = 0.13$ Hz (we modified the original design to allow the P300 ERP to reside within the pass-band of the filter).

Next, the signal undergoes a second-stage amplification. The gain of a 3140 opamp is set to approximately 65. The output signal is then filtered using an RC low-pass filter with a cut-off frequency of approximately 48 Hz. This frequency was chosen to preserve the low-frequency content of the EEG signal, while removing 50-60 Hz power line noise from the signal. We ordered parts from Digi-Key and samples from Analog Devices, Texas Instruments, and Maxim Semiconductor. We also sampled silver-plated passive EEG electrodes. From there, we were able to amplify a 125 $\mu\mathrm{V}\;\mathrm{V}_\mathrm{pp}$, 10 Hz square wave calibration signal to well within the ADC reference voltage range (0-1.1 $V$) and plot it on a PC in real-time. We constructed this prototype circuit on a breadboard. A schematic diagram of the amplifier board follows:



Figure: EEG Amplifier Board Schematic

Microcontroller Board Design

The microcontroller board contains a voltage divider that outputs a 125 $\mu\mathrm{V}\;\mathrm{V}_\mathrm{pp}$, 10 Hz square wave calibration signal from Pin D.3 of the microcontroller. Moreover, it contains a $+6 V$ DC battery-power supply that provides DC power to the microcontroller and the amplifiers. A schematic diagram of the microcontroller board follows:




Figure: Microcontroller Power Supply Schematic

Opto-Isolated UART over USB Design

We constructed an isolated +6 VDC power supply using 4 AA batteries and connected it the microcontroller using the PCB target board. We cut the ground trace connecting the microcontroller ground to the USB ground using a dremel tool. An illustration of the cut that we performed follows:



Figure: Ground Isolation Dremel Cut on Target Board

We used a Fairchild Semiconductor 6N137 optoisolator to isolate the USB power from the microcontroller power. The line of isolation is between the microcontroller UART RX and TX pins (Pin D.0 and Pin D.1) and the FTDI chip's RX and TX pins. A schematic diagram of the isolation circuit follows:




Figure: USB UART Opto-Isolation Schematic

Electrode Cap Design

We constructed an EEG helmet consisting of an old baseball cap modified to contain EEG electrodes. We followed the International 10-20 System of Electrode Placement by including electrodes at the designated locations on the scalp: Occipital lobe ($O$), Central lobe ($F_z$, $P_z$, $C_3$, $C_4$, $C_z$), and Frontal lobe ($F_{p_1}$, $F_{p_2}$, $G$). A diagram of the 10-20 system of electrode placement follows:



Figure: EEG Electrode Placement Diagram
Source: http://www.immrama.org/eeg/electrode.html


Software Design


MATLAB Serial Code

The primary function of the MATLAB serial code is to acquire digital EEG signal data from the microcontroller over the serial port. We wrote some code to plot the signal onto the screen and to perform rudimentary signal processing tasks (FFT and filtering). The MATLAB code consists of three files: plot_samples.m, plot_samples_rt.m, and serial_test.m




Figure: Calibration Signal Plotted within MATLAB


The serial_test.m script opens the serial port and displays an almost real-time plot of the serial data. It parses the serial data in a while loop via fscanf and str2num. Additionally, it updates the plot window contents using MATLAB's drawnow command. The loop terminates if the user closes the plot window, causing the script to clean up the file handle and close the serial port.




Figure: Actual EEG Signal Plotted within MATLAB (observe the beta waves and eye blink artifact)


The plot_samples.m script opens the serial port, and reads exactly $N=200$ samples of EEG data (around 2 seconds). It then closes and cleans up the serial port. Next, a 60 Hz notch filter is applied to the signal to remove powerline noise via the iirnotch and filter, and the DC offset (mean value) is subtracted from the signal. Finally, the time-domain signal is displayed in a plot window, along with the single-sided amplitude spectrum computed via MATLAB's fft function.

The plot_samples_rt.m script performs exactly the same operations as the plot_samples.m script, except it performs them in a loop. The script operates by sampling $N=200$ samples repeatedly until the user closes the plot window. As an effect, the signal plot and the frequency spectrum are refreshed every 2 seconds giving psuedo real-time operation.

OpenGL Plotter

A real-time serial plotter was written in OpenGL and C++. OpenGL is a high-performance graphics library, while C++ is much faster than MATLAB since it runs as native code. The program displays the real-time wave form from the serial port, along with its FFT. We use the FFTW (Fastest FFT in the West) library for computing the FFT using fast algorithms. Moreover, extensions were later added to the plotting code allow Pong to be played using brain waves, along with a P300 ERP detector. A data logging feature was added to allow us to record our EEG data while asleep to a file. The SDL library is used to collect user input from the keyboard and output the OpenGL graphics buffer to an on-screen window.




Figure: Calibration Signal Plotted within OpenGL


Initialization and Event Loop

The main() function initializes the SDL library, clears the ADC buffers, initializes OpenGL, initializes SDL TTF (for True Type font rendering), opens the serial port, opens the log file, and initializes the FFTW library. From there, the program enters the main event loop, an infinite while loop which checks for and handles key presses, along with drawing the screen via OpenGL.

The Quit() function cleans up SDL, closes the serial port, de-initializes FFTW frees the font, and quits the program using the exit UNIX system call.

The resizeWindow() function changes the window size of the screen by changing the viewport dimensions and the projection matrix. For this project, we use an orthographic projection with ranges $x \in [0, \mathrm{X\_SIZE}]$ and $y \in [0, \mathrm{ADC\_RESOLUTION}]$. The handleKeyPress() function intercepts key presses to quit the game (via the [ESCAPE] key) and to toggle full screen mode (using the [F1] key).

The initGL() function initializes OpenGL by setting the shading model, the clear color and depth, and the depth buffer.


Drawing Code

The drawGLScene() function is called once per frame to update the screen. It clears the screen, sets the model view matrix to the identity matrix $I_4$, shifts the oscilloscope buffer running_buffer to the left by one, and fills the new spot in the buffer with the newest sample from the serial port. This sample is obtained by calling the readSerialValue() function in serial.cpp. This function also contains logic to perform the FFT and draw it onto the screen as well. The sample is sent to the P300 module and logged to the output file. Moreover, the power spectrum of the FFT is computed using rfftw_one() and by squaring the frequency amplitudes.

The FFT bars and the oscilloscope points are plotted using GL_LINES and GL_POINTS respectively. Moreover, lines join adjacent oscilloscope points. The frequency ranges corresponding to each brain wave classification (alpha, beta, etc.) are calculated, along with their relative powers. Moreover, BCI code is executed which will be discussed in the OpenGL Pong sub-section of this report. The pong and P300 modules are then drawn on the screen, along with status text using TTF fonts and the framebuffer is swapped. Lastly, the frame rate is calculated using SDL_GetTicks(), an SDL library function.




Figure: Actual EEG Signal Plotted within OpenGL


Serial Code

The serial code handles serial communications over USB and is located in serial.cpp. Important parameters such as the buffer size BAUD_RATE, the port name PORT_NAME, and the baud rate B57600 are stored as pre-processor directives. The openSerial() function opens the serial port and sets the baud rate to 57600 (56k). The readSerialValue() function reads and parses one 10-bit ASCII ADC value from the serial port by scanning for a new-line terminator and using sscanf (readByte() is unused). Lastly, the closeSerial() function closes the serial port device. The UNIX system calls open, read, and close are used to carry out serial I/O, along with the GNU C library's system() function. The device file name for USB serial in UNIX is /dev/ttyUSB0.

Note that if the NULL_SERIAL pre-processor directive is set, then mock serial data is used rather than actually collecting data from the serial device. This functionality is useful for testing purposes.


Configuration

The config.h file contains pre-processor directives that can be used to configure the OpenGL plotting application during compile-time. Important parameters include the screen width and height (SCREEN_WIDTH and SCREEN_HEIGHT respectively), the screen bit depth (SCREEN_BPP), ADC_RESOLUTION the ADC resolution (the number of y-values) set to $2^{10}=1024$, and the log file name LOG_FILENAME (defaulting to "eeg_log.csv").


Debugging

Useful utility functions for debugging purposes are found in debug.cpp. The FileExists() function returns a boolean indicating whether the given file exists. The OpenLog() and CloseLog() functions are useful for writing log files with time and date stamps. The log_out() function can be passed a format string which is written to the log file, along with a time stamp. The format() function takes a format string and returns the formatted result. It uses the vformat utility function to generate the format string based on the arguments passed to the function. The dump() function dumps regions of memory to the screen in a human-friendly hexadecimal format.


Font Rendering

TTF font rendering support was implemented using NeHe sample code (located in the References section of this report). The glBegin2D() and glEnd2D() functions set up an orthographic screen projection for font rendering. Meanwhile, power_of_two() is a utility function that calculates the next largest power-of-two of the given integer input (useful for computing OpenGL texture dimensions which must be powers of two). The SDL_GL_LoadTexture() function converts an SDL_Surface structure in memory to an OpenGL texture object. This is useful because the SDL_TTF library only returns an SDL_Surface, but we are rendering the fonts using OpenGL. The InitFont() and FreeFont() functions use the SDL_TTF library functions to load and free fonts respectively. Lastly, glPrint() acts like an OpenGL implementation of printf by printing strings to the screen using TTF font textures in SDL_TTF.


OpenGL Pong




Figure: Screenshot of OpenGL Pong game


The Pong game logic and drawing code is located in the pong.cpp source file. Moreover, a set of constant "glyphs" is located in the glyphs.h header file. These glyphs are 2D arrays of boolean values corresponding to pixel on and pixel off for each scoreboard number that is displayed on the screen. The blocky font used for each glyph lends a retro syling to the game (which suits the Pong aesthetics much more than the TTF font renderer would).


Pong Game Logic

The game logic defines several variables including the paddle width and height (stored in PADDLE_WIDTH and PADDLE_HEIGHT respectively), the positions of both paddles, the $x$ and $y$ positions and velocities of the ball (ballposx, ballposy, ballvelx, ballvely). The sprite[] array includes a boolean bitmap storing the ball sprite.

The updateball() function performances rudimentary numerical integration to get the new ball position by addition. Bounds checks is performed and the ball velocity is negated to reflect the ball direction if the ball collides with the screen edges. Collision detection is performed with both paddles for all four sides also using bounds checking. Moreover, the ball's velocity is reflected and its speed is increased slightly by the constant value stored in vel_threshold upon colliding. Lastly, if the ball collides with the left or right side of the screen, the appropriate player's score variable (score1 or score2) is incremented depending on which side the ball collides on. The ResetVelocity() function is invoked to reset the ball's speed back to the initial speed.

The pongInit() function determines a random direction of the pong ball for the initial "face off" by calling a utility function randomNum() which generates a random integer in the set ${0, 1}$. The pongHandleKeyDown and pongHandleKeyUp update the game state based on keyboard presses and depresses. The [UP] key moves the right paddle up, and the [DOWN] key moves the right paddle down. The left paddle is controlled by the user's brain waves. The updatepaddle1() and updatepaddle2() functions perform numerical integration using addition to update the paddle $y$ velocities (paddle1yvel and paddle2yvel respectively).

Pong Game Drawing Code

The pongUpdateAndRender() function is invoked from the OpenGL plotter drawing code. It invokes all of the updating functions (updateball(), updatepaddle1(), updatepaddle2()) and all of the drawing functions (drawsprite(), drawpaddlesprite(), drawscore(), and drawline()). The drawscore() function draws both score glyphs onto the screen at the proper positions by invoking drawglyph(). The drawglyph() function takes an integer number, an $x$ and $y$ position, and an $(r,g,b)$ floating point color value. It uses GL_POINTS to draw each pixel of the glyph onto the screen. The drawsprite() function takes an $(x,y)$ position and an $(r,g,b)$ color and draws the ball sprite at that location. The drawline() function draws a striped line with spacing defined by LINE_SPACING across the center of the screen depicting the center of the field. Finally, the drawpaddlesprite() function takes an $(x,y)$ position and an $(r,g,b)$ color and draws the rectangular paddle sprite using GL_POINTS.


Pong Brain Wave Control and Brain-Computer Interface (BCI) Code

The most important part of the Pong software is the code snippet in main.cpp that updates the left paddle velocity based on the user's brain waves. We provide two modes of control. The first, alpha rhythm (8-13 Hz) modulation, provides proportional control based on the user's alpha waves. The EEG electrodes are placed on the user's forehead (near frontal lobes) during alpha rhythm measurement. The user concentrates to move the paddle down and relaxes to move the paddle up.

The other method of control is based on mu rhythm (8-13 Hz) suppression. The user imagines moving their feet up and down (or actually moves them) to move the paddle, and if mu suppression reaches a threshold, then the paddle moves down; otherwise, it moves up. The user places the electrodes on the top of the scalp near the sensorimotor cortex (10-20 locations $C_3$ and $C_4$) during mu suppression measurement. Although both methods worked equally well, we found during user testing that users preferred the alpha modulation control method over the mu rhythm suppression method.


Alpha Rhythm Modulation

The alpha rhythm modulation control is determined by two boundary variables ALPHA_MIN and ALPHA_MAX. The paddle's $y$ position posy is proportional to the alpha rhythm's relative power spectrum's value within this range $[\mathrm{ALPHA\_MIN}, \mathrm{ALPHA\_MAX}]$. From our testing, we found that values of $0.01$ and $0.04$ worked best for ALPHA_MIN and ALPHA_MAX respectively. Users were able to control the paddle position quite accurately after some practice.


Mu Rhythm Suppression

The mu rhythm supression control is determined by one threshold value MU_THRESHOLD. The paddle's $y$ velocity paddle1yvel is set to a value of $0.1$ if the mu rhythm's relative power spectrum is below MU_THRESHOLD (indicating mu suppression, or movement visualization in the user). The paddle1yvel is set to $-0.1$ if the mu rhythm's relative power spectrum exceeds or matches MU_THRESHOLD. Users were also able to use this method of control albeit less successfully due to the weaker signal received from placing the saline electrodes on the user's hair rather than their forehead. Nevertheless, mu rhythm suppression was also a viable control scheme.


Neurofeedback and Cursor Control

The alpha modulatiom and mu supression control schemes have diverse applications beyond simply playing the game Pong in brain-computer interfaces. Wheelchair and cursor control (both 1D and 2D) have been accomplished by mu rhythm suppression. In one instance, users controlled a cursor in 2D by imagining clenching either their left hand, their right hand, or moving their feet. This control scheme requires three channels measure three locations of the sensorimotor cortex near the top of the scalp: user's left side ($C_3$), center ($C_z$), and user's right side ($C_4$). Even though we had one channel, we could easily extend this to support 2D cursor control, along with detecting eye blinking artifacts for "clicking" the mouse. One could imagine applying this technology to allow users with special needs to control computer mouse movement.

The other application is in the field of neurofeedback. Neurofeedback creates a feedback loop for users attempting to meditate or treat ADHD disorder. The user visually sees or audibly hears the power of their alpha waves and is able to manipulate their alpha intensity. This neurofeedback has applications in the military and aircraft control as well, as users can be trained to focus and are alerted if they lose concentration. The Pong game can be viewed as a neurofeedback device since the user's concentration level is visually depicted on the screen as the position of the left paddle. Thus, the Brain-Computer Interface component of this project has diverse applications that go far beyond playing a simple computer game with one's brain waves.




Figure: Photograph of User Interaction with OpenGL Pong game

P300 Detector

The P300 detection code was an attempt to detect which color a user is thinking of from a discrete set of randomly flashed colors displayed on the screen. The software used machine learning algorithms for support vector machines (SVMs) provided by the libSVM C library. This attempt was not successful. In a training set of 50 trials, we were unable to obtain classification accuracy beyond 64%. Nonetheless, we document our code here and provide some suggestions for improvements and future work.

The P300 code uses a finite state machine (FSM) to display colors randomly on the screen in either a training mode or a testing mode. The colors are chosen randomly from the set {red, green, blue, yellow, purple}, and after each color in the set has been displayed exactly once, a trial is considered to have been completed. Five trials are performed during both testing and training, and the recorded brain waves are averaged. The idea is to attempt to classify one-second sets of brain data as either containing a P300 potential (target) or not (non-target) using the SVM. The target set corresponds with the color that the user is thinking of. While the colors are flashed on the computer screen, the user is instructed to count the number of times the target appears.


Code Structure

The code contains pre-processor definitions (TRAINING_DATA_FILENAME and TESTING_DATA_FILENAME) for the data file names, along with configuration variables for the number of colors NUM_COLORS, the number of trials (NUM_TRIALS), and the buffer size (BUFFER_SIZE).

The color_choices array contains $(r,g,b)$ floating-point tuples for each of the NUM_COLORS colors. The color_names array contains the human-readable names of each color (for text display). The trialBuffer is a 3D array indexed by trial number, color index, and buffer position containing the EEG waves recorded for each color and each trial. Moreover, targetBuffer and nonTargetBuffer contain the averaged target and non-target EEG waves respectively in training mode (to provide an additional target and non-target training instance to the SVM). Meanwhile, testBuffers is a 2D array of EEG wave buffers for each color used during testing mode (to test each color individually using the SVM).

State variables include current_color, the current color being displayed, and trainingTarget, the randomly-selected target color to be "chosen" by the user during training mode. The bufferPtr variable contains an index into the current trialBuffer. It gets incremented as additional samples are received from the serial port. The current_trial variable contains the index corresponding to the current color trial $\in [0, \mathrm{NUM\_TRIALS})$. The color_histogram variable contains an array of booleans signifying whether color $i$ has been displayed in the current trial yet or not.

The p300init() function initializes the P300 module by clearing the histograms and the buffers. It initially sets the training target and calls a placeholder function to train the SVM. For our purposes, we used libSVM's provided scripts to process the data, rather than directly integrating it within our code. This worked because our code merely generates data that can be used by libSVM offline for training and testing the support vector machine. The data files are updated whenever a new training or testing instance is provided, and then they can be used later by the user with libSVM.

The p300AddSample() function adds a new sample collected from the serial port to the current trial buffer. If the buffer has been filled, it starts a new trial, and if the last trial has finished, then the state is reset to the initial state. The background color is reset to black, and p300_state is set to P300_READY).

In both training and testing mode, the p300StartTrial() function is used to start a new trial. If the clear parameter is set, then this trial is considered to be the first trial, and state variables are reset accordingly. Otherwise, we check if all colors have been displayed. If they have been, then we increment current_trial and clear the histogram. We then choose the next random color and set its value in the histogram to true. Lastly, we update the OpenGL screen clear color to set the background color to the new randomly-chosen color.


User Interaction Code

The p300UpdateAndRender() function updates the screen to contain status text corresponding to the state of the P300 FSM (ready, training, testing mode), along with the training target. The SDL_ttf library functions from font.cpp are used to render the text onto the screen using OpenGL.

The p300HandleKeyDown() function checks for key presses and handles them accordingly. If the [F2] key is pressed, then the P300 FSM is switched to "training mode." If the [F3] key is pressed, then the P300 FSM is switched to "testing mode." Note that if a trial is already running, then nothing happens.


SVM Training Code

The p300setTrainingTarget() function sets the next training target of the training session. Initially, we randomly chose a color. However, we found that user fatigue is greater if the color changes during training, so we instead fixed the color index to correspond with the yellow color for all training sessions. Because training instances are only distinguished by their label (i.e. target or non-target), this does not have an effect on the training procedure, other than the fact that it is easier to concentrate on a single color throughout the entire training session rather than a randomly changing color.

The p300addTrainingExample() function constructs a targetBuffer and the nonTargetBuffer from the trialBuffer by averaging the target and non-target buffers. The data is then scaled for improved SVM performance in the range of $[-1,+1]$. Last, the training instance is appended to the TRAINING_DATA_FILENAME file using the ASCII format specified in the libSVM README file.


SVM Testing Code

The p300testandReport() function clears the testBuffers and stores the average EEG wave for each color throughout all trials from the trialBuffer array in each testBuffer. Then, the data is written to a testing data file specified by TESTING_DATA_FILENAME using the libSVM data format.


AVR Firmware

The AVR firmware was developed using the Atmel AVR Studio 5 integrated development environment (IDE) software program. The latest stable release of the WinAVR compiler was used, along with the ATAVRISP2 programmer.


Initialization Code

The firmware's C source code initializes the microcontroller ADC to use a $1.1 V$ reference voltage by setting the REFS1 bit of the ADMUX register. Next, the ADEN and ADIE bits are set in the ADCSRA register, enabling the ADC and the ADC interrupt respectively. A prescaler value of 125,000 is used. The LED port Pin D.2 and Pin D.3 are set to outputs. Next, Timer 0 is set to a 1 ms time base by setting a prescaler of 64 and a compare match ISR with an OCR0A value of 249 (implying 250 time ticks).

ADC sleep is enabled by setting the SM0 bit in the SMCR sleep mode control register. The UART is then initialized by calling uart_init(), sleep is enabled via sleep_enable(), and interrupts are enabled using sei().

From there, the firmware enters an infinite while loop. The CPU is put to sleep via sleep_cpu(), which automatically starts an ADC conversion. When the CPU wakes up, the current value of the ADC is sent via UART using fprintf(). A delay loop waits for the UART transmission to finish by delaying 1 ms until the UDRE0 bit is set in the UCSR0A register.


ADC Interrupt Handler

The ADC interrupt handler ADC_vect reads a 10-bit ADC value by storing the contents of the ADCL register inside a temporary variable, then reading ADCH, and computing ADCL+ADCH*256. Note that the register reads must be performed in this order, otherwise the ADC will lock up (see the Mega644 datasheet for details).


Timer 0 Compare ISR Interrupt Handler

The Timer 0 compare ISR vector TIMER0_COMPA_vect generates a test square wave on Pin D.2 and Pin D.3. It executes at a rate of 10 Hz, toggling the values of both pins in the PORTD register. Note that this task runs once every 10 ticks because of the 1 ms time base. We later disabled this feature because it was introducing extraneous noise in the EEG signal.


Results



Figure: Mengxiang using the EEG Device

Speed of Execution

The MATLAB serial plotter was not fast enough for real-time; however, the OpenGL plotter runs at 200-400 frames per second (FPS). Code was written in the OpenGL software to measure and output the FPS rate to the terminal console. This benchmark includes running Pong and doing a real-time FFT on the incoming serial data concurrently. The microcontroller samples the ADC at a rate of 200 Hz, and it sends the serial data to the PC over USB/UART at 57,600 (56k) baud.


Accuracy

The ADC code achieves 10-bit accuracy, utilizing ADC sleep to power down the CPU clock in between successive analog voltage measurements in order to reduce clock noise. Accuracy is highly dependent on electrode placement and electrical connection to the user's head. We discuss these measurement accuracy issues further in the Interference sub-section. Unfortunately, we were unable to attain acceptable accuracy performing SVM classification of P300 target wave forms, but one-dimensional BCI cursor control in Pong works perfectly for both alpha rhythm modulation (concentrating and relaxing) and mu suppression (visualizing motor movement).


Safety

Because this is a device attached to the user's head, safety was our utmost priority. The microcontroller was only ever powered by four 1.5V AA batteries, rather than through an AC power supply connected to mains. Moreover, serial communication to a PC over USB was isolated from the USB using optocouplers, which we tested extensively using a multimeter to ensure that both ground loops were separated. Only laptops running off of battery power supplies (no AC adapters connected) were ever connected to the microcontroller over USB. As a corollary, the microcontroller was never connected to a user's head while the programmer cable was connected to a PC. We promised that 120 VAC power will never be connected to this project directly or indirectly. As a result, users were never allowed to touch any other electrical devices while wearing the EEG helmet. We took safety very seriously throughout the development of the project.


Interference



Figure: 60 Hz Noise Corrupting an EEG Signal

Our device is very susceptible to interference from outside sources. We constructed DIY shielded electrode cables using aluminum foil, but we still encountered problems with 60 Hz power-line noise occasionally. Moreover, when the electrodes do not make sufficient electrical contact with the user's scalp, galvanic voltages show up that corrupt the signal. These galvanic voltages are less of an issue while playing Pong because an additional 50 Hz or 60 Hz frequency band does not affect the relevant frequency content of the signal. However, for time-domain analysis of the P300 ERP, these galvanic voltages can severely degrade the accuracy of this process.




Figure: DIY Shielded Aluminum Foil Electrode Cables

We also noticed that our 10 Hz calibration test square wave signal was introducing noise into the EEG signal even when it was not connected to the instrumentation amplifier, so we disabled this feature in the firmware. This solution eliminated the noise.


Usability and Special Needs Considerations

This project will have great societal impact because it is specifically designed for users with special needs. Using brain-computer interfaces, users with special needs will be able to interact in computers in ways that were not previously possible merely by using their brain waves. In addition, patients with sleep apnea will be able to collect and analyze their own EEG data while asleep without having to participate in expensive overnight sleep studies at hospitals. They will be able to see their data; rather than being shielded from it by a medical doctor. Because our budget is less than $75, they will be able to do this at a very low cost.


Conclusions

Expectations Met

When we were initially brainstorming for the project, we wanted a full sleep apnea diagnostic machine. This included not only an EEG, but also heart rate monitor, blood oxgen level monitor, temperature, etc. However, we immediately realized that doing just the EEG was a very challenging task by itself and so decided to focus on only that. When we actually created the proposal for the project, we had many ambitious ideas for what the EEG would do. Some of the ideas were fully implemented such as reading the Alpha Waves for checking how relaxed or focused you are, and the Mu Waves which correspond with thoughts of motion. Some of the ideas were somewhat implemented such as the reading of the P300 signal for identifying color. And some of the ideas were not implemented at all such as mouse cursor movement, due to time constraints and technical limitations such as only having 10-bit resolution and one channel, while most commercial EEGs have much more channels and better resolution. If we were to do this project again, we might add more channels and try to write a better P300 training and prediction program.


Conformed Standards

Although there is not an IEEE Standard for EEG operation, there are medical guidelines for specifications of what the EEG must have. The guidelines can be found at ftp://ansuk.org/pub/clinical_governance/dig_eeg.pdf

The first specification is that the EEG must have a minimum of 25 channels, preferably 32.
We did not meet this specification because we could not afford that many.

The input impedance also must be greater than 10 mega-ohms.
Our AD620 input impedance is 10||2 giga-ohms, which meets the specification.

The Noise must be below 2 microVpp from 0.16 to 100 Hz.
We used function generator input and we had noise around 3 microVolt, so we did not meet the specification.

The Common Mode Rejection Ratio must be greater than 80-100 dB.
The AD620 Common-Mode Rejection Ratio has a minimum of 110 dB for the gain and the range of frequencies measured, so the specification is met.

The sampling rate must be a minimum of 200 Hz, preferably 250-400 Hz.
Ours was 200 Hz, so we did meet this specification.

The dynamic range must be better than 2 microVolts.
We calculated our dynamic range to be around 0.7 microVolts, therefore meeting specification.

There must also be low pass filters of 15, 30, 50, 70, 100 Hz, and high pass filters of 0.16, 0.5, 1.6, 5, 10 Hz.
We had one low pass filter with 36 Hz, and a high pass filter of 0.13 Hz, therefore not enough filters to meet specification.

A Notch Filter with attenuation ratio 1:20 at 50 and 60 Hz is also required.
We did not have a hardware Notch Filter although we implemented one in software on MATLAB.

The electrode impedances must be displayed for each electrode and saved for review and playback.
We do not measure and display the electrode impedances so we did not meet this specification.

The ADC resolution must be equal to or greater than 12 bits.
Our ADC only has 10 bit accuracy, therefore not meeting this standard.

There must be enough memory storage, at least 1 GB.
Our computers have more than that so we meet the specification.

Inputs must be safely isolated and complies with the IEC 601-1/EN 60601-1 Type BF UL 544 Isolated standard.
We do not believe we conform to these standards; however, the inputs are isolated.

There must be a color monitor with at least 17" with minimum of 1280 x 1024 resolution.
We have a 22" color monitor with 1680 x 1050 resolution.

There must be a printer that prints at paper speed of 1, 2, 5, 10, 15, 30, 60 mm/sec.
We do not have such a printer, so we do not meet this specification.

Although we did not meet the majority of these specifications, we did not intend this device for medical use as we are simply students creating a microcontroller project in a month. Therefore, it is difficult if not impossible for us to create a product with the same amount of quality required for medical use.


Intellectual Property Considerations

We reused some code from the web to help us write the programs needed for analyzing the EEG.

The following are all public domain and open source therefore free to use and modify as needed. The ADC Sleep code for the microcontroller comes from the ECE 4760 Lab 2 example code by Professor Land. This code was used to increase the resolution of the ADC to 10-bits by putting the microcontroller to sleep and thus reduce noise. UART code, which was used in all the previous labs, is written by Joerg Wunsch. This was used to help debug the hardware by allowing the microcontroller to send messages to the computer using the usb. The code has a beerware license, which requires only that the license header remain on the file when used. When intially developing software for the EEG, we used MATLAB, and one of the first problems we faced was significant 60Hz noise. We looked online and found a MATLAB notch filter program which reduced the noise. When we found out that MATLAB was running too slow for real time, we turned to C instead. For plotting in C, we used OpenGL, which has the open source license. We used the examples from the Neon Helium (NeHe) OpenGL website as guides for creating programs in OpenGL. For writing text on the OpenGL window we used a sprintf program on the flipcode website.

We also used GNU licensed code in our project. The license states that we cannot modify the code, and we most provide the full source code if we distribute our code. As such, the full source code of our project can be found on GitHub at https://github.com/TheChuckster/EEG_BCI

One of the GNU licensed code is FFTW, a Fast Fourier Transform algorithm for C, which allowed as to quickly figure out the frequencies present in our EEG while plotting in OpenGL. Another is the GNU C Library (glibc), which provides the standard libary of functions for using C. The GNU Compiler Collection is the compiler system we used for our C programs. SDL_TTF is another GNU licensed code which was used for displaying truetype font.

We did not reverse engineering any design nor did we sign any non-disclosure agreements for any parts.


Ethical Considerations

When we started the project, we understood that we were hooking up electrical devices up to a person's head. If we did not carefully make the device and made some mistakes, the device could potentially harm or even kill a person. Therefore, safety was one of our main concerns when starting the project. Our first objective was to make sure that even if we made some mistake, there was still no way to harm a person. This was done by optoisolating the electrode circuit from the microcontroller, and also providing the microcontroller's power source with AA batteries instead of an outlet. This way, even if we shorted out the circuit somehow, the most voltage the person will receive will be 6V from 4 AA batteries.

When we were making the project and testing out our EEG, we got various signals. Although sometimes we thought the signals we got were actual brain waves, we also understood that we could be wrong. By careful calculations and procedures, we were able to distinguish noise from what could be brain waves. Even then, we were not sure if the signal we got was indeed brain waves, so we also asked a researcher at the Cornell Sleep Lab to verify that indeed we were right.


Legal Considerations


NOTE: This project is made by students and not medical professionals, therefore it is not intended for medical use!

We also made a PONG-like game, and the original PONG game was made by Atari. However, we are simply using the game as a demo of our EEG and not for commercial gain, therefore there should be no legal implications for having PONG.




Figure: Hard at work inside Phillips ECE Lab

Appendix

Sleep Recording


Because this device was originally intended to diagnose sleep-related disorders, we used to EEG to record Mengxiang taking a nap and discovered several interesting sleep-related waveforms. I took several screenshots of the OpenGL plotting application. Sleep spindles and k-complexes occur during the onset of stage 2 sleep and aid in maintaining muscle memory and sleep-related relaxation (See http://en.wikipedia.org/wiki/Sleep_spindle). We also observed theta waves during Stage 1 and Stage 2 sleep (http://en.wikipedia.org/wiki/Stage_2_sleep), and REM sleep (http://en.wikipedia.org/wiki/REM_sleep). I woke Mengxiang up during the observed REM waves on the OpenGL plot, and he confirmed that he was dreaming.

We included a series of screenshots (the green bars show the Fourier power spectrum):




Figure: Sleep Spindle



Figure: K Complex ($\frac{2}{3}$ of the way to the right)



Figure: Another K Complex (near the left, notice the strong deflections in voltage)



Figure: Theta Waves (Lower frequency)



Figure: REM sleep (Higher frequency)

Source Code Listings

You may find a copy of our source code on Git Hub at https://github.com/TheChuckster/EEG_BCI.


AVR Microcontroller Firmware


lab5.c
#include <inttypes.h>
#include <avr/io.h>
#include <avr/interrupt.h>
#include <avr/sleep.h>
#include <stdio.h>
#include <stdlib.h>
#include <util/delay.h>
#include "uart.h"

// UART file descriptor
// putchar and getchar are in uart.c
FILE uart_str = FDEV_SETUP_STREAM(uart_putchar, uart_getchar, _FDEV_SETUP_RW);

// timeout values for each task
#define t1 200
#define t2 5

volatile unsigned char time1=0, time2=0; // timeout counter
unsigned char led; // light states

volatile int Ain, AinLow; // raw A to D number

ISR (ADC_vect)
{
	//program ONLY gets here when ADC done flag is set
   //when reading 10-bit values 
   //you MUST read the low byte first 
   AinLow = (int)ADCL;
   Ain = (int)ADCH*256; 
   Ain = Ain + AinLow;
}

// timer 0 compare ISR
ISR (TIMER0_COMPA_vect)
{
	if (time1 > 0) --time1;
	
	if (time1 == 0)
	{
		time1 = t1;

		// toggle the LED and the square wave test output
		led ^= 1;
		PORTD = (led << PORTD2) | (led << PORTD3);
	}
	
	if (time2 > 0) --time2;
}

int main()
{
	// init the A to D converter 
	//channel zero/ right adj /ext Aref
	//!!!DO NOT CONNECT Aref jumper!!!!
	ADMUX = (1<<REFS1);  // 1.1 V ref
	//enable ADC and set prescaler to 1/127*16MHz=125,000
	//and set int enable
	ADCSRA = (1<<ADEN) | (1<<ADIE) + 7 ; 

	// set up the LED port
  	DDRD = (1 << PORTD2) | (1 << PORTD3); // PORT D.2 is an ouput

  	// set up timer 0 for 1 mSec timebase
  	TIMSK0 = 1 << OCIE0A; // turn on timer 0 cmp match ISR
  	OCR0A = 249; // set the compare register to 250 time ticks
  	TCCR0B = 3; // set prescalar to divide by 64
  	TCCR0A = 1 << WGM01; // turn on clear-on-match
	  
   SMCR = (1<<SM0) ; // sleep -- choose ADC mode 

  	 	
  	// set up timer for PWM
  	//pwm_init();

  	led = 0x00; // init the LED status

  	time1 = t1; // init the task timer
  	time2 = t2;
   
	// init the UART -- uart_init() is in uart.c
	uart_init();
	stdout = stdin = stderr = &uart_str;
	fprintf(stdout,"\n\rStarting ADC ISR demo...\n\r"); 

	// Need the next two statments so that the USART finishes
	// BEFORE the cpu goes to sleep.
	while (!(UCSR0A & (1<<UDRE0))) ; 
	_delay_ms(1);

	sleep_enable();                             
	sei();

	// measure and display loop
	while (1)
	{
	    sleep_cpu();
		fprintf(stdout, "%d\n\r", Ain);
		while (!(UCSR0A & (1<<UDRE0))) ;
		_delay_ms(1);

		/*if (time2 == 0)
		{
			time2 = t2;	
			fprintf(stdout, "%d\n\r", Ain);
		}
		
		if (ADCSRA & (1<<ADSC)) continue; // skip processing if not ready to do ADC -- now we don't need an interrupt

		AinLow = (int)ADCL;
		Ain = (int)ADCH*256; 
		Ain = Ain + AinLow;
			
		// start another conversion
		ADCSRA |= (1<<ADSC); */
	}
}
                 
uart.c
/*
 * ----------------------------------------------------------------------------
 * "THE BEER-WARE LICENSE" (Revision 42):
 * <joerg@FreeBSD.ORG> wrote this file.  As long as you retain this notice you
 * can do whatever you want with this stuff. If we meet some day, and you think
 * this stuff is worth it, you can buy me a beer in return.        Joerg Wunsch
 * ----------------------------------------------------------------------------
 *
 * Stdio demo, UART implementation
 *
 * $Id: uart.c,v 1.1 2005/12/28 21:38:59 joerg_wunsch Exp $
 *
 * Mod for mega644 BRL Jan2009
 */


/* CPU frequency */
#define F_CPU 16000000UL

/* UART baud rate */
#define UART_BAUD  57600


#include <stdint.h>
#include <stdio.h>

#include <avr/io.h>

#include "uart.h"

/*
 * Initialize the UART to 9600 Bd, tx/rx, 8N1.
 */
void
uart_init(void)
{
#if F_CPU < 2000000UL && defined(U2X)
  UCSR0A = _BV(U2X);             /* improve baud rate error by using 2x clk */
  UBRR0L = (F_CPU / (8UL * UART_BAUD)) - 1;
#else
  UBRR0L = (F_CPU / (16UL * UART_BAUD)) - 1;
#endif
  UCSR0B = _BV(TXEN0) | _BV(RXEN0); /* tx/rx enable */
}

/*
 * Send character c down the UART Tx, wait until tx holding register
 * is empty.
 */
int
uart_putchar(char c, FILE *stream)
{

  if (c == '\a')
    {
      fputs("*ring*\n", stderr);
      return 0;
    }

  if (c == '\n')
    uart_putchar('\r', stream);
  loop_until_bit_is_set(UCSR0A, UDRE0);
  UDR0 = c;

  return 0;
}

/*
 * Receive a character from the UART Rx.
 *
 * This features a simple line-editor that allows to delete and
 * re-edit the characters entered, until either CR or NL is entered.
 * Printable characters entered will be echoed using uart_putchar().
 *
 * Editing characters:
 *
 * . \b (BS) or \177 (DEL) delete the previous character
 * . ^u kills the entire input buffer
 * . ^w deletes the previous word
 * . ^r sends a CR, and then reprints the buffer
 * . \t will be replaced by a single space
 *
 * All other control characters will be ignored.
 *
 * The internal line buffer is RX_BUFSIZE (80) characters long, which
 * includes the terminating \n (but no terminating \0).  If the buffer
 * is full (i. e., at RX_BUFSIZE-1 characters in order to keep space for
 * the trailing \n), any further input attempts will send a \a to
 * uart_putchar() (BEL character), although line editing is still
 * allowed.
 *
 * Input errors while talking to the UART will cause an immediate
 * return of -1 (error indication).  Notably, this will be caused by a
 * framing error (e. g. serial line "break" condition), by an input
 * overrun, and by a parity error (if parity was enabled and automatic
 * parity recognition is supported by hardware).
 *
 * Successive calls to uart_getchar() will be satisfied from the
 * internal buffer until that buffer is emptied again.
 */
int
uart_getchar(FILE *stream)
{
  uint8_t c;
  char *cp, *cp2;
  static char b[RX_BUFSIZE];
  static char *rxp;

  if (rxp == 0)
    for (cp = b;;)
      {
	loop_until_bit_is_set(UCSR0A, RXC0);
	if (UCSR0A & _BV(FE0))
	  return _FDEV_EOF;
	if (UCSR0A & _BV(DOR0))
	  return _FDEV_ERR;
	c = UDR0;
	/* behaviour similar to Unix stty ICRNL */
	if (c == '\r')
	  c = '\n';
	if (c == '\n')
	  {
	    *cp = c;
	    uart_putchar(c, stream);
	    rxp = b;
	    break;
	  }
	else if (c == '\t')
	  c = ' ';

	if ((c >= (uint8_t)' ' && c <= (uint8_t)'\x7e') ||
	    c >= (uint8_t)'\xa0')
	  {
	    if (cp == b + RX_BUFSIZE - 1)
	      uart_putchar('\a', stream);
	    else
	      {
		*cp++ = c;
		uart_putchar(c, stream);
	      }
	    continue;
	  }

	switch (c)
	  {
	  case 'c' & 0x1f:
	    return -1;

	  case '\b':
	  case '\x7f':
	    if (cp > b)
	      {
		uart_putchar('\b', stream);
		uart_putchar(' ', stream);
		uart_putchar('\b', stream);
		cp--;
	      }
	    break;

	  case 'r' & 0x1f:
	    uart_putchar('\r', stream);
	    for (cp2 = b; cp2 < cp; cp2++)
	      uart_putchar(*cp2, stream);
	    break;

	  case 'u' & 0x1f:
	    while (cp > b)
	      {
		uart_putchar('\b', stream);
		uart_putchar(' ', stream);
		uart_putchar('\b', stream);
		cp--;
	      }
	    break;

	  case 'w' & 0x1f:
	    while (cp > b && cp[-1] != ' ')
	      {
		uart_putchar('\b', stream);
		uart_putchar(' ', stream);
		uart_putchar('\b', stream);
		cp--;
	      }
	    break;
	  }
      }

  c = *rxp++;
  if (c == '\n')
    rxp = 0;

  return c;
}
                 
uart.h
/*
 * ----------------------------------------------------------------------------
 * "THE BEER-WARE LICENSE" (Revision 42):
 * <joerg@FreeBSD.ORG> wrote this file.  As long as you retain this notice you
 * can do whatever you want with this stuff. If we meet some day, and you think
 * this stuff is worth it, you can buy me a beer in return.        Joerg Wunsch
 * ----------------------------------------------------------------------------
 *
 * Stdio demo, UART declarations
 *
 * $Id: uart.h,v 1.1 2005/12/28 21:38:59 joerg_wunsch Exp $
 */

/*
 * Perform UART startup initialization.
 */
void	uart_init(void);

/*
 * Send one character to the UART.
 */
int	uart_putchar(char c, FILE *stream);

/*
 * Size of internal line buffer used by uart_getchar().
 */
#define RX_BUFSIZE 80

/*
 * Receive one character from the UART.  The actual reception is
 * line-buffered, and one character is returned from the buffer at
 * each invokation.
 */
int	uart_getchar(FILE *stream);
                 

OpenGL Pong/P300


config.h
#ifndef CONFIG_H_
#define CONFIG_H_

#define SCREEN_WIDTH  1024
#define SCREEN_HEIGHT  768
#define SCREEN_BPP      16

#define ADC_RESOLUTION 1024
#define FFT_SCALE_FACTOR (1024*1024*10)
#define X_SIZE 1024

#define TRUE  1
#define FALSE 0

#define LOG_FILENAME "eeg_log.csv"

#endif
                 
debug.cpp
#include "debug.h"
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <time.h>
using namespace std;

ofstream log_file;

#define BPERL 16 // byte/line for dump

bool FileExists(std::string filename)
{
	std::fstream fin;
	fin.open(filename.c_str(), std::ios::in);

	if(fin.is_open())
	{
		fin.close();
		return true;
	}

	fin.close();
	return false;
}

void OpenLog()
{
	string filename = "";

	int i=0;
	bool taken=true;

	do {
		i++;
		filename = format("logs/log%d.txt", i);
		if (FileExists(filename))
			taken = true;
		else
			taken = false;
	} while(taken == true);

	printf("Using log file %s ...\n", filename.c_str());
	log_file.open(filename.c_str(), ofstream::out | ofstream::trunc);
}

void CloseLog()
{
	log_file.close();
	printf("Log file closed.\n");
}

void log_out(const char* fmt ...)
{
	va_list argList;
	va_start(argList, fmt);
	std::string result = vformat(fmt, argList);
	va_end(argList);

	char timebuf[52];
	long timeval;

	time(&timeval);
	strftime(timebuf,32,"%I:%M:%S",localtime(&timeval));

	cout << "[" << timebuf << "] " << result;

	if (log_file.is_open()) log_file << "[" << timebuf << "] " << result;

	return;
}

std::string format(const char* fmt ...)
{
	va_list argList;
	va_start(argList, fmt);
	std::string result = vformat(fmt, argList);
	va_end(argList);

	return result;
}

std::string vformat(const char *fmt, va_list argPtr)
{
	const int maxSize = 1000000;
	const int bufSize = 161;
	char stackBuffer[bufSize];

	int attemptedSize = bufSize - 1;

	int numChars = vsnprintf(stackBuffer, attemptedSize, fmt, argPtr);

	if (numChars >= 0)
		return std::string(stackBuffer);	// Got it on the first try.

	char* heapBuffer = NULL;	// Now use the heap.

	while ((numChars == -1) && (attemptedSize < maxSize))
	{
		attemptedSize *= 2;	// Try a bigger size
		heapBuffer = (char*)realloc(heapBuffer, attemptedSize + 1);
		numChars = vsnprintf(heapBuffer, attemptedSize, fmt, argPtr);
	}

	std::string result = std::string(heapBuffer);
	free(heapBuffer);

	return result;
}

void dump(unsigned char *data, unsigned count)
{
	unsigned byte1, byte2;

	while(count != 0)
	{
		for(byte1 = 0; byte1 < BPERL; byte1++)
		{
			if(count == 0)
				break;

			printf("%02X ", data[byte1]);
			count--;
		}
		printf("\t");
		for(byte2 = 0; byte2 < byte1; byte2++)
		{
/*			if(data[byte2] < ' ')
				printf(".");
			else
				printf("%c", data[byte2]);*/
			if(data[byte2] >= '!' && data[byte2] <= '}')
				printf("%c", data[byte2]);
			else
				printf(".");
		}
		printf("\n");
		data += BPERL;
	}
}
                 
debug.h
#ifndef DEBUG_H_
#define DEBUG_H_

#include <string>
#include <stdio.h>
#include <cstdarg>
#include <stdarg.h>

bool FileExists(std::string filename);
void OpenLog();
void CloseLog();

std::string format(const char* fmt ...);
std::string vformat(const char *fmt, va_list argPtr);
void dump(unsigned char *data, unsigned count);
void log_out(const char* fmt ...);

#endif
                 
font.cpp
#include <math.h>
#include "font.h"

#define FONT_SIZE 16
SDL_Color red = {255,0,0}, blue = {0,0,255};

void glBegin2D()
{
	glViewport(0, 0, SCREEN_WIDTH, SCREEN_HEIGHT);

	glMatrixMode(GL_PROJECTION);
	glPushMatrix();
	glLoadIdentity();

    glOrtho(0.0, (GLdouble)SCREEN_WIDTH, (GLdouble)SCREEN_HEIGHT, 0.0, 0.0, 1.0);

	glMatrixMode(GL_MODELVIEW);
	glPushMatrix();
	glLoadIdentity();

	glPushAttrib(GL_ENABLE_BIT);
	glDisable(GL_DEPTH_TEST);
	glDisable(GL_CULL_FACE);
	glEnable(GL_TEXTURE_2D);

	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
}

void glEnd2D()
{
	glMatrixMode(GL_PROJECTION);
	glPopMatrix();
	glMatrixMode(GL_MODELVIEW);
	glPopMatrix();
	glPopAttrib();
}

static int power_of_two(int input)
{
	int value = 1;

	while (value < input)
	{
		value <<= 1;
	}
	return value;
}

GLuint SDL_GL_LoadTexture(SDL_Surface *surface, GLfloat *texcoord)
{
	GLuint texture;
	int w, h;
	SDL_Surface *image;
	SDL_Rect area;
	Uint32 saved_flags;
	Uint8  saved_alpha;

	w = power_of_two(surface->w);
	h = power_of_two(surface->h);

	texcoord[0] = 0.0f;			/* Min X */
	texcoord[1] = 0.0f;			/* Min Y */
	texcoord[2] = (GLfloat)surface->w / w;	/* Max X */
	texcoord[3] = (GLfloat)surface->h / h;	/* Max Y */

	image = SDL_CreateRGBSurface(SDL_SWSURFACE, w, h, 32,
#if SDL_BYTEORDER == SDL_LIL_ENDIAN /* OpenGL RGBA masks */
			0x000000FF,
			0x0000FF00,
			0x00FF0000,
			0xFF000000
#else
			0xFF000000,
			0x00FF0000,
			0x0000FF00,
			0x000000FF
#endif
		       );

	if (image == NULL) return 0;

	saved_flags = surface->flags&(SDL_SRCALPHA|SDL_RLEACCELOK);
	saved_alpha = surface->format->alpha;
	if ((saved_flags & SDL_SRCALPHA) == SDL_SRCALPHA) SDL_SetAlpha(surface, 0, 0);

	/* Copy the surface into the GL texture image */
	area.x = 0;
	area.y = 0;
	area.w = surface->w;
	area.h = surface->h;
	SDL_BlitSurface(surface, &area, image, &area);

	/* Restore the alpha blending attributes */
	if ((saved_flags & SDL_SRCALPHA) == SDL_SRCALPHA) SDL_SetAlpha(surface, saved_flags, saved_alpha);

	/* Create an OpenGL texture for the image */
	glGenTextures(1, &texture);
	glBindTexture(GL_TEXTURE_2D, texture);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, w, h, 0, GL_RGBA, GL_UNSIGNED_BYTE, image->pixels);
	SDL_FreeSurface(image); /* No longer needed */

	return texture;
}

TTF_Font *font;

void InitFont()
{
	font = TTF_OpenFont("fonts/arial.ttf", FONT_SIZE);
	if (font == NULL)
	{
		printf("TTF_OpenFont: %s\n", SDL_GetError());
		return;
	}
}

void FreeFont()
{
	TTF_CloseFont(font);
}

void glPrint(std::string message, int x, int y, SDL_Color color)
{
	int w, h;
	GLfloat texcoord[4];
	GLfloat texMinX, texMinY;
	GLfloat texMaxX, texMaxY;
	SDL_Surface *text;
	GLuint fonttexture;

	text = TTF_RenderText_Blended(font, message.c_str(), color);
	fonttexture = SDL_GL_LoadTexture(text, texcoord);
	SDL_FreeSurface(text);

	w = text->w;
	h = text->h;

	texMinX = texcoord[0];
	texMinY = texcoord[1];
	texMaxX = texcoord[2];
	texMaxY = texcoord[3];

	glBegin2D();
	glDisable(GL_LIGHTING);
	glBindTexture(GL_TEXTURE_2D, fonttexture);
	glColor3f(1.0,1.0,1.0);
	glBegin(GL_TRIANGLE_STRIP);
		glTexCoord2f(texMinX, texMinY); glVertex2i(x,   y  );
		glTexCoord2f(texMaxX, texMinY); glVertex2i(x+w, y  );
		glTexCoord2f(texMinX, texMaxY); glVertex2i(x,   y+h);
		glTexCoord2f(texMaxX, texMaxY); glVertex2i(x+w, y+h);
	glEnd();
	glEnable(GL_LIGHTING);
	glEnd2D();

	glDeleteTextures(1, &fonttexture);
}
                 
font.h
#ifndef FONT_H_
#define FONT_H_

#include <string>
#include <GL/gl.h>
#include <GL/glu.h>
#include "SDL.h"
#include "SDL_ttf.h"
#include "config.h"

void glPrint(std::string message, int x, int y, SDL_Color color);
void InitFont();
void FreeFont();

void glBegin2D();
void glEnd2D();

extern SDL_Color red, blue;

#endif
                 
glyphs.h
const unsigned char num0[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1
};

const unsigned char num1[] = 
{
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1
};

const unsigned char num2[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1
};

const unsigned char num3[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1
};

const unsigned char num4[] = 
{
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1
};

const unsigned char num5[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1
};

const unsigned char num6[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,0,0,0,0,0,0,0,0,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1
};

const unsigned char num7[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1
};

const unsigned char num8[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1
};

const unsigned char num9[] = 
{
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,0,0,0,0,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1,
0,0,0,0,0,0,0,0,1,1,1,1
};
                 
main.cpp
#include <stdio.h>
#include <stdlib.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <time.h>
#include <rfftw.h>
#include "SDL.h"
#include "font.h"
#include "debug.h"
#include "serial.h"
#include "config.h"
#include "pong.h"
#include "p300.h"

// drawing stuff
SDL_Surface *surface;
unsigned int point_buffer[X_SIZE], running_buffer[X_SIZE], xline=0;

// logging stuff
FILE* logFile = NULL;

// FFT stuff
fftw_real in[X_SIZE], out[X_SIZE], power_spectrum[X_SIZE/2+1];
rfftw_plan p;

void Quit(int returnCode)
{
	SDL_ShowCursor(SDL_ENABLE);

    printf("Quiting...\n");
    closeSerial();
    if (logFile) fclose(logFile);
    rfftw_destroy_plan(p);

    FreeFont();
    SDL_Quit();
    exit(returnCode);
}

int resizeWindow(int width, int height)
{
    GLfloat ratio;

    if (height == 0) height = 1;

    ratio = (GLfloat)width / (GLfloat)height;

    glViewport(0, 0, (GLsizei)width, (GLsizei)height);

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

    //gluPerspective(45.0f, ratio, 0.1f, 100.0f);
    gluOrtho2D(0, X_SIZE, 0, ADC_RESOLUTION);

    glMatrixMode(GL_MODELVIEW);

    glLoadIdentity();

    return TRUE;
}

void handleKeyPress(SDL_keysym *keysym)
{
    switch (keysym->sym)
	{
	case SDLK_ESCAPE:
	    Quit(0);
	    break;
	case SDLK_F1:
	    SDL_WM_ToggleFullScreen(surface);
	    break;
	default:
	    break;
	}

    return;
}

int initGL(void)
{
    glShadeModel(GL_SMOOTH);
    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
    glClearDepth(1.0f);
    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LEQUAL);
    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);

    return TRUE;
}

int drawGLScene(void)
{
    static GLint T0 = 0, Frames = 0;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glLoadIdentity();
    glTranslatef(0.0f, 0.0f, 0.0f);

    unsigned int s = point_buffer[xline] = readSerialValue();
    for (int i=1; i<X_SIZE; i++)
    {
        running_buffer[i-1] = running_buffer[i];
    }
    running_buffer[X_SIZE-1] = point_buffer[xline];

    // log the output to a file
    fprintf(logFile, "%d,%d\n", time(NULL), point_buffer[xline]);

    // send sample to P300 buffer
    p300AddSample(s);

    xline++;
    xline %= X_SIZE;

        for (int i=0; i<X_SIZE; i++)
        {
            in[i] = running_buffer[i] - ADC_RESOLUTION/2;
        }

        rfftw_one(p, in, out);
        power_spectrum[0] = out[0]*out[0];  /* DC component */
        for (int i=1; i < (X_SIZE+1)/2; i++)  /* (k < N/2 rounded up) */
        {
            power_spectrum[i] = out[i]*out[i] + out[X_SIZE-i]*out[X_SIZE-i];
            //printf("i=%d power_spectrum[i]=%f\n", i, power_spectrum[i]);
        }

        if (X_SIZE % 2 == 0) /* N is even */
        {
            power_spectrum[X_SIZE/2] = out[X_SIZE/2]*out[X_SIZE/2];  /* Nyquist freq. */
        }

    // draw FFT
    glBegin(GL_LINES);
        for (int i=0; i<=X_SIZE/2; i++)
        {
            glColor4f(0,1,0,1);
            glVertex2i(i*2, 0);
            glVertex2i(i*2, power_spectrum[i]/FFT_SCALE_FACTOR);
            glVertex2i(i*2+1, 0);
            glVertex2i(i*2+1, power_spectrum[i]/FFT_SCALE_FACTOR);
        }
    glEnd();

    // draw points here
    glBegin(GL_POINTS);
        for (int i=0; i<X_SIZE; i++)
        {
            glColor4f(1,0,0,1);
            glVertex2f(i, running_buffer[i]*((float)SCREEN_HEIGHT/ADC_RESOLUTION));
        }
    glEnd();

    glBegin(GL_LINES);
        for (int i=1; i<X_SIZE; i++)
        {
            glColor4f(1,0,0,1);
            glVertex2f(i, running_buffer[i-1]*((float)SCREEN_HEIGHT/ADC_RESOLUTION));
            glVertex2f(i, running_buffer[i]*((float)SCREEN_HEIGHT/ADC_RESOLUTION));
        }
    glEnd();

    float delta=0, theta=0, alpha=0, beta=0, gamma=0, mu=0, total=0;

    for (int i=0; i<(X_SIZE/2); i++) total += power_spectrum[i];
    for (int i=0; i<4; i++) delta += power_spectrum[i];
    for (int i=4; i<=8; i++) theta += power_spectrum[i];
    for (int i=8; i<=13; i++) alpha += power_spectrum[i];
    for (int i=14; i<=30; i++) beta += power_spectrum[i];
    for (int i=30; i<=100; i++) gamma += power_spectrum[i];
    for (int i=8; i<=13; i++) mu += power_spectrum[i];
    delta /= total; theta /= total; alpha /= total; beta /= total; gamma /= total; mu /= total;

    // do BCI paddle control
    const float THETA_MIN=0.01f, THETA_MAX=0.04f;
    posy = (SCREEN_HEIGHT - 64) * (theta - THETA_MIN) / (THETA_MAX - THETA_MIN);
    // TODO: low-pass filter

/*    // 0.03 BELOW THAT MOVE PADDLE DOWN, OTHERWISE MOVE UP
    const float MU_THRESHOLD = 0.03f;
    if (mu < MU_THRESHOLD)
        paddle1yvel = 0.1f;
    else
        paddle1yvel = -0.1f; */

    // update and draw pong game state
	pongUpdateAndRender();

    // update and draw P300 state machine
    p300UpdateAndRender();

    // draw text
    glPrint(format("Delta=%f; Theta=%f; Alpha=%f", delta, theta, alpha), 10, 10, blue);
    glPrint(format("Beta=%f; Gamma=%f; Mu=%f", beta, gamma, mu), 10, 30, blue);

    SDL_GL_SwapBuffers();

    Frames++;
    {
        GLint t = SDL_GetTicks();
        if (t - T0 >= 5000)
        {
            GLfloat seconds = (t - T0) / 1000.0;
            GLfloat fps = Frames / seconds;
            printf("%d frames in %g seconds = %g FPS\n", Frames, seconds, fps);
            T0 = t;
            Frames = 0;
        }
    }

    return TRUE;
}

int main(int argc, char **argv)
{
    for (int i=0; i<X_SIZE; i++)
    {
        point_buffer[i] = ADC_RESOLUTION/2; //0;
        running_buffer[i] = ADC_RESOLUTION/2; //0;
    }

    int videoFlags;
    int done = FALSE;
    SDL_Event event;
    const SDL_VideoInfo *videoInfo;

    if (SDL_Init(SDL_INIT_VIDEO) < 0)
	{
	    fprintf(stderr, "Video initialization failed: %s\n", SDL_GetError());
	    Quit(1);
	}

#ifdef WIN32
	freopen("CON", "w", stdout);
    freopen("CON", "w", stderr);
#endif

    videoInfo = SDL_GetVideoInfo();

    if (!videoInfo)
	{
	    fprintf(stderr, "Video query failed: %s\n", SDL_GetError());
	    Quit(1);
	}

    videoFlags  = SDL_OPENGL;
    videoFlags |= SDL_GL_DOUBLEBUFFER;
    videoFlags |= SDL_HWPALETTE;
    videoFlags |= SDL_RESIZABLE;

    if (videoInfo->hw_available)
        videoFlags |= SDL_HWSURFACE;
    else
        videoFlags |= SDL_SWSURFACE;

    if (videoInfo->blit_hw)
        videoFlags |= SDL_HWACCEL;

    SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);

    surface = SDL_SetVideoMode(SCREEN_WIDTH, SCREEN_HEIGHT, SCREEN_BPP, videoFlags);

    if (!surface)
	{
	    fprintf(stderr,  "Video mode set failed: %s\n", SDL_GetError());
	    Quit(1);
	}

    initGL();

    resizeWindow(SCREEN_WIDTH, SCREEN_HEIGHT);

    if (TTF_Init() == -1) 
    {
        printf("Unable to initialize SDL_ttf: %s \n", TTF_GetError());
        Quit(1);
    }

    // init font
    InitFont();

    // open serial
    openSerial();

    // open log file
    logFile = fopen(LOG_FILENAME, "w");

    if (logFile == NULL)
    {
        fprintf(stderr, "gl_plot main(): Can't open log output file %s!\n", LOG_FILENAME);
        Quit(1);
    }

    // init FFT
    printf("Initializing FFTW 2...\n");
    p = rfftw_create_plan(X_SIZE, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
    printf("Done!\n");

    // init pong
    printf("Initializing pong...\n");
	srand((unsigned)time(0)); // Init Random Number Generator
    pongInit();

    // init P300
    p300init();

    while (!done)
	{
	    while (SDL_PollEvent(&event))
		{
		    switch(event.type)
			{
			case SDL_VIDEORESIZE:
			    surface = SDL_SetVideoMode(event.resize.w, event.resize.h, SCREEN_BPP, videoFlags);
			    if (!surface)
				{
				    fprintf(stderr, "Could not get a surface after resize: %s\n", SDL_GetError());
				    Quit(1);
				}
			    resizeWindow(event.resize.w, event.resize.h);
			    break;
			case SDL_KEYDOWN:
				handleKeyPress(&event.key.keysym);
				pongHandleKeyDown(&event.key.keysym);
				p300HandleKeyDown(&event.key.keysym);
				break;
			case SDL_KEYUP:
				pongHandleKeyUp(&event.key.keysym);
				p300HandleKeyUp(&event.key.keysym);
				break;
			case SDL_QUIT:
			    done = TRUE;
			    break;
			default:
			    break;
			}
		}

		drawGLScene();
	}

    Quit(0);

    return(0);
}
                 
Makefile
CXX = g++

CFLAGS := -g -I/usr/include/SDL -O2 -pipe -Wall -ansi
LDFLAGS := -lrfftw -lfftw -lm -lGLEW -lGL -lGLU -lSDL_image -lSDL_ttf `sdl-config --cflags --libs`

PROGRAM = gl_plot

SRC=\
	serial.cpp\
	font.cpp\
	main.cpp\
	debug.cpp\
	pong.cpp\
	p300.cpp\

OBJ = $(SRC:%.cpp=%.o)

.cpp.o :
	$(CXX) -c $< $(CFLAGS) -o $@

all:	$(PROGRAM)

$(PROGRAM):	$(OBJ)
	$(CXX) -o $(PROGRAM) $(OBJ) $(LDFLAGS)

.PHONY : clean
clean:
	rm -f $(OBJ)                 
                 
p300.cpp
#include <stdio.h>
#include <stdlib.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <time.h>
#include <string>
//#include <rfftw.h>
#include "font.h"
#include "debug.h"
#include "config.h"

#include "p300.h"

#define TRAINING_DATA_FILENAME "p300_data"
#define TESTING_DATA_FILENAME "p300_data.t"

#define P300_READY 0
#define P300_TRAINING 1
#define P300_TESTING 2

unsigned int p300_state = P300_READY;

#define NUM_COLORS 5
#define NUM_TRIALS 3
#define BUFFER_SIZE 512

const float color_choices[NUM_COLORS][3] = { {1,0,0}, {0,1,0}, {0,0,1}, {1,1,0}, {1,0,1} };
const std::string color_names[NUM_COLORS] = { "Red", "Green", "Blue", "Yellow", "Purple" };
unsigned int trialBuffer[NUM_TRIALS][NUM_COLORS][BUFFER_SIZE];
float targetBuffer[BUFFER_SIZE], nonTargetBuffer[BUFFER_SIZE], testBuffers[NUM_COLORS][BUFFER_SIZE];

int current_color = -1, trainingTarget = -1; // black
unsigned int bufferPtr = 0;
unsigned int current_trial = 0;
bool color_histogram[NUM_COLORS], classified_histogram[NUM_COLORS];
bool displayBuffers = false;

bool p300StartTrial(bool clear);
void p300trainSVM();

void p300setTrainingTarget()
{
    trainingTarget = 3;//rand() % NUM_COLORS;
}

void p300init()
{
    for (int i=0; i<NUM_COLORS; i++)
    {
        color_histogram[i] = classified_histogram[i] = false;
        for (int j=0; j<BUFFER_SIZE; j++)
        {
            testBuffers[i][j] = 0;
        }
    }

    for (int i=0; i<BUFFER_SIZE; i++)
    {
        targetBuffer[i] = nonTargetBuffer[i] = 0;
    }

    p300setTrainingTarget();
    p300trainSVM();
}

void p300trainSVM()
{
    printf("TODO: Training SVM!\n");

    // load binary data file, use it to train the SVM and get coefficients
}

void p300addTrainingExample()
{
    printf("Adding training example buffer!\n");

    // so we know trainingTarget, so iterate through all trials and add the buffers where training target is the color to targetBuffer and add the buffers where training target is NOT the color to non-target buffer, then divide targetBuffer by NUM_TRIALS (since each trial has exactly one target) and divide nonTargetBuffer by (NUM_TRIALS*(NUM_COLORS-1))

    // clear buffers first!
    for (int i=0; i<BUFFER_SIZE; i++)
    {
        targetBuffer[i] = nonTargetBuffer[i] = 0;
    }

    for (int i=0; i<NUM_TRIALS; i++)
    {
        for (int j=0; j<NUM_COLORS; j++)
        {
            for (int k=0; k<BUFFER_SIZE; k++)
            {
                if (j == trainingTarget)
                {
                    targetBuffer[k] += trialBuffer[i][j][k];
                } else {
                    nonTargetBuffer[k] += trialBuffer[i][j][k];
                }
            }
        }
    }

    for (int i=0; i<BUFFER_SIZE; i++)
    {
        targetBuffer[i] /= NUM_TRIALS;
        nonTargetBuffer[i] /= (NUM_TRIALS*(NUM_COLORS-1));
    }

    // then scale the data
    for (int i=0; i<BUFFER_SIZE; i++) // goes from 0 to ADC_RESOLUTION initially so subtract ADC_RESOLUTION/2 then divide by ADC_RESOLUTION/2 to get it in [-1, +1] range
    {
        targetBuffer[i] = (targetBuffer[i]-ADC_RESOLUTION/2)/(ADC_RESOLUTION/2);
        nonTargetBuffer[i] = (nonTargetBuffer[i]-ADC_RESOLUTION/2)/(ADC_RESOLUTION/2);
    }

    // then append it to a training SVM text file
    FILE *fp = fopen(TRAINING_DATA_FILENAME, "a");
    if (fp == NULL)
    {
        printf("p300addTrainingExample(): Unable to open training data file %s!\n", TRAINING_DATA_FILENAME);
        return;
    }

    // first do non-target buffer
    fprintf(fp, "-1 ");
    for (int i=0; i<BUFFER_SIZE; i++)
    {
        fprintf(fp, "%d:%f ", i+1, nonTargetBuffer[i]);
    }
    fprintf(fp, "\n");

    // now do target buffer
    fprintf(fp, "+1 ");
    for (int i=0; i<BUFFER_SIZE; i++)
    {
        fprintf(fp, "%d:%f ", i+1, targetBuffer[i]);
    }
    fprintf(fp, "\n");

    fclose(fp);

    // average the trials, append the trial to a binary data file, also call p300trainSVM() ?
    // TODO: update targetBuffer, nonTargetBuffer
}

void p300testAndReport()
{
    printf("Testing and reporting!\n");

    // so we have a bunch of colors with multiple trials each, so we need to go through each color and then each trial and add all the data for each color into the appropriate test buffer then scale all of the test buffers by NUM_TRIALS, then spit each one into the test data file (after scaling the data, of course)

    // clear buffers first!
    for (int i=0; i<NUM_COLORS; i++)
    {
        for (int j=0; j<BUFFER_SIZE; j++)
        {
            testBuffers[i][j] = 0;
        }
    }

    for (int i=0; i<NUM_TRIALS; i++)
    {
        for (int j=0; j<NUM_COLORS; j++)
        {
            for (int k=0; k<BUFFER_SIZE; k++)
            {
                testBuffers[j][k] += trialBuffer[i][j][k];
            }
        }
    }

    for (int i=0; i<NUM_COLORS; i++)
    {
        for (int j=0; j<BUFFER_SIZE; j++)
        {
            testBuffers[i][j] /= NUM_TRIALS;
        }
    }

    // then scale the data
    for (int i=0; i<NUM_COLORS; i++)
    {
        for (int j=0; j<BUFFER_SIZE; j++)
        {
            testBuffers[i][j] = (testBuffers[i][j]-ADC_RESOLUTION/2)/(ADC_RESOLUTION/2);
        }
    }

    // then write it to a testing SVM text file
    FILE *fp = fopen(TESTING_DATA_FILENAME, "w");
    if (fp == NULL)
    {
        printf("p300testAndReport(): Unable to open testing data file %s!\n", TESTING_DATA_FILENAME);
        return;
    }

    // do all test buffers with 0 for classification label
    for (int i=0; i<NUM_COLORS; i++)
    {
        fprintf(fp, "0 ");
        for (int j=0; j<BUFFER_SIZE; j++)
        {
            fprintf(fp, "%d:%f ", j+1, testBuffers[i][j]);
        }
        fprintf(fp, "\n");
    }

    fclose(fp);

    // so take the current trial and classify it using libSVM, then set the classification histogram
    // TODO: update targetBuffer, nonTargetBuffer (???)
}

void p300AddSample(unsigned int s)
{
    if (p300_state != P300_READY)
    {
        if (bufferPtr < BUFFER_SIZE) // still room left in the buffer
        {
            trialBuffer[current_trial][current_color][bufferPtr++] = s;
        } else { // buffer filled up
            // start a new trial unless we are finished then set p300_state to READY
            if (p300StartTrial(false)) // finished
            {
                // do SVM stuff here
                if (p300_state == P300_TRAINING)
                {
                    p300addTrainingExample();
                    p300setTrainingTarget(); // done training, set a new target
                } else if (p300_state == P300_TESTING)
                    p300testAndReport();
                else
                    printf("p300AddSample(): Unknown state %d\n", p300_state);

                // reset state
                p300_state = P300_READY;
                glClearColor(0.0f, 0.0f, 0.0f, 0.0f);

                // might as well reset the state variables just to be safe
                current_color = -1;
                bufferPtr = 0;
                current_trial = 0;
            }
        }
    }
}

// returns true if all finished
bool p300StartTrial(bool clear)
{
    // reset buffer pointer
    bufferPtr = 0;

    if (clear)
    {
        // set trial to zero and clear histogram
        current_trial = 0;
        for (int i=0; i<NUM_COLORS; i++) color_histogram[i] = false;
    }

    // check if done with trial
    bool increment_trial = true;
    for (int i=0; i<NUM_COLORS; i++)
    {
        if (color_histogram[i] == false)
        {
            increment_trial = false;
            break;
        }
    }

    // done with trial
    if (increment_trial)
    {
        if (current_trial < NUM_TRIALS - 1) // check if finished all together
            current_trial++;
        else
            return true;

        for (int i=0; i<NUM_COLORS; i++) color_histogram[i] = false;
    }

    // pick a random color, set it in the histograph
    bool picked = false;
    while (!picked)
    {
        unsigned int next_color = rand() % NUM_COLORS;
        if (!color_histogram[next_color] && next_color != current_color)
        {
            current_color = next_color;
            color_histogram[current_color] = true;
            picked = true;
        }
    }

    glClearColor(color_choices[current_color][0], color_choices[current_color][1], color_choices[current_color][2], 0.0f);
    return false;
}

void p300UpdateAndRender()
{
    if (displayBuffers) // TODO: test and fix!!!
    {
        glBegin(GL_LINES);
            for (int i=1; i<BUFFER_SIZE; i++)
            {
                glColor4f(0,1,0,1);
                glVertex2f(i, targetBuffer[i-1]/2.0f*((float)SCREEN_HEIGHT/ADC_RESOLUTION));
                glVertex2f(i, targetBuffer[i]/2.0f*((float)SCREEN_HEIGHT/ADC_RESOLUTION));
            }
        glEnd();

        glBegin(GL_LINES);
            for (int i=1; i<BUFFER_SIZE; i++)
            {
                glColor4f(0,0,1,1);
                glVertex2f(i, nonTargetBuffer[i-1]/2.0f*((float)SCREEN_HEIGHT/ADC_RESOLUTION)+SCREEN_HEIGHT/2.0f);
                glVertex2f(i, nonTargetBuffer[i]/2.0f*((float)SCREEN_HEIGHT/ADC_RESOLUTION)+SCREEN_HEIGHT/2.0f);
            }
        glEnd();
    }

    // draw status text:
    std::string classified_str = "Classified: { ";
    // Classified: { Red Blue }
    for (int i=0; i<NUM_COLORS; i++)
    {
        if (classified_histogram[i])
        {
            classified_str += color_names[i] + " ";
        }
    }
    classified_str += " }";

    // add mode string to status text
    switch (p300_state)
    {
    case P300_READY:
        classified_str = "P300 Ready: " + classified_str;
        break;
    case P300_TRAINING:
        classified_str = "P300 Training Mode: " + classified_str;
        break;
    case P300_TESTING:
        classified_str = "P300 Testing Mode: " + classified_str;
        break;
    default:
        classified_str = "P300 Unknown State! " + classified_str;
    }

    // display training target too
    classified_str += " - Training Target: " + color_names[trainingTarget];

    glPrint(classified_str.c_str(), 10, 50, blue);
}

void p300HandleKeyDown(SDL_keysym *keysym)
{
	switch(keysym->sym)
	{
    case SDLK_F2: // switch to "training mode"
        if (p300_state == P300_READY)
        {
            p300_state = P300_TRAINING;
            p300StartTrial(true);
        }

        break;
    case SDLK_F3: // switch to "testing mode"
        if (p300_state == P300_READY)
        {
            p300_state = P300_TESTING;
            p300StartTrial(true);
        }

        break;
    case SDLK_F5:
        displayBuffers = !displayBuffers;
        break;
    default:
        break;
	}
}

void p300HandleKeyUp(SDL_keysym *keysym)
{
/*	switch(keysym->sym)
	{
	case SDLK_UP:
		paddle2yvel = 0;
		break;
	case SDLK_DOWN:
		paddle2yvel = 0;
		break;
	case SDLK_a:
		paddle1yvel = 0;
		break;
	case SDLK_z:
		paddle1yvel = 0;
		break;
	} */
}
                 
p300.h
#ifndef P300_H_
#define P300_H_

#include "SDL.h"

// state variables
extern unsigned int p300_state;

// methods
void p300UpdateAndRender();
void p300HandleKeyUp(SDL_keysym *keysym);
void p300HandleKeyDown(SDL_keysym *keysym);
void p300AddSample(unsigned int s);
void p300init();

#endif
                 
pong.cpp
#include <SDL/SDL.h>
#include <GL/gl.h>
#include <GL/glu.h>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <ctime>
#include <malloc.h>
#include <string>
#include "config.h"
#include "glyphs.h"

void drawsprite(int x, int y, float r, float g, float b);
void drawpaddlesprite(int x, int y, float r, float g, float b);
void updatepaddle1();
void updatepaddle2();
void drawglyph(int num, int x, int y, float r, float g, float b);
void drawline();
void drawscore();
void ResetVelocity();
int randomNum();

#define PADDLE_WIDTH 16
#define PADDLE_HEIGHT 128

unsigned int score1=0, score2=0;
float posy=SCREEN_HEIGHT/2, paddle1yvel=0;
float pos2y=SCREEN_HEIGHT/2, paddle2yvel=0;

float ballposx=SCREEN_WIDTH/2-8, ballposy=SCREEN_HEIGHT/2-8;
float ballvelx=0, ballvely=0;

const float vel_threshold = 0.1f;

const unsigned char sprite[] =
{
0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,
0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,
0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,
0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0
};

void updateball()
{
	ballposx += ballvelx;
	ballposy += ballvely;

	if (ballposy <= 0)
	{
		ballvely = -ballvely;
		ballposy = 0;
	} else if (ballposy >= SCREEN_HEIGHT - 16) {
		ballvely = -ballvely;
		ballposy = SCREEN_HEIGHT - 16;
	}

	if (ballposx <= 0)
	{
		ballvelx = -ballvelx;
		ballposx = 0;
	} else if (ballposx >= SCREEN_WIDTH - 16) {
		ballvelx = -ballvelx;
		ballposx = SCREEN_WIDTH - 16;
	}

	// COLLISION DETECTION
	// front checking for both paddles
	if ((ballposx >= 75) && (ballposx <= 75+PADDLE_WIDTH) && (ballposy >= posy-(PADDLE_HEIGHT-1)) && (ballposy <= posy+PADDLE_HEIGHT-1))
	{
		ballposx = 75+(PADDLE_WIDTH-1);
		ballvelx = -ballvelx;

		if (ballvelx > 1)
			ballvelx += vel_threshold;
		else if (ballvelx < 1)
			ballvelx -= vel_threshold;

		if (ballvely > 1)
			ballvely += vel_threshold;
		else if (ballvely < 1)
			ballvely -= vel_threshold;
	}

	if ((ballposx <= (SCREEN_WIDTH-91)) && (ballposx >= (SCREEN_WIDTH-91)-PADDLE_WIDTH) && (ballposy >= pos2y-(PADDLE_WIDTH-1)) && (ballposy <= pos2y+PADDLE_HEIGHT-1))
	{
		ballposx = (SCREEN_WIDTH-91)-PADDLE_WIDTH;
		ballvelx = -ballvelx;

		if (ballvelx > 1)
			ballvelx += vel_threshold;
		else if (ballvelx < 1)
			ballvelx -= vel_threshold;

		if (ballvely > 1)
			ballvely += vel_threshold;
		else if (ballvely < 1)
			ballvely -= vel_threshold;
	}

	// back checking for both paddles
	if ((ballposx == 75-PADDLE_WIDTH) && (ballposy >= posy-(PADDLE_WIDTH-1)) && (ballposy <= posy+PADDLE_HEIGHT-1))
		ballvelx = -ballvelx;

	if ((ballposx == (SCREEN_WIDTH-91)+PADDLE_WIDTH) && (ballposy >= pos2y-(PADDLE_WIDTH-1)) && (ballposy <= pos2y+PADDLE_HEIGHT-1))
		ballvelx = -ballvelx;

	// Top and bottom checking for paddle1
	if ((ballposy == posy-PADDLE_WIDTH) && (ballposx > 75-(PADDLE_WIDTH-1)) && (ballposx < 75+(PADDLE_WIDTH-1)))
		ballvely = -ballvely;

	if ((ballposy == posy+PADDLE_HEIGHT) && (ballposx > 75-(PADDLE_WIDTH-1)) && (ballposx < 75+(PADDLE_WIDTH-1)))
		ballvely = -ballvely;

	// Top and bottom checking for paddle 2
	if ((ballposy == pos2y-PADDLE_WIDTH) && (ballposx > (SCREEN_WIDTH-91)-(PADDLE_WIDTH-1)) && (ballposx < (SCREEN_WIDTH-91)+(PADDLE_WIDTH-1)))
		ballvely = -ballvely;

	if ((ballposy == pos2y+PADDLE_HEIGHT) && (ballposx > (SCREEN_WIDTH-91)-(PADDLE_WIDTH-1)) && (ballposx < (SCREEN_WIDTH-91)+(PADDLE_WIDTH-1)))
		ballvely = -ballvely;

	if (ballposx == 0)
	{
		score1++;

		ResetVelocity();
	}

	if (ballposx == SCREEN_WIDTH-16)
	{
		score2++;

		ResetVelocity();
	}
}

void ResetVelocity()
{
	if (ballvelx > 1)
		ballvelx = 1;
	else if (ballvelx < 1)
		ballvelx = -1;

	if (ballvely > 1)
		ballvely = 1;
	else if (ballvely < 1)
		ballvely = -1;
}

void pongUpdateAndRender()
{
	int tick = SDL_GetTicks();

	updateball();
	updatepaddle1();
	updatepaddle2();

	drawsprite(ballposx, ballposy, 0, 1, 0);

	drawpaddlesprite(75, posy, 1, 0, 0);
	drawpaddlesprite((SCREEN_WIDTH-91), pos2y, 0, 0, 1);

	drawscore();
	drawline();
}

void drawscore()
{
	drawglyph(score1, SCREEN_WIDTH/2+10, 25, 1, 1, 0);
	drawglyph(score2, SCREEN_WIDTH/2-12-10, 25, 1, 1, 0);
}

void drawglyph(int num, int x, int y, float r, float g, float b)
{
	const unsigned char *glyph;

	switch (num)
	{
	case 0:
		glyph = num0;
		break;
	case 1:
		glyph = num1;
		break;
	case 2:
		glyph = num2;
		break;
	case 3:
		glyph = num3;
		break;
	case 4:
		glyph = num4;
		break;
	case 5:
		glyph = num5;
		break;
	case 6:
		glyph = num6;
		break;
	case 7:
		glyph = num7;
		break;
	case 8:
		glyph = num8;
		break;
	case 9:
		glyph = num9;
		break;
	default:
		glyph = num0;
		break;
	}

	glMatrixMode(GL_PROJECTION);
	glPushMatrix();
	glLoadIdentity();
	gluOrtho2D(0, SCREEN_WIDTH, SCREEN_HEIGHT, 0);

	glColor4f(r,g,b,1.0f);
	glBegin(GL_POINTS);
		for (int i = 0, c = 0; i < 20; i++)
		{
			for (int j = 0; j < 12; j++, c++)
			{
				if (glyph[c])
				{
					glVertex2i(x+j,y+i);
				}
			}
		}
	glPopMatrix();
}

void drawsprite(int x, int y, float r, float g, float b)
{
	glBegin(GL_POINTS);
	glColor4f(r,g,b,1.0f);

	for (int i = 0, c = 0; i < 16; i++)
	{
		for (int j = 0; j < 16; j++, c++)
		{
			if (sprite[c])
			{
				glVertex2i(x+j,y+i);
			}
		}
	}
	glEnd();
}

#define LINE_SPACING 5
void drawline()
{
	glBegin(GL_LINES);

	glColor4f(1.0f, 1.0f, 1.0f, 1.0f); // white divider line in the middle
	for (int i=0; i<SCREEN_HEIGHT; i+=LINE_SPACING*2)
	{
		glVertex2i(SCREEN_WIDTH/2, i);
		glVertex2i(SCREEN_WIDTH/2, i+LINE_SPACING);
	}

	glEnd();
}

void pongInit()
{
	printf("SCREEN_HEIGHT=%d\n", SCREEN_HEIGHT);
	int i = randomNum();

	if (i == 0)
		ballvelx = -1;
	else
		ballvelx = 1;

	i = randomNum();

	if (i == 0)
		ballvely = -1;
	else
		ballvely = 1;
}

void drawpaddlesprite(int x, int y, float r, float g, float b)
{
	glBegin(GL_POINTS);
	glColor4f(r,g,b,1.0f);
	for (int i = 0, c = 0; i < PADDLE_HEIGHT; i++)
	{
		for (int j = 0; j < PADDLE_WIDTH; j++, c++)
		{
			glVertex2i(x+j,y+i);
		}
	}
	glEnd();
}

void pongHandleKeyDown(SDL_keysym *keysym)
{
	switch(keysym->sym)
	{
	case SDLK_UP:
		paddle2yvel = -1;
		break;
	case SDLK_DOWN:
		paddle2yvel = 1;
		break;
	/*case SDLK_a:
		paddle1yvel = -1;
		break;
	case SDLK_z:
		paddle1yvel = 1;
		break; */
	default:
		break;
	}
}

void pongHandleKeyUp(SDL_keysym *keysym)
{
	switch(keysym->sym)
	{
	case SDLK_UP:
		paddle2yvel = 0;
		break;
	case SDLK_DOWN:
		paddle2yvel = 0;
		break;
/*	case SDLK_a:
		paddle1yvel = 0;
		break;
	case SDLK_z:
		paddle1yvel = 0;
		break; */
    default:
        break;
	}
}

void updatepaddle1()
{
    //printf("updating...\n");
	posy += paddle1yvel;

	if (posy < 0)
		posy = 0;
	else if (posy > (SCREEN_HEIGHT - 64))
		posy = SCREEN_HEIGHT - 64;
}

void updatepaddle2()
{
	pos2y += paddle2yvel;

	if (pos2y < 0)
		pos2y = 0;
	else if (pos2y > (SCREEN_HEIGHT - 64))
		pos2y = SCREEN_HEIGHT - 64;
}

int randomNum() // Obtain a random integer between defined range
{
	int range=(1-0)+1; // Calculate range

	// Use rand()
	int retval=0+int(range*rand()/(RAND_MAX + 1.0));

	return retval; // Return the number
}
                 
pong.h
#ifndef PONG_H_
#define PONG_H_

void pongInit();
void pongUpdateAndRender();

void pongHandleKeyDown(SDL_keysym *keysym);
void pongHandleKeyUp(SDL_keysym *keysym);

extern float posy, paddle1yvel;

#endif
                 
serial.cpp
#include <stdio.h>   // Standard input/output definitions
#include <stdlib.h>
#include <string.h>  // String function definitions
#include <unistd.h>  // UNIX standard function definitions
#include <fcntl.h>   // File control definitions
#include <errno.h>   // Error number definitions
#include <termios.h> // POSIX terminal control definitions

#include "serial.h"

#define BUFFER_SIZE 80

char lastError[1024], serialBuffer[BUFFER_SIZE];

int usbdev = 0;

#define PORT_NAME "/dev/ttyUSB0"
#define BAUD_RATE B57600

//#define NULL_SERIAL

bool openSerial()
{
#ifndef NULL_SERIAL
    system("stty -F /dev/ttyUSB0 57600 cs8 -cstopb -parity -icanon min 1 time 1");
    usbdev = open("/dev/ttyUSB0", O_RDWR);
    return (usbdev != NULL);
#else
    return true;
#endif
}

char* readByte()
{
    return NULL;
}

unsigned int readSerialValue()
{
#ifndef NULL_SERIAL
    int ptr=0;
    unsigned char last_read = NULL;
    while (ptr < BUFFER_SIZE - 1 && last_read != '\n')
    {
        read(usbdev, &serialBuffer[ptr], 1);
        last_read = serialBuffer[ptr];
        ptr++;
    }

    int val = 0;
    sscanf(serialBuffer, "%d\n", &val);
    return val;
#else
    return 0;
#endif
}

void closeSerial()
{
#ifndef NULL_SERIAL
    close(usbdev);
#endif
}
                 
serial.h
#ifndef SERIAL_H_INCLUDED
#define SERIAL_H_INCLUDED

bool openSerial();
void closeSerial();

unsigned int readSerialValue();

#endif // SERIAL_H_INCLUDED
                 

MATLAB Plotting Source Code


plot_samples.m
% settings
SerialPort='com3'; % serial port
N=200;
Fs=200;

m=zeros(1,N);

s = serial(SerialPort);
set(s,'BaudRate',57600);
fopen(s);

for i = 1:N
    datum = fscanf(s, '%s');
    fprintf('%s\n', datum);

    if (length(datum) > 0)
        m(i) = str2num(datum);
    else
        m(i) = 0;
    end
end

% Clean up the serial port
fclose(s);
delete(s);
clear s;

% Filter m with 60 Hz notch
Wo = 60/(Fs/2);  BW = Wo/35;
[b,a] = iirnotch(Wo,BW);
m = filter(b,a,m);

% Remove DC offset
mu = mean(m);
m = m - mu + 1024/2;

figure(1)
hLine = plot(m);
ylim([0 1024]);
set(hLine,'YData',m);

figure(2)
L=length(m);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(m,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
                 
plot_samples_rt.m
% settings
SerialPort='com3'; % serial port
N=200;
Fs=200;

m=zeros(1,N);

s = serial(SerialPort);
set(s,'BaudRate',57600);
fopen(s);

for i = 1:N
    datum = fscanf(s, '%s');
    fprintf('%s\n', datum);

    if (length(datum) > 0)
        m(i) = str2num(datum);
    else
        m(i) = 0;
    end
end

% Clean up the serial port
fclose(s);
delete(s);
clear s;

% Filter m with 60 Hz notch
Wo = 60/(Fs/2);  BW = Wo/35;
[b,a] = iirnotch(Wo,BW);
m = filter(b,a,m);

% Remove DC offset
mu = mean(m);
m = m - mu + 1024/2;

figure(1)
hLine = plot(m);
ylim([0 1024]);
set(hLine,'YData',m);

figure(2)
L=length(m);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(m,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|') 
                 
serial_test.m
SerialPort='com8'; % serial port
N=1000;

m=zeros(1,1000);
figure(1)
hLine = plot(m);
ylim([0 1024]);

KeepRunning = 1;
while KeepRunning
    s = serial(SerialPort);
    set(s,'BaudRate',57600);
    fopen(s);
    for i = 1:N
        datum = fscanf(s, '%s');
        %fprintf('%s\n', datum);

        if (length(datum) > 0)
            m(i) = str2num(datum);
        else
            m(i) = 0;
        end

        if (ishandle(1))
            set(hLine,'YData',m);
        else
            KeepRunning = 0;
            break;
        end
        drawnow
    end
    
    % Clean up the serial port
    fclose(s);
    delete(s);
    clear s;
end
                 

Schematics



Figure: Project Schematics


Cost Details and Parts


Power Supply Circuit

QuantityPart NumberPart NameCost
1SparkFun: COM-00102SPDT Mini Power Switch$1.50
1LABGreen LEDLAB
2LAB1 MΩ ResistorLAB
3LAB1 kΩ ResistorLAB
1Pre-Owned100 Ω ResistorPre-Owned (Resistor Kit)
2LABJumper Wire2*1.00=$2.00
1Digi-Key: BC4AAW-NDBattery Holder 4-AA Cells$1.33
4Panasonic Akaline PlusAA BatteryPre-Owned
UART Opto-Isolation

QuantityPart NumberPart NameCost
1Digi-Key: 160-1791-NDFairchild Semiconductor 6N137 Opto-Isolator$1.00
1LAB0.01 μF CapacitorLAB
1Pre-Owned220 Ω ResistorPre-Owned (Resistor Kit)
1Pre-Owned1.1 kΩ ResistorPre-Owned (Resistor Kit)
1LABJumper Wire$1.00
1LABDIP Socket$0.50
Amplifier Board

QuantityPart NumberPart NameCost
1Digi-Key: AD620ANZ-NDIC AMP INST LP LN 18MA 8DIP AD620 Instrumentation Amplifier$8.27
2Digi-Key: CA3140EZ-NDIC OP AMP 4.5MHZ BIMOS 8-DIP 3140 Operational Amplifier2*1.89=$3.78
1Digi-Key: 3386F-105LF-ND1 MΩ Trim Potentiometer$1.12
2LAB0.01 μF CapacitorLAB
1Digi-Key: 490-5362-ND0.0068 μF Capacitor$0.28
1Digi-Key: 445-2851-ND1 μF Capacitor$0.38
2Pre-Owned.22 μF CapacitorPre-Owned (Capacitor Kit)
2Pre-Owned2.2 kΩ ResistorPre-Owned (Resistor Kit)
3LAB1 MΩ ResistorLAB
1LAB10 kΩ ResistorLAB
2Pre-Owned15 kΩ ResistorPre-Owned (Resistor Kit)
Misc/Other

QuantityPart NumberPart NameCost
1LABMega644$6.00
1LABFTDI USB Comms Board$4.00
1LABTarget Board$4.00
1LABTarget Board SocketLAB
7LABFemale Wire Headers7*0.05=$0.35
7LABFemale Wire Sockets7*0.05=$0.35
1LABWhite Board$6.00
1LABSolder Board (6 inch)$2.50
1Contec Medical Systems: eBayTen Pcs Silver Plating Electrodes$20.00
1Pre-OwnedBaseball CapPre-Owned

The total cost of the project was $64.36 which is rather costly (owing to the medical-grade electrodes that we used) but still well under our budget.


Tasks and Workload Division


The tasks were divided among ourselves as follows:

Charles
Software development
Wrote the MATLAB code, AVR firmware, OpenGL C code
BCI research, project idea
Mengxiang
Hardware design, construction, circuit debugging
Part selection for hardware (Digi-Key ordering)
Testing the amplifier board

Overall, we felt this was a fair, equal division of labor.


References


Background sites


Electroencephalography (EEG) article, Wikipedia.
http://en.wikipedia.org/wiki/Electroencephalography

LibSVM -- A Library for Support Vector Machines
http://www.csie.ntu.edu.tw/~cjlin/libsvm/

Simple Direct Media Layer (SDL)
http://www.libsdl.org/

FFTW
http://www.fftw.org/

MATLAB Documentation
http://www.mathworks.com/help/techdoc/

OpenEEG Project
http://openeeg.sourceforge.net/


Papers


Matsuoka, G. and Sugi, T. and Kawana, F. and Nakamura, M. Automatic detection of apnea and EEG arousals for sleep apnea syndrome.
In ICCAS-SICE, 2009, pages 4651-4654.

F. Lotte. A review of classification algorithms for EEG-based brain computer interfaces.
In Journal of Neural Engineering, 2007.

Hazrati MKh and Efranian A. An online EEG-based brain-computer interface for controlling hand grasp using an adaptive probabilisitc neural network.
In PubMed, 2010.

Jorge Baztarrica Ochoa. EEG Signal Classification for Brain Computer Interface Applications.
Thesis, École Polytechnique Federale de Lausanne, 2002.

Kouhyar Tavakolian, Faratash Vasefi, Kaveh Naziripour, and Siamak Reazei. Mental task classification for brain computer interface applications.
In First Canadian Student Conference on Biomedical Programming.

Yuanqing Li, Chuanchu Wang, Haihong Zhang, and Cuntai Guan. An EEG-based BCI System for 2D Cursor Control.

Ali S. AlMejrad. Human Emotions Detection using Brain Wave Signals.
In European Journal of Scientific Research, 2010, pages 640-659.

Akinari Onishi, Yu Zhang, Qibin Zhao, Andrzej Cichocki. Fast and Reliable P300-Based BCI with Facial Images.

Christoph Guger, et al. Rapid Prototyping of an EEG-based Brain-Computer Interface (BCI).
University of Technology Graz, Austria.

Schoresch Presentation Slides: Neurofeedback - applied neuroscience
Kompetenzzentrum für Neurofeedback.

Hideaki Touyama. EEG-Based Personal Identification
Toyama Perfectural University, Japan.

Chih-Wei Hsu, Chih-Chung Chang, and Chih-Jen Lin. A Practical Guide to Support Vector Classification
Department of Computer Science, National Taiwan University, Taiwan.

Muhammad Bilal Khalid, et al. Think. Done! A Brain Computer Interface (BCI)
Thesis, National University of Sciences and Technology, Rawalpindi.

Kana Omori, Tomonari Yamaguchi, and Katsuhiro Inoue. Feature Extraction from EEG Signals in P300 Spelling System
ICROS-SICE International Joint Conferent 2009, Fukuoka International Congress Center, Japan.

Dandan Huang, et al. Decoding human motor activity from EEG single trials for a discrete two-dimensional cursor control.
IOP Publishing, Journal of Neural Engineering, 2009.

Ulrich Hoffmann, et al. An efficient P300-based brain-computer interface for disabled subjects.
Elsevier, Journal of Neuroscience Methods 2008, pages 115-125.


Code/designs borrowed from others


Cornell ECE 4760 Lab 2 ADC Sleep Example
http://people.ece.cornell.edu/land/courses/ece4760/labs/s2012/lab2.html

chipstein: Homebrew DIY EEG, EKG, and EMG
https://sites.google.com/site/chipstein/home-page/eeg-with-an-arduino

Joerg Wunsch's UART AVR C Code
www.nongnu.org/avr-libc/examples/stdiodemo/uart.c

MATLAB Notch Filter Implementation Example
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292960

NeHe Productions: OpenGL Tutorials (Template code for OpenGL font rendering and "OpenGL with SDL" example)
http://nehe.gamedev.net/

flipcode - Safe sprintf Example
http://www.flipcode.com/archives/Safe_sprintf.shtml


Data sheets


Analog Devices AD 620 Low Cost Low Power Instrumentation Amplifier
http://www.analog.com/static/imported-files/data_sheets/AD620.pdf

Fairchild Semiconductor Single-Channel: 6N137 High Speed 10MBit/s Logic Gate Optocoupler
http://www.fairchildsemi.com/ds/6N/6N137.pdf

Intersil CA3140, CA3140A: 4.5MHz, BiMOS Operational Amplifier with MOSFET Input/Bipolar Output
http://www.intersil.com/data/fn/fn957.pdf

Atmel ATmega644/V: 8-bit Atmel Microcontroller with 64K Bytes In-System Programmable Flash
http://www.atmel.com/Images/doc2593.pdf


Vendor sites


Contec Medical Systems Co., LTD
http://www.contecmed.com/main/Default.asp

Analog Devices
http://www.analog.com/en/index.html

Digi-Key
http://www.digikey.com/

SparkFun Electronics
http://www.sparkfun.com/

Atmel
http://www.atmel.com/

Texas Instruments
http://www.ti.com/

Fairchild Semiconductor
http://www.fairchildsemi.com/

Maxim IC
http://www.maxim-ic.com/

Intersil
http://www.intersil.com/cda/home/



Contact Us

You may reach Charles Moyes via e-mail at cwm55@cornell.edu and Mengxiang Jiang via mj294@cornell.edu