0001 function classify_ic_eeg(p)
0002
0003
0004
0005
0006
0007 disp(mfilename);
0008
0009
0010 method_list = {'kurtosis', 'entropy'};
0011 Nmethod = length(method_list);
0012 ica_program = 'apply_ica_eeg';
0013
0014
0015 data_file = fullfile(p.proj_root, p.eeg_dirname, p.trial_dirname, 'br_s1.eeg.mat');
0016 eeg = vb_load_meg_data(data_file);
0017 [Nch, Nt, Ntr] = size(eeg);
0018
0019
0020 ica_dir = fullfile(p.proj_root, p.eeg_dirname, p.ica_dirname);
0021 ica_result = fullfile(ica_dir, [ica_program '.mat']);
0022 load(ica_result)
0023
0024
0025 ic = zeros(Nch,Nt,Ntr);
0026 for tr = 1:Ntr
0027 ic(:, :, tr) = unmix*eeg(:, :, tr);
0028 end
0029
0030
0031 brain_ic_each_method = cell(Nmethod, 1);
0032 for method = 1:Nmethod
0033
0034 value = zeros(Nch,Ntr);
0035 for comp = 1:Nch
0036 for tr = 1:Ntr
0037 eval(['value(comp,tr)=' method_list{method} '(ic(comp,:,tr));']);
0038 end
0039 end
0040 mean_val = mean(value(:));
0041 std_val = std(value(:));
0042
0043 fprintf('Noise ICs [%s]: \n', method_list{method});
0044 brain_ic_each_method{method} = [];
0045 for comp = 1:Nch
0046
0047 Noutlier = 0;
0048 for tr = 1:Ntr
0049 value(comp, tr) = (value(comp,tr)-mean_val)/std_val;
0050 if abs(value(comp, tr)) > p.threshold_ourlier
0051 Noutlier = Noutlier+1;
0052 end
0053 end
0054
0055
0056 if Noutlier/Ntr <= p.threshold_num_of_outlier
0057 brain_ic_each_method{method} = [brain_ic_each_method{method} comp];
0058 else
0059 fprintf('[%s]', num2str(comp));
0060 end
0061 end
0062 fprintf('\n');
0063 end
0064
0065
0066 brain_ic = brain_ic_each_method{1};
0067 for method = 2:Nmethod
0068 brain_ic = intersect(brain_ic,brain_ic_each_method{method});
0069 end
0070 noise_ic = 1:Nch;
0071 noise_ic(brain_ic) = [];
0072
0073
0074 output_file = fullfile(ica_dir, [mfilename '.mat']);
0075 load(ica_result, 'w', 's', 'v', 'unmix', 'mix')
0076 save(output_file,'w', 's', 'v', 'unmix', 'mix', 'brain_ic', 'noise_ic')
0077
0078
0079
0080 function h = entropy(x, n)
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091 if nargin < 2
0092 n = 256;
0093 end
0094
0095 x = double(x);
0096 xh = hist(x(:), n);
0097 xh = xh / sum(xh(:));
0098
0099
0100 i = find(xh);
0101
0102 h = -sum(xh(i) .* log2(xh(i)));
0103
0104