apply supper covariance (bci competition 3 dataset2 p300 speller)

nadia-aghili commented 4 years ago

1. > clc;
2. > 
3. > clear;
4. > close all;
5. > 
6. > 
7. > 
8. > TargetChar=[];
9. > StimulusType=[];
10. > 
11. > load 'Subject_A_Train.mat' % load data file
12. > window=160; % window after stimulus (1s)
13. > Signal1=double(Signal);
14. > number_train=1;
15. > 
16. > 
17. > % convert to double precision
18. > Signal1=double(Signal);
19. > Signal2=  Signal1(:,:,[11 34 48 51 54 56 60 2 6 9 13]);
20. > % Signal2=  Signal1(:,:,:);
21. > 
22. > Flashing=double(Flashing);
23. > StimulusCode=double(StimulusCode);
24. > StimulusType=double(StimulusType);
25. > 
26. > 
27. > 
28. > %% digital filtering
29. > N = 7794  ; f1 = 0.1  ; f2 = 10  ; Fs = 240 ;
30. > K1 = round((N*f1)/Fs)  ;    K2 = round((N*f2)/Fs) ;
31. > cc=11;
32. > FFT_X = fft(Signal2, [], 2);
33. > Filterd_FFT_X = zeros(85 , 7794 , cc);
34. > Filterd_FFT_X(: ,[K1:K2 , N-K2:N-K1] ,: ) = FFT_X (: , [K1:K2 , N-K2:N-K1] , : );
35. > Filterd_X = ifft(Filterd_FFT_X , [] , 2 ,'symmetric');
36. > Signal1=Filterd_X;
37. > %%
38. > 
39. > g=50;
40. > % for each character epoch
41. > b=0;
42. > for epoch=1:g   
43. >     epoch
44. >     % get reponse samples at start of each Flash
45. >     for n=2:size(Signal1,2)
46. >         if Flashing(epoch,n)==0 & Flashing(epoch,n-1)==1
47. >             b=b+1;
48. >             rowcol=StimulusCode(epoch,n-1); %%inja mifahmeh kodam satro sotoon toshan shode faghat
49. >             responses(:,:,b)=Signal1(epoch,n-24:n+window-25,:);  %%signal bad az har tahriko joda
50. >             if isempty(StimulusType)==0
51. >                 label(b,:)=unique(StimulusCode(epoch,n-1).*StimulusType(epoch,n-1));
52. >                 if label(b,1)~= 0
53. >                     label(b,1)=1;
54. >                 end
55. >             end
56. >             
57. >         end
58. >     end
59. >     
60. > end
61. > 
62. > X_train2= responses;
63. > xx = permute(X_train2 , [3 , 1 , 2]);   %Rearrange dimensions of N-D array  %%%180*g*160*11
64. > % X_train4 = reshape( xx, 180*g , 160*cc) ;
65. >  
66. > X_trains= X_train4;
67. > 
68. > %%% ebteda class ha ra joda kardam ta supper covariance besazam %%%%
69. > c=0;
70. > b=0;
71. > for i = 1:size(X_trains,1)
72. >     if label(i,:)==1
73. >         b=b+1;
74. >         x1(b,:,:)=X_trains(i,:,:);
75. >     else
76. >         c=c+1;
77. >         x2(c,:,:)=X_trains(i,:,:);
78. >     end
79. >     
80. > end
81. > 
82. > %% apply super cov
83. > 
84. > x1_mean = mean(x1,1);
85. > x1_mean = squeeze(x1_mean);
86. > 
87. > x2_mean = mean(x2,1);
88. > x2_mean = squeeze(x2_mean);
89. > %class1 for each trial
90. > for i_sup = 1 : size(x1,1)
91. >     i_sup
92. >     X_1 =squeeze( x1(i_sup,:,:));
93. >    XRAW_mean1 = [x1_mean';X_1'];
94. >    XRAW_mean11(i_sup,:,:) = cov(XRAW_mean1');
95. > end
96. > %class2 for each trial
97. > for i_sup = 1 : size(x2,1)
98. >     i_sup
99. >     X_2 =squeeze( x2(i_sup,:,:));
100. >    XRAW_mean2 = [x2_mean';X_2'];
101. >    XRAW_mean22(i_sup,:,:) = cov(XRAW_mean2');
102. > end
103. > 
104. > x_new_sup = cat(1,XRAW_mean11,XRAW_mean22);   % new response
105. > 
106. > X_train_sup= x_new_sup;
107. > xx_sup = permute(X_train_sup , [3 , 1 , 2]);
108. > X_train4_sup = reshape( xx_sup, 180*g , 4*cc*cc) ;
109. > 
110. > for d=1:size(X_train4_sup,1)
111. > dec_X_train_sup(d,:) = downsample( X_train4_sup(d,:) , 10);
112. > end
113. > 
114. >  for d=1:size(xx,1)
115. > dec_X_train_sup(d,:,:) = downsample( X_train4_sup(d,:) , 10);
116. >  end
117. > X_trains_sup= dec_X_train_sup;
118. > 
119. > c=0;
120. > b=0;
121. > for i = 1:size(X_trains_sup,1)
122. >     if label(i,:)==1
123. >         b=b+1;
124. >         x1_sup(b,:)=X_trains_sup(i,:);
125. >     else
126. >         c=c+1;
127. >         x2_sup(c,:)=X_trains_sup(i,:);
128. >     end
129. >     
130. > end
131. > 
132. > 
133. > 
134. > c=size(x2_sup,1);
135. > x_tr11_sup = [x1_sup; x1_sup ; x2_sup ;x1_sup ; x1_sup; x1_sup];
136. > c1= size(x1_sup,1);
137. > label1= ones(c1,1);
138. > label2= zeros(c,1);
139. > labeltr_sup = [label1; label1; label2 ; label1 ; label1  ; label1 ];
140. > 
141. > %% LDA train
142. > ldamodel1 = fitcdiscr(x_tr11_sup,labeltr_sup);
143. > %% SVM train
144. > SVMModel= fitcsvm(x_tr11_sup,labeltr_sup,'Standardize',...
145. >     true,'KernelFunction','rbf','BoxConstraint',5, 'KernelScale',20);
146. > 
147. >  %%
148. >  
149. > %%%%%%%%%%% test %%%%%%%%
150. > 
151. > Signal1 = Filterd_X(51:85,:,:);
152. > 
153. > screen = char('A','B','C','D','E','F',...
154. >     'G','H','I','J','K','L',...
155. >     'M','N','O','P','Q','R',...
156. >     'S','T','U','V','W','X',...
157. >     'Y','Z','1','2','3','4',...
158. >     '5','6','7','8','9','_');
159. > 
160. > b=0;
161. > t=35; %% TEDAD CHAR TEST
162. > event=[];
163. > t1=0;
164. > for epoch=51:85
165. >  %%%%%% event
166. >  t1=t1+1;
167. >     h(1)=0;
168. >     l=0;
169. >     for i=2:size(StimulusCode(epoch,:),2)
170. >         h(i)=StimulusCode(1,i);
171. >         if h(i)~= h(i-1) && h(i)~=0             
172. >             l=l+1;
173. >             eventtest(t1,l)=StimulusCode(epoch,i);          
174. >         end
175. >     end
176. >      %%%%%%%%
177. >     for n=2:size(Signal1,2)
178. >         if Flashing(epoch,n)==0 & Flashing(epoch,n-1)==1
179. >             b=b+1;
180. >             rowcol=StimulusCode(epoch,n-1); 
181. >             responsestest(:,:,b)=Signal1(t1,n-24:n+window-25,:); 
182. >            
183. >         end
184. >     end
185. > end
186. > 
187. > 
188. > X_testsa= responsestest;
189. > xtt = permute(X_testsa , [3 , 1 , 2]);   %Rearrange dimensions of N-D array  %%%180*g*160*11
190. > 
191. > x3_mean = mean(xtt,1);
192. > x3_mean = squeeze(x3_mean);
193. > 
194. > for i_sup = 1 : size(X_testsa,3)
195. >     i_sup
196. >     X_3 =squeeze( xtt(i_sup,:,:));
197. >    XRAW_mean3 = [x3_mean';X_3'];
198. >    XRAW_mean33(i_sup,:,:) = cov(XRAW_mean3');
199. > end
200. > 
201. > X_testsa_sup= XRAW_mean33;
202. > % xtt_sup = permute(X_testsa_sup , [3 , 1 , 2]);   %Rearrange dimensions of N-D array  %%%180*g*160*11
203. > X_testsa_sup = reshape(X_testsa_sup, 12*15*t , 4*cc*cc) ;
204. > % 
205. > for d=1:size(X_testsa_sup,1)
206. > dec_X_test_sup(d,:) = downsample( X_testsa_sup(d,:) , 10);
207. > end
208. > 
209. > X_testsa_supper = dec_X_test_sup  ;
210. > 
211. > [aa1,classified_train1 ]= predict(ldamodel1,X_testsa_supper);
212. > Score_pos1 = classified_train1(:,2);
213. > Score_pos1 = normalize(Score_pos1);
214. > 
215. > Score_pos=Score_pos1;
216. > for i=1:t
217. >     for j=1:12
218. >         c=find(eventtest(i,:)==j);
219. >         mean_scor((i-1)*12+j) = mean(Score_pos(180*(i-1)+c));  %%%barasase j chideh
220. >     end
221. > end
222. > 
223. > %% target detection
224. > for k = 1:2*t
225. >     [~ , ind(k)] = max(mean_scor((k-1)*6+1:(k*6)));
226. > end
227. > 
228. > col = ind(1:2:end)';
229. > row = ind(2:2:end)';
230. > row_col = col+(row-1)*6;
231. > target_predict = screen(row_col);
232. > test_str = TargetChar(51:85);
233. > target_ACC = sum(test_str(1:t) == target_predict')/length(target_predict);
234. > disp(target_ACC)
235. > 
236. > %%
237. > %%
238. > [SVMlabel,SVMscoresvm] = predict(SVMModel,X_testsa_supper);
239. > score_target=SVMscoresvm(:,2);
240. > 
241. > % sadafscore=SVMscores;
242. > % sadaflabel=SVMlabel;
243. > 
244. > 
245. > Score_possvm = SVMscoresvm(:,2);
246. > 
247. > for i=1:t
248. >     for j=1:12
249. >         c=find(eventtest(i,:)==j);
250. >         mean_scorsvm((i-1)*12+j) = mean(Score_possvm(180*(i-1)+c));  %%%barasase j chideh
251. >     end
252. > end
253. > 
254. > %% target detection
255. > for k = 1:2*t
256. >     [~ , ind(k)] = max(mean_scorsvm((k-1)*6+1:(k*6)));
257. > end
258. > 
259. > col = ind(1:2:end)';
260. > row = ind(2:2:end)';
261. > row_col = col+(row-1)*6;
262. > target_predictsvm = screen(row_col);
263. > 
264. > target_ACCsvm = sum(test_str(1:t) == target_predictsvm')/length(target_predictsvm);
265. > disp(target_ACCsvm)

nadia-aghili commented 4 years ago

this is my code i don't have a good classification accuracy with this method, can you help me what's my problem?

lkorczowski commented 4 years ago

28. > %% digital filtering
29. > N = 7794  ; f1 = 0.1  ; f2 = 10  ; Fs = 240 ;
30. > K1 = round((N*f1)/Fs)  ;    K2 = round((N*f2)/Fs) ;
31. > cc=11;
32. > FFT_X = fft(Signal2, [], 2);
33. > Filterd_FFT_X = zeros(85 , 7794 , cc);
34. > Filterd_FFT_X(: ,[K1:K2 , N-K2:N-K1] ,: ) = FFT_X (: , [K1:K2 , N-K2:N-K1] , : );
35. > Filterd_X = ifft(Filterd_FFT_X , [] , 2 ,'symmetric');
36. > Signal1=Filterd_X;

here is the issue I see in the preprocessing steps:

avoid fft filtering as it can distort the temporal signal in an unexpected way (filter design is hard ^^). Use instead common butterworth filtering (back-forward is better, use filtfilt in matlab).
for filter design, as parameters, here is one mistake: f1 = 0.1 ; f2 = 10won't work. As mentioned in several papers, lower frequencies will make poorly conditionned covariance matrix. You want a minimal of f1 = 1 for P300 as the usual windows length is 1 second. Lower values of f1 means your signal is NOT zero mean. Zero mean is a necessary condition. f2 = 10 is a bit too low but it should be fine (I use 16 or 25 for better results).
remember that fft is sampled on half the sample rate, therefore the correct line is: K1 = round((N*f1)/Fs/2) ; K2 = round((N*f2)/Fs/2) ; % don't forget to divide by 2

lkorczowski commented 4 years ago

What do Nand cc and why are they hardcoded ?

Usually, you want to avoid any hardcoded values. If N is the total number of samples, just read it from your data.

lkorczowski commented 4 years ago

69. > c=0;
70. > b=0;
71. > for i = 1:size(X_trains,1)
72. >     if label(i,:)==1
73. >         b=b+1;
74. >         x1(b,:,:)=X_trains(i,:,:);
75. >     else
76. >         c=c+1;
77. >         x2(c,:,:)=X_trains(i,:,:);
78. >     end
79. >     
80. > end

instead, use something like
x1 = X_trains(label == 1,:,:) x2 = X_trains(label ~= 1,:,:)

lkorczowski commented 4 years ago

If you want an example of how to do the super covariance matrix, just check that simple snippet (you can skip generate training-test set if you do it manually)

https://github.com/lkorczowski/GIPSEEG/blob/b3306ed1d8c5ddbb964c53f3e0c44d0a835dbede/scripts/BCI/script_rob_easy.m#L41-L50

lkorczowski / GIPSEEG

apply supper covariance (bci competition 3 dataset2 p300 speller) #1