-
Notifications
You must be signed in to change notification settings - Fork 1
/
online.m
125 lines (99 loc) · 3.96 KB
/
online.m
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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
function [Trials,SmoothedTotal]=online(Session,Mu,Sigma,classifier,FeaturesIndeces,alpha)
%{
Inputs: whole signal 4th run (s{1,4}=Session), mu, sigma, classifier, FeaturesIndeces
%}
global SubjectID
f=4:2:40;
%%Divide the run into the (30) trials
TrialStart=Session.event.position(Session.event.name==400);%400=MI initiation
TrialEnd=Session.event.position(Session.event.name==700);%700=Relax phase
Trials=[];
for i=1:(length(TrialStart))
Trials{i}=Session.data(TrialStart(i):TrialEnd(i),:)';
end
%%Divide each trial into 1-second windows shifted by 0.0625s
%(1 second=512 points; 0.0625 s = 32 points)
Aux=Trials;
Trials=[];
for i=1:length(Aux)
j = 1;
for t = 512:32:(size(Aux{i},2)-mod(size(Aux{i},2),32)) %we remove the last window (not of 1 second)
for idxChannels = 1:16
Trials{i}.Windows(idxChannels,:,j)=Aux{i}(idxChannels,t-512+1:t);
end
j= j+1;
end
end
%% Real pseudoonline
Evidence=[];
for i=1:length(Trials)
Evidence{i}=zeros(1,size(Trials{i}.Windows,3)+1);
Evidence{1,i}(1)=0.5; %for each trial we start by 0.5
Trials{i}.Features=zeros(size(Trials{i}.Windows,3),304);
for j=1:size(Trials{i}.Windows,3)
%%CAR filtering
mean_channels=[];
mean_channels = mean(Trials{i}.Windows(:,:,j));
Trials{i}.Windows(:,:,j) = Trials{i}.Windows(:,:,j)- mean_channels;
%%pwelch (already reshaped)
for idxChannels=1:16
%[psd,freqgrid] = pwelch(Trials{i}.Windows(idxChannels,:,j),0.5*512,0.4375*512,[],512);
%psd=log(psd);
%[freqs, idfreqs] = intersect(freqgrid,f);
%psd = psd(idfreqs);
%Trials{i}.Features(j,(19*(idxChannels-1)+1):19*idxChannels)=psd;
Trials{i}.Features(j,(19*(idxChannels-1)+1):19*idxChannels)=pwelch(Trials{i}.Windows(idxChannels,:,j),0.5*512,0.25*512,f,512);%0.4375*512
end
for idxChannels=1:16
Trials{i}.Features(j,19*(idxChannels-1)+1:19*idxChannels) = log(Trials{i}.Features(j,19*(idxChannels-1)+1:19*idxChannels));
end
%%normalization
Trials{i}.Features(j,:)=(Trials{i}.Features(j,:)-Mu)./Sigma;
%%predict (posterior prob)
[Trials{i}.Predicted(j),Trials{i}.PosteriorProb(j,:),~]=predict(classifier,Trials{i}.Features(j,FeaturesIndeces));
%%smoothing
Evidence{1,i}(j+1)=Evidence{1,i}(j)*alpha+(1-alpha)*Trials{i}.PosteriorProb(j,2);
fprintf('Trial %d Window %d \n',i,j);
end
end
%% Concatenate all the smoothed probabilities
SmoothedTotal=[];
for i=1:length(Trials)
SmoothedTotal=cat(2,SmoothedTotal,Evidence{1,i});
end
figure
sz = 4;
c = [0,0,1];
scatter((1-0.5*length(Evidence{1,1})*30/length(SmoothedTotal)):(30/length(SmoothedTotal)):(31-30/length(SmoothedTotal)-0.5*length(Evidence{1,1})*30/length(SmoothedTotal)),SmoothedTotal,sz,c,'filled')
hold on
cont=1-0.5*length(Evidence{1,1})*30/length(SmoothedTotal);
for i=1:length(Trials)
if cont~=1-0.5*length(Evidence{1,1})*30/length(SmoothedTotal)
vline(cont,'r','');
% else
% vline(cont,'k','');
end
cont=cont+length(Evidence{1,i})*30/length(SmoothedTotal);
end
xlabel('Trials')
ylabel('Smoothed probability')
title(sprintf('Accumulated Evidence - %s', SubjectID))
xlim([1-0.5*length(Evidence{1,1})*30/length(SmoothedTotal) 31-30/length(SmoothedTotal)-0.5*length(Evidence{1,1})*30/length(SmoothedTotal)])
%set(gca,'TickLength',[0 0]);
set(gca,'xtick',5:5:30);
%% Examples of plots
if SubjectID == "ak6"
figure
plot([0:(length(Evidence{1,10})-1)]/16,Evidence{1,10},'b')%otherwise 6 or 7
xlabel('Time [s]')
ylabel('Accumulated Evidence')
set(gca,'YLim',[0 1])
title('Good Trial')
figure
plot([0:(length(Evidence{1,22})-1)]/16,Evidence{1,22},'b')
xlabel('Time [s]')
ylabel('Accumulated Evidence')
set(gca,'YLim',[0 1])
title('Bad Trial')
end
end