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