Comment on page

# Code Examples

We have provided a couple of code samples below for matlab to get you started. Firstly you need to export the data from EmotivPRO as CSV. The data output format in CSV is based around the EDF export format, it contains one row for each datapoint recorded. The first row of the file contains information on file name (A1), when the file was recorded(A2), sample rate (A3) and the columns headers (A5). If you expand A5 using text to columns in excel this will help identify which are important to you. For this example we are only interested in EEG data in columns 3:16, the column mapping is shown in the table below for reference:

Channel | Column |

3 | AF3 |

4 | F7 |

5 | F3 |

6 | FC5 |

7 | T7 |

8 | P7 |

9 | O1 |

10 | O2 |

11 | P8 |

12 | T8 |

13 | FC6 |

14 | F4 |

15 | F8 |

16 | AF4 |

To import the data into matlab we use the following command:

% Import raw CSV file

input_data = csvread('EPOC+datafile.csv',2,0); % This excludes the first two rows

eegcols = 3:16; % EEG Columns.

eeg.raw = input_data(:, eegcols);

The above matlab code imports the csv data and places only the EEG data into the eeg struct;

To further analyse your data, you can convert it from the time domain to the frequency domain using an FFT, but before performing an Fast Fourier Transform (FFT) it is necessary to remove the DC offset from the data.

The simplest removal method and least accurate is to subtract the average value from the entire data channel. This method can result in increased low frequency noise. Ideally you should also apply a high-pass filter which matches the characteristics of the electronics - that is, you should use a 0.16Hz first order high-pass filter to remove the background signal (this also removes any longer term drift, which is not achieved by the average subtraction method).

% median value removal

med = median(eeg.raw,2); % calculate median of each sample

eeg.raw = eeg.raw - repmat(med, 1, 14); % remove it

% limit slew rate

for j=2:size(eeg.raw,1)

del = eeg.raw(j,:) - eeg.raw(j-1,:);

del = min(del, ones(1,14)*15);

del = max(del, -ones(1,14)*15);

eeg.raw(j,:) = eeg.raw(j-1,:) + del;

end

% High pass filter

a = 0.06; % HPF filter coeffs

b = 0.94;

preVal = zeros(1,14);

eeg.filt = zeros(size(eeg.raw));

for j=2:size(eeg.raw,1)

preVal = a * eeg.raw(j,:) + b * preVal;

eeg.firfilt(j,:) = eeg.raw(j,:) - preVal;

end

Another method is to use an IIR filter to track the background level and subtract it, as is shown below.

In Matlab pseudocode, assuming the first row has been removed from the array input_data():

IIR_TC = 256; % 2 second time constant- adjust as required

EEG_data = input_data( : ,3:16); % import raw data

[rows, columns] = size(EEG_data); % rows = number of data samples, columns = 14

AC_EEG_data = zeros(rows, columns); % create an empty array

back = EEG_data( 1, : ); % copy first row into background

% run IIR filter

for r = 2 : rows

back = (back * ( IIR_TC- 1 ) + EEG_data( r,:)) / IIR_TC;

eeg.iirfilt(r,:) = EEG_data( r,:) - back;

end

This demonstration code is not efficient in memory and assumes the entire file is available. It is quite straightforward to modify the code to replace the data in the source array rather than making a separate AC-coupled array, and also to run the IIR filter in open-ended form for processing in real time.

Note the vectorised form of the background recalculation at each iteration - each individual channel background average is retained in the relevant column of “back”. At each step the running average is re-estimated using the new input value. Note also that the first IIR _TC samples are biased towards the initial value but this settles down after about 2 * IIR_TC samples.

It is very important to remove the background signal before performing an FFT - you should also apply a tapered window function such as a HANNING transform before executing the FFT to ensure there are no wrapping artefacts where the FFT treats the data as an infinitely repeating sequence, and any mismatch between the first and last samples appears as a STEP FUNCTION in the analysis, injecting noise across the spectrum.

The following snippet explains how to run the FFT and then display the result in a graph:

% Variables

fftlength = 1024; % Set FFT length

stepsize = 32; % Set step size

samples = fftlength:stepsize:length(eeg.raw); % create an index array

c = 1:fftlength;

channel = 1; % AF3

%Specify hanning window

hanning = [1:fftlength]';

hanning_in = 2* pi() * (hanning - (fftlength+1)/2)/(fftlength+1);

%rescaled x-axis to match sample length?

hanning = (sin(hanning_in)./hanning_in).^2; % sinc^2?

hanning = repmat(hanning, 1, size(eeg.raw,2)); % match to number of channels

f=[128/fftlength:128/fftlength:128]; % frequency index for the spectral array

for kk = 1:length(samples)

k = samples(kk);

% step through every quarter second starting at first possible sample

spectrum = fft(eeg.iirfilt(k-fftlength+1:k,:) .* hanning)/fftlength; % apply window to HP filtered data

spectrum = 2 * (sqrt(spectrum .* conj(spectrum))); % get magnitude

%plotting the time domain and IIR FFT data

subplot(211)

% time domain

plot(c,eeg.iirfilt(k-fftlength+1:k, channel));

title('Amplitude of Channels');

xlabel('Samples');

ylabel('Amplitude');

xlim([0 fftlength]);

%IIR fft

subplot(212)

plot(f,spectrum(:, channel));

%add graph details.

title('Amplitude Spectrum of Channels');

xlabel('Frequency (Hz)');

ylabel('|Y(f)|');

xlim([0 50]);

% ylim([0 80]);

pause(0.5);

end

The FFT is carried out on the last fftlength samples with the window moving forward every stepsize. In the plot shown we used in the above analysis we looked at channel AF3 which is column 1 if our data eeg.iirfilt. The value displayed is set in the variable channel. It would be easy to modify the above file to display the FIR filtered results