Upload files to 'features'

features/features/ICL_feature_extration_full.m

@@ -0,0 +1,291 @@
+function [featuresmat, feature_label] = ICL_feature_extration_full(EEG, study_ID, dataset_ID)
+
+if ~exist('study_ID', 'var') || isempty(study_ID)
+    study_ID = nan;
+end
+if ~exist('dataset_ID', 'var') || isempty(dataset_ID)
+    dataset_ID = nan;
+end
+
+% =========================================================================
+% EEG stuff 
+% =========================================================================
+
+% % for testing only
+% lower2 = 2;
+% EEG = pop_subcomp(EEG, setdiff(1:size(EEG.icaweights,1), ...
+%     randperm(size(EEG.icaweights,1), min(lower2,size(EEG.icaweights,1)))),0);
+
+assert(~isempty(EEG.icawinv), ...
+    'No ica weights in the current EEG dataset! Compute ICA on your data first.')
+n_ic= size(EEG.icawinv,2); % ncps is number of components
+
+% =========================================================================
+% Convert to average reference
+% =========================================================================
+
+pop_reref(EEG,[],'exclude',setdiff(1:EEG.nbchan,EEG.icachansind));
+
+% calculate ica activations if missing
+if isempty(EEG.icaact)
+    EEG.icaact = eeg_getdatact(EEG,'component',1:n_ic);
+end
+
+% =========================================================================
+% Computing dipole
+% =========================================================================
+
+% 1 dipole
+if isfield(EEG, 'dipfit') && ~isempty(EEG.dipfit) ...
+        && isfield(EEG.dipfit, 'model') && ~isempty(EEG.dipfit.model) ...
+        && length(EEG.dipfit.model) == n_ic
+    % reuse existing dipoles
+    ind2dp = cellfun(@(c) size(c, 1), {EEG.dipfit.model.posxyz}) == 2;
+    
+    % 2 dipole
+    EEG.dipfit2 = getfield(pop_multifit( ...
+        EEG, find(~ind2dp), 'threshold', 100, 'dipoles', 2), 'dipfit');
+    
+    % 1 dipole
+    for it = find(ind2dp)
+        EEG.dipfit.model(it).posxyz(2, :) = [];
+        EEG.dipfit.model(it).momxyz(2, :) = [];
+    end
+    EEG = pop_multifit(EEG, find(ind2dp), 'threshold', 100, 'dipoles', 1);
+    
+else
+    % recompute all
+    try
+        EEG = pop_multifit(EEG,[],'threshold',100);
+    catch
+        EEG = pop_dipfit_settings(EEG);
+        EEG = pop_multifit(EEG,[],'threshold',100);
+    end
+    
+    % 2 dipole
+    EEG.dipfit2 = getfield(pop_multifit( ...
+        EEG,[],'threshold',100,'dipoles',2),'dipfit');
+end
+
+
+% =========================================================================
+% Creating structure for sasica 
+% =========================================================================
+
+optvals = struct('FASTER_blinkchans',     '',...
+                'chancorr_channames',     '',...
+                'chancorr_corthresh',     0.6,...
+                'EOGcorr_Heogchannames',  [] ,...
+                'EOGcorr_corthreshH',     'auto',  ...
+                'EOGcorr_Veogchannames',  '',...
+                'EOGcorr_corthreshV',     'auto',...
+                'resvar_thresh',          15,...
+                'SNR_snrcut',             1,...
+                'SNR_snrBL',              [-Inf,0],...
+                'SNR_snrPOI',             [0  Inf],...
+                'trialfoc_focaltrialout', 'auto',...
+                'focalcomp_focalICAout',  'auto',...
+                'autocorr_autocorrint',    20,...
+                'autocorr_dropautocorr',  'auto');            
+            
+
+% Enable/disable fields
+fenable = struct('MARA_enable',      0,...
+                 'FASTER_enable',    1,...
+                 'ADJUST_enable',    1,...
+                 'chancorr_enable',  1,...
+                 'EOGcorr_enable',   0,...
+                 'resvar_enable',    0,... % !!! changed!
+                 'SNR_enable',       1,...
+                 'trialfoc_enable',  0,...
+                 'focalcomp_enable', 1,...
+                 'autocorr_enable',  1);
+        
+fields = regexp(fieldnames(fenable),'[^_]*','match')';
+
+for i_f = 1:length(fields)
+    instr.(fields{i_f}{1}).(fields{i_f}{2}) = ...
+        fenable.([fields{i_f}{1} '_' fields{i_f}{2}]);
+end
+
+fields2 = regexp(fieldnames(optvals),'[^_]*','match')';
+
+for i_f = 1:numel(fields2)
+    if ~isempty(strfind(fields2{i_f}{2},'channame'))
+        if ischar(optvals.([fields2{i_f}{1} '_' fields2{i_f}{2}]))
+            instr.(fields2{i_f}{1}).(fields2{i_f}{2}) = ...
+                eval(['[chnb(''' optvals.([fields2{i_f}{1} '_' fields2{i_f}{2}]) ''')]']);
+        else
+            instr.(fields2{i_f}{1}).(fields2{i_f}{2}) = ...
+                eval(['[chnb([' optvals.([fields2{i_f}{1} '_' fields2{i_f}{2}])  '])]']);
+        end
+    else
+        try
+            instr.(fields2{i_f}{1}).(fields2{i_f}{2}) = ...
+                eval(['[' optvals.([fields2{i_f}{1} '_' fields2{i_f}{2}]) ']']);
+        catch
+            instr.(fields2{i_f}{1}).(fields2{i_f}{2}) = ...
+                optvals.([fields2{i_f}{1} '_' fields2{i_f}{2}]) ;
+        end
+
+    end
+end
+instr.opts.noplot = 1;
+
+
+% =========================================================================
+% Running sasica
+% =========================================================================
+
+[EEG] = myeeg_SASICA(EEG,instr);
+disp 'leaving sasica'
+
+
+% =========================================================================
+% Generating topos
+% =========================================================================
+
+topo = zeros(n_ic,740);
+plotrad = zeros(n_ic,1);
+
+disp 'saving topomap'
+for i = 1:n_ic
+    scalpmap_norm = EEG.icawinv(:,i)/sqrt(sum(EEG.icawinv(:,i).^2));
+    [~,Zi,plotrad(i)] = topoplotFast( scalpmap_norm, EEG.chanlocs, 'chaninfo', EEG.chaninfo, ...
+        'shading', 'interp', 'numcontour', 3,'electrodes','on','noplot','on');
+    topo(i,:) = Zi(~isnan(Zi));
+end
+disp 'finished saving topomap'
+EEG.reject.SASICA.topo = topo;
+EEG.reject.SASICA.plotrad = plotrad;
+
+
+% =========================================================================
+% Generating psds
+% =========================================================================
+
+[EEG.reject.SASICA.spec, EEG.reject.SASICA.specvar, ...
+    EEG.reject.SASICA.speckrt] = eeg_rpsd(EEG, 100);
+
+
+% =========================================================================
+% Generating autocorrs
+% =========================================================================
+
+if EEG.trials == 1
+    if EEG.pnts / EEG.srate > 5
+        EEG.reject.SASICA.autocorr = eeg_autocorr_welch(EEG);
+    else
+        EEG.reject.SASICA.autocorr = eeg_autocorr(EEG);
+    end
+else
+    EEG.reject.SASICA.autocorr = eeg_autocorr_fftw(EEG);
+end
+
+
+% =========================================================================
+% Generate and save feature vectors
+% =========================================================================
+
+[~, feature_labels] = genfvect;
+featuresmat = nan(n_ic,length(feature_labels));
+
+disp 'saving features'
+for i = 1:n_ic
+    [featuresmat(i,:), feature_label] = genfvect(EEG,study_ID,dataset_ID,i);
+end
+
+           
+% =========================================================================
+% =========================================================================
+% =========================================================================
+function [nb,channame,strnames] = chnb(channame, varargin)
+
+% chnb() - return channel number corresponding to channel names in an EEG
+%           structure
+%
+% Usage:
+%   >> [nb]                 = chnb(channameornb);
+%   >> [nb,names]           = chnb(channameornb,...);
+%   >> [nb,names,strnames]  = chnb(channameornb,...);
+%   >> [nb]                 = chnb(channameornb, labels);
+%
+% Input:
+%   channameornb  - If a string or cell array of strings, it is assumed to
+%                   be (part of) the name of channels to search. Either a
+%                   string with space separated channel names, or a cell
+%                   array of strings.
+%                   Note that regular expressions can be used to match
+%                   several channels. See regexp.
+%                   If only one channame pattern is given and the string
+%                   'inv' is attached to it, the channels NOT matching the
+%                   pattern are returned.
+%   labels        - Channel names as found in {EEG.chanlocs.labels}.
+%
+% Output:
+%   nb            - Channel numbers in labels, or in the EEG structure
+%                   found in the caller workspace (i.e. where the function
+%                   is called from) or in the base workspace, if no EEG
+%                   structure exists in the caller workspace.
+%   names         - Channel names, cell array of strings.
+%   strnames      - Channel names, one line character array.
+narginchk(1, 2);
+if nargin == 2
+    labels = varargin{1};
+else
+
+    try
+        EEG = evalin('caller','EEG');
+    catch
+        try
+            EEG = evalin('base','EEG');
+        catch
+            error('Could not find EEG structure');
+        end
+    end
+    if not(isfield(EEG,'chanlocs'))
+        error('No channel list found');
+    end
+    EEG = EEG(1);
+    labels = {EEG.chanlocs.labels};
+end
+if iscell(channame) || ischar(channame)
+
+    if ischar(channame) || iscellstr(channame)
+        if iscellstr(channame) && numel(channame) == 1 && isempty(channame{1})
+            channame = '';
+        end
+        tmp = regexp(channame,'(\S*) ?','tokens');
+        channame = {};
+        for i = 1:numel(tmp)
+            if iscellstr(tmp{i}{1})
+                channame{i} = tmp{i}{1}{1};
+            else
+                channame{i} = tmp{i}{1};
+            end
+        end
+        if isempty(channame)
+            nb = [];
+            return
+        end
+    end
+    if numel(channame) == 1 && not(isempty(strmatch('inv',channame{1})))
+        cmd = 'exactinv';
+        channame{1} = strrep(channame{1},'inv','');
+    else
+        channame{1} = channame{1};
+        cmd = 'exact';
+    end
+    nb = regexpcell(labels,channame,[cmd 'ignorecase']);
+
+elseif isnumeric(channame)
+    nb = channame;
+    if nb > numel(labels)
+        nb = [];
+    end
+end
+channame = labels(nb);
+strnames = sprintf('%s ',channame{:});
+if not(isempty(strnames))
+    strnames(end) = [];
+end

features/features/eeg_autocorr.m

@@ -0,0 +1,21 @@
+function resamp = eeg_autocorr(EEG, pct_data)
+
+if ~exist('pct_data', 'var') || isempty(pct_data)
+    pct_data = 100;
+end
+
+% calc autocorrelation
+X = fft(EEG.icaact, 2^nextpow2(2*EEG.pnts-1), 2);
+c = ifft(mean(abs(X).^2, 3), [], 2);
+if EEG.pnts < EEG.srate
+    ac = [c(:, 1:EEG.pnts, :) zeros(size(c, 1), EEG.srate - EEG.pnts + 1)];
+else
+    ac = c(:, 1:EEG.srate + 1, :);
+end
+
+% normalize by 0-tap autocorrelation
+ac = bsxfun(@rdivide, ac, ac(:, 1)); 
+
+% resample to 1 second at 100 samples/sec
+resamp = resample(ac', 100, EEG.srate)';
+resamp(:, 1) = [];

features/features/eeg_autocorr_fftw.m

@@ -0,0 +1,28 @@
+function resamp = eeg_autocorr_fftw(EEG, pct_data)
+
+if ~exist('pct_data', 'var') || isempty(pct_data)
+    pct_data = 100;
+end
+
+nfft = 2^nextpow2(2*EEG.pnts-1);
+
+% calc autocorrelation
+fftw('planner', 'hybrid');
+ac = zeros(size(EEG.icaact, 1), nfft, EEG.trials);
+for it = 1:size(EEG.icaact, 1)
+    X = fft(EEG.icaact(it, :, :), nfft, 2);
+    ac(it, :, :) = abs(X).^2;
+end
+ac = ifft(mean(ac, 3), [], 2);
+if EEG.pnts < EEG.srate
+    ac = [ac(:, 1:EEG.pnts, :) zeros(size(ac, 1), EEG.srate - EEG.pnts + 1)];
+else
+    ac = ac(:, 1:EEG.srate + 1, :);
+end
+
+% normalize by 0-tap autocorrelation
+ac = bsxfun(@rdivide, ac(:, 1:EEG.srate + 1, :), ac(:, 1)); 
+
+% resample to 1 second at 100 samples/sec
+resamp = resample(ac', 100, EEG.srate)';
+resamp(:, 1) = [];

features/features/eeg_autocorr_welch.m

@@ -0,0 +1,43 @@
+function ac = eeg_autocorr_welch(EEG, pct_data)
+
+% clean input cutoff freq
+if ~exist('pct_data', 'var') || isempty(pct_data)
+    pct_data = 100;
+end
+
+% setup constants
+ncomp = size(EEG.icaweights, 1);
+n_points = min(EEG.pnts, EEG.srate * 3);
+nfft = 2^nextpow2(n_points * 2 - 1);
+cutoff = floor(EEG.pnts / n_points) * n_points;
+index = bsxfun(@plus, ceil(0:n_points / 2:cutoff - n_points), (1:n_points)');
+
+% separate data segments
+if pct_data ~=100
+    rng(0)
+    n_seg = size(index, 2) * EEG.trials;
+    subset = randperm(n_seg, ceil(n_seg * pct_data / 100)); % need to find a better way to take subset
+    rng('shuffle') % remove duplicate data first (from 50% overlap)
+    temp = reshape(EEG.icaact(:, index, :), [ncomp size(index) .* [1 EEG.trials]]);
+    segments = temp(:, :, subset);
+else
+    segments = reshape(EEG.icaact(:, index, :), [ncomp size(index) .* [1 EEG.trials]]);
+end
+
+%% calc autocorrelation
+fftpow = abs(fft(segments, nfft, 2)).^2;
+ac = ifft(mean(fftpow, 3), [], 2);
+
+% normalize
+if EEG.pnts < EEG.srate
+    ac = [ac(:, 1:EEG.pnts, :) ./ (ac(:, 1) * (n_points:-1:1) / (n_points)) ...
+        zeros(ncomp, EEG.srate - n_points + 1)];
+else
+    ac = ac(:, 1:EEG.srate + 1, :) ./ ...
+        (ac(:, 1) * [(n_points:-1:n_points - EEG.srate + 1) ...
+        max(1, n_points - EEG.srate)] / (n_points));
+end
+
+% resample to 1 second at 100 samples/sec
+ac = resample(ac', 100, EEG.srate)';
+ac(:, 1) = [];

features/features/eeg_rpsd.m

@@ -0,0 +1,61 @@
+function [psdmed, psdvar, psdkrt] = eeg_rpsd(EEG, nfreqs, pct_data)
+
+% clean input cutoff freq
+nyquist = floor(EEG.srate / 2);
+if ~exist('nfreqs', 'var') || isempty(nfreqs)
+    nfreqs = nyquist;
+elseif nfreqs > nyquist
+    nfreqs = nyquist;
+end
+if ~exist('pct_data', 'var') || isempty(pct_data)
+    pct_data = 100;
+end
+
+% setup constants
+ncomp = size(EEG.icaweights, 1);
+n_points = min(EEG.pnts, EEG.srate);
+window = hamming(n_points)';
+cutoff = floor(EEG.pnts / n_points) * n_points;
+index = bsxfun(@plus, ceil(0:n_points / 2:cutoff - n_points), (1:n_points)');
+rng(0)
+n_seg = size(index, 2) * EEG.trials;
+subset = randperm(n_seg, ceil(n_seg * pct_data / 100)); % need to improve this
+rng('shuffle')
+
+try
+    %% slightly faster but high memory use
+    % calculate windowed spectrums
+    temp = reshape(EEG.icaact(:, index, :), [ncomp size(index) .* [1 EEG.trials]]);
+    temp = bsxfun(@times, temp(:, :, subset), window);
+    temp = fft(temp, n_points, 2);
+    temp = abs(temp).^2;
+    temp = temp(:, 2:nfreqs + 1, :) * 2 / (EEG.srate*sum(window.^2));
+    if nfreqs == nyquist
+        temp(:, end, :) = temp(:, end, :) / 2; end
+
+    % calculate outputs
+    psdmed = db(median(temp, 3));
+    psdvar = db(var(temp, [], 3), 'power');
+    psdkrt = db(kurtosis(temp, [], 3)/4);
+
+catch
+    %% lower memory but slightly slower
+    psdmed = zeros(ncomp, nfreqs);
+    psdvar = zeros(ncomp, nfreqs);
+    psdkrt = zeros(ncomp, nfreqs);
+    for it = 1:ncomp
+        % calculate windowed spectrums
+        temp = reshape(EEG.icaact(it, index, :), [1 size(index) .* [1 EEG.trials]]);
+        temp = bsxfun(@times, temp(:, :, subset), window);
+        temp = fft(temp, n_points, 2);
+        temp = temp .* conj(temp);
+        temp = temp(:, 2:nfreqs + 1, :) * 2 / (EEG.srate*sum(window.^2));
+        if nfreqs == nyquist
+            temp(:, end, :) = temp(:, end, :) / 2; end
+
+        psdmed(it, :) = db(median(temp, 3));
+        psdvar(it, :) = db(var(temp, [], 3), 'power');
+        psdkrt(it, :) = db(kurtosis(temp, [], 3)/4);
+    end
+end
+end

features/features/genfvect.m

@@ -0,0 +1,124 @@
+function [vec_feature, vec_label] = genfvect(EEG, study_ID, dataset_ID, ICnum)
+% check input
+if ~exist('EEG', 'var') || isempty(EEG);
+    study_ID = 1;
+    dataset_ID = 1;
+    ICnum = 1;
+    EEG = gen_eeg;
+end
+
+% Getting a couple more features
+dip_single = [EEG.dipfit.model(ICnum).posxyz EEG.dipfit.model(ICnum).momxyz];
+dip_double1 = [EEG.dipfit2.model(ICnum).posxyz(1, :) EEG.dipfit2.model(ICnum).momxyz(1, :)];
+dip_double2 = [EEG.dipfit2.model(ICnum).posxyz(2, :) EEG.dipfit2.model(ICnum).momxyz(2, :)];
+
+vec_feature = [];
+vec_label = {};
+index = 0;
+
+% book keeping
+update_fvec([study_ID, dataset_ID, ICnum], {'Study ID', 'Dataset ID', 'IC#'}, 3)
+
+% scalp topography
+update_fvec(EEG.reject.SASICA.topo(ICnum, :), 'topo image', 740)
+
+% SASICA
+update_fvec(EEG.reject.SASICA.icaautocorr(ICnum),'SASICA autocorrelation', 1)
+update_fvec(EEG.reject.SASICA.icafocalcomp(ICnum),'SASICA focal topo', 1)
+update_fvec(EEG.reject.SASICA.icaSNR(ICnum),'SASICA snr', 1)
+update_fvec(EEG.reject.SASICA.var(ICnum),'SASICA ic variance', 1)
+
+% ADJUST
+update_fvec(EEG.reject.SASICA.icaADJUST.diff_var(ICnum),'ADJUST diff_var', 1)
+update_fvec(EEG.reject.SASICA.icaADJUST.meanK(ICnum),'ADJUST Temporal Kurtosis', 1)
+update_fvec(EEG.reject.SASICA.icaADJUST.SED(ICnum),'ADJUST Spatial Eye Difference', 1)
+update_fvec(EEG.reject.SASICA.icaADJUST.SAD(ICnum),'ADJUST spatial average difference', 1)
+update_fvec(EEG.reject.SASICA.icaADJUST.GDSF(ICnum),'ADJUST General Discontinuity Spatial Feature', 1)
+update_fvec(EEG.reject.SASICA.icaADJUST.nuovaV(ICnum),'ADJUST maxvar/meanvar', 1) % should remove
+
+% FASTER
+% 1 Median gradient value, for high frequency stuff
+% 2 (NOT USED) Mean slope around the LPF band (spectral)
+% 3 Kurtosis of spatial map
+% 4 Hurst exponent
+% 5 (NOT USED) Eyeblink correlations
+label = {'FASTER Median gradient value','FASTER Kurtosis of spatial map','FASTER Hurst exponent'};
+update_fvec(EEG.reject.SASICA.icaFASTER.listprops(ICnum,[1 3 4]),label, 3)
+
+
+% basic info
+update_fvec(EEG.nbchan,'number of channels', 1)
+nic = size(EEG.icaact,1);
+update_fvec(nic,'number of ICs', 1)
+update_fvec(EEG.reject.SASICA.plotrad(ICnum),'topoplot plot radius', 1)
+update_fvec(EEG.trials>1,'epoched?', 1)
+update_fvec(EEG.srate,'sampling rate', 1)
+update_fvec(EEG.pnts,'number of data points', 1)
+
+% dipole info
+update_fvec(EEG.dipfit.model(ICnum).rv,'dip1 rv (SASICA)', 1)
+label = {'dip1 posx', 'dip1 posy', 'dip1 posz', 'dip1 momx', 'dip1 momy', 'dip1 momz'};
+update_fvec(dip_single,label, 6)
+update_fvec(EEG.dipfit2.model(ICnum).rv,'dip2 rv', 1)
+label = {'dip2_1 posx', 'dip2_1 posy', 'dip2_1 posz', 'dip2_1 momx', 'dip2_1 momy', 'dip2_1 momz'};
+update_fvec(dip_double1,label, 6)
+label = {'dip2_2 posx', 'dip2_2 posy', 'dip2_2 posz', 'dip2_2 momx', 'dip2_2 momy', 'dip2_2 momz'};
+update_fvec(dip_double2,label, 6)
+
+% spectrum info
+update_fvec(EEG.reject.SASICA.spec(ICnum, :), 'psd_med', 100)
+update_fvec(EEG.reject.SASICA.specvar(ICnum, :), 'psd_var', 100)
+update_fvec(EEG.reject.SASICA.speckrt(ICnum, :), 'psd_kurt', 100)
+
+% autocorrelation function
+update_fvec(EEG.reject.SASICA.autocorr(ICnum, :), 'autocorr', 100)
+
+
+function update_fvec(val, label, len)
+    vec_feature(index + (1:len)) = val;
+    if iscellstr(label)
+        vec_label(index + (1:len)) = label;
+    elseif ischar(label)
+        [vec_label{index + (1:len)}] = deal(label);
+    end
+    index = index + len;
+end
+
+end
+
+function EEG = gen_eeg
+EEG = eeg_emptyset;
+
+EEG.reject.SASICA.icaresvar = 1;
+EEG.reject.SASICA.icaautocorr = 1;
+EEG.reject.SASICA.icafocalcomp = 1;
+EEG.reject.SASICA.icaSNR = 1;
+EEG.reject.SASICA.var = 1;
+EEG.reject.SASICA.icaADJUST.diff_var = 1;
+EEG.reject.SASICA.icaADJUST.meanK = 1;
+EEG.reject.SASICA.icaADJUST.SED = 1;
+EEG.reject.SASICA.icaADJUST.SAD = 1;
+EEG.reject.SASICA.icaADJUST.GDSF = 1;
+EEG.reject.SASICA.icaADJUST.nuovaV = 1;
+
+EEG.dipfit.model.rv = 0;
+EEG.dipfit.model.posxyz = zeros(1,3);
+EEG.dipfit.model.momxyz = zeros(1,3);
+EEG.dipfit2.model.rv = 0;
+EEG.dipfit2.model.posxyz(1,:) = zeros(1,3);
+EEG.dipfit2.model.momxyz(1,:) = zeros(1,3);
+EEG.dipfit2.model.posxyz(2,:) = zeros(1,3);
+EEG.dipfit2.model.momxyz(2,:) = zeros(1,3);
+
+EEG.reject.SASICA.topo = zeros(1,740);
+EEG.reject.SASICA.spec = zeros(1,100);
+EEG.reject.SASICA.specvar = zeros(1,100);
+EEG.reject.SASICA.speckrt = zeros(1,100);
+EEG.reject.SASICA.autocorr = zeros(1,100);
+EEG.reject.SASICA.icaFASTER.listprops = zeros(1,5);
+
+EEG.nbchan = 1;
+EEG.icaact = ones(2);
+EEG.reject.SASICA.plotrad = 1;
+EEG.pnts = 1;
+end