Search…
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:
1
% Import raw CSV file
2
input_data = csvread('EPOC+datafile.csv',2,0); % This excludes the first two rows
3
eegcols = 3:16; % EEG Columns.
4
eeg.raw = input_data(:, eegcols);
Copied!
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).
1
% median value removal
2
med = median(eeg.raw,2); % calculate median of each sample
3
eeg.raw = eeg.raw - repmat(med, 1, 14); % remove it
4
5
% limit slew rate
6
for j=2:size(eeg.raw,1)
7
del = eeg.raw(j,:) - eeg.raw(j-1,:);
8
del = min(del, ones(1,14)*15);
9
del = max(del, -ones(1,14)*15);
10
eeg.raw(j,:) = eeg.raw(j-1,:) + del;
11
end
12
13
% High pass filter
14
a = 0.06; % HPF filter coeffs
15
b = 0.94;
16
17
preVal = zeros(1,14);
18
eeg.filt = zeros(size(eeg.raw));
19
for j=2:size(eeg.raw,1)
20
preVal = a * eeg.raw(j,:) + b * preVal;
21
eeg.firfilt(j,:) = eeg.raw(j,:) - preVal;
22
end
Copied!
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():
1
IIR_TC = 256; % 2 second time constant- adjust as required
2
EEG_data = input_data( : ,3:16); % import raw data
3
[rows, columns] = size(EEG_data); % rows = number of data samples, columns = 14
4
AC_EEG_data = zeros(rows, columns); % create an empty array
5
back = EEG_data( 1, : ); % copy first row into background
6
7
% run IIR filter
8
for r = 2 : rows
9
back = (back * ( IIR_TC- 1 ) + EEG_data( r,:)) / IIR_TC;
10
eeg.iirfilt(r,:) = EEG_data( r,:) - back;
11
end
Copied!
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:
1
% Variables
2
fftlength = 1024; % Set FFT length
3
stepsize = 32; % Set step size
4
samples = fftlength:stepsize:length(eeg.raw); % create an index array
5
c = 1:fftlength;
6
channel = 1; % AF3
7
8
%Specify hanning window
9
hanning = [1:fftlength]';
10
hanning_in = 2* pi() * (hanning - (fftlength+1)/2)/(fftlength+1);
11
%rescaled x-axis to match sample length?
12
hanning = (sin(hanning_in)./hanning_in).^2; % sinc^2?
13
hanning = repmat(hanning, 1, size(eeg.raw,2)); % match to number of channels
14
f=[128/fftlength:128/fftlength:128]; % frequency index for the spectral array
15
16
for kk = 1:length(samples)
17
18
k = samples(kk);
19
% step through every quarter second starting at first possible sample
20
spectrum = fft(eeg.iirfilt(k-fftlength+1:k,:) .* hanning)/fftlength; % apply window to HP filtered data
21
spectrum = 2 * (sqrt(spectrum .* conj(spectrum))); % get magnitude
22
23
%plotting the time domain and IIR FFT data
24
subplot(211)
25
% time domain
26
plot(c,eeg.iirfilt(k-fftlength+1:k, channel));
27
title('Amplitude of Channels');
28
xlabel('Samples');
29
ylabel('Amplitude');
30
xlim([0 fftlength]);
31
%IIR fft
32
subplot(212)
33
plot(f,spectrum(:, channel));
34
%add graph details.
35
title('Amplitude Spectrum of Channels');
36
xlabel('Frequency (Hz)');
37
ylabel('|Y(f)|');
38
xlim([0 50]);
39
% ylim([0 80]);
40
pause(0.5);
41
42
end
Copied!
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
Last modified 2yr ago
Export as PDF
Copy link