EEG_Pipeline/Emotiv_EPOC_GMAC_script.m at master · mdtdev/EEG_Pipeline · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
% Script for extracting power in each frequency range from EEG data block with
% rows being the time course and columns being the electrode axis. Have not
% checked to see if this works yet.
%
% From: https://emotiv.com/forum/forum4/topic18214/messages/
%   Posted by: GMAC Administrator @ emotiv
%
% Assumes: 128 Hz Sampling Freq; matrix has EEG data only; no corrections!
%
% MDT
% 2016.01.18

fftlength  = 256;                    % make the window for sample length fftlength
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, 14); % match to number of channels

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

thetaIndex    = find(f>=4 & f<8);
alphaIndex    = find(f>=8 & f<12);
lowBetaIndex  = find(f>=12 & f<16);
highBetaIndex = find(f>=16 & f<25);
gammaIndex    = find(f>=25 & f<40);
totIndex      = find(f>=6 & f<=40);
outdata       = [];

med     = median(eeg.raw,2);             % remove median of each sample
eeg.raw = eeg.raw - repmat(med, 1, 14);

for j=2:size(eeg.raw,1)                  % limit slew rate
    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 section

a = 0.0078125;  % HPF filter coefs
b = 0.9921875;

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.filt(j,:) = eeg.raw(j,:) - preVal;
end                                          % end HPF

% Add more vectors to the "eeg" data structure!

eeg.theta    = [];
eeg.alpha    = [];
eeg.lowBeta  = [];
eeg.highBeta = [];
eeg.gamma    = [];
eeg.tot      = [];
eeg.totmed   = [];

% Main loop!

for k = fftlength:32:size(eeg.filt,1)                            % step through every quarter second starting at first possible sample
    spectrum     = fft(eeg.filt(k-fftlength+1:k,:) .* hanning);  % apply window to HP filtered data
    spectrum     = sqrt(spectrum .* conj(spectrum));             % get magnitude
    eeg.theta    = [eeg.theta; k sum(spectrum(thetaIndex,:))];   % append total spectral power in band, including sample index k
    eeg.alpha    = [eeg.alpha; k sum(spectrum(alphaIndex,:))];
    eeg.lowBeta  = [eeg.lowBeta; k sum(spectrum(lowBetaIndex,:))];
    eeg.highBeta = [eeg.highBeta; k sum(spectrum(highBetaIndex,:))];
    eeg.gamma    = [eeg.gamma; k sum(spectrum(gammaIndex,:))];
    eeg.tot      = [eeg.tot; k sum(spectrum(totIndex,:))];
end

% PLOT

timeRange = 30:2500;
channel   = 15;

figure;
plot(eeg.theta(timeRange,1), eeg.theta(timeRange,channel), 'b');
hold on
plot(eeg.alpha(timeRange,1), eeg.alpha(timeRange,channel), 'r');
title('AF3')   % Hopefully AF3 is 2 -- need to check
legend('theta','alpha');
xlabel('Sample Number (Start of FFT; 32 sample steps)')
ylabel('Amplitude (squared?)');
ylim([0 2500]);

a2tr = eeg.alpha(:,channel)./eeg.theta(:,channel);

figure;
plot(eeg.alpha(timeRange,1), a2tr(timeRange), 'k');
legend('alpha to theta ratio');
xlabel('Sample Number (Start of FFT; 32 sample steps)')
ylabel('Ratio (Alpha to Theta)');