-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsa_freqanalyzespindles.m
123 lines (100 loc) · 3.9 KB
/
sa_freqanalyzespindles.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% classifyspindles
% by Til Ole Bergmann 2013
% last modified 2016/11/22 by TOB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% analyzes time-frequency-representation of spindle events based on candidate epochs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [spindle_eventdata, spindle_freq] = sa_freqanalyzespindles(gcfg, spindle_eventdata)
display('Classifying spindles...');
% select raw channels only
procchannels = {'_bp','_rms'};
goodchan =[];
for i = 1:length(spindle_eventdata.label)
for j = 1:length(procchannels) % loop over channels postfixes to exclude
k = strfind(spindle_eventdata.label{i},procchannels{j});
goodchan(i,j) = isempty(k);
end
end
goodchan = logical(prod(goodchan,2));
rawchannels = spindle_eventdata.label(logical(goodchan));
if ismember('all', gcfg.nnmfchannel)
gcfg.nnmfchannel = rawchannels;
end
%% TFR
windowlength = (size(spindle_eventdata.trial{1},2)-1)/spindle_eventdata.fsample; % length of time window
% standard
cfg = [];
cfg.method = 'mtmconvol';
cfg.taper = 'hanning';
cfg.channel = rawchannels;
cfg.keeptrials = 'yes';
cfg.toi = [-windowlength/2:0.005:windowlength/2];
maxfoi = min(200,spindle_eventdata.fsample/5);
cfg.foi = [5:1:maxfoi];
for i = 1:maxfoi, cycles(i) = floor(100/(1000/i)); end % keep windows about 100ms with integer number of cycles
cycles(cycles < 5) = 5; % but at least 5 cycles!
cfg.t_ftimwin = cycles(cfg.foi(1):cfg.foi(end))./cfg.foi;
% cfg.t_ftimwin = 5./cfg.foi;
cfg.output = 'pow';
cfg.polyremoval = 1; % 0 = mean (default), 1 = linear
spindle_freq = ft_freqanalysis(cfg, spindle_eventdata);
% cfg = [];
% cfg.method = 'wavelet';
% cfg.taper = 'hanning';
% cfg.channel = rawchannels;
% cfg.keeptrials = 'yes';
% cfg.toi = [-windowlength/2:0.01:windowlength/2];
% cfg.foi = [4:1:200];
% cfg.width = 9;
% cfg.gwidth = 3;
% cfg.output = 'pow';
% cfg.polyremoval = 0; % 0 = mean (default), 1 = linear
% spindle_freq = ft_freqanalysis(cfg, spindle_eventdata);
% % low frequencies (with hanning tapers)
% cfg = [];
% cfg.method = 'mtmconvol';
% cfg.taper = 'hanning';
% cfg.channel = rawchannels;
% cfg.keeptrials = 'yes';
% cfg.foi = [4:1:29];
% cfg.toi = [-windowlength/2:0.01:windowlength/2];
% cfg.t_ftimwin = 5./cfg.foi;
% cfg.output = 'pow';
% spindle_freq = ft_freqanalysis(cfg, spindle_eventdata);
% % high frequencies (with multitapers)
% cfg = [];
% cfg.method = 'mtmconvol';
% cfg.taper = 'dpss';
% cfg.channel = rawchannels;
% cfg.keeptrials = 'yes';
% cfg.foi = [30:5:200];
% cfg.toi = [-windowlength/2:0.01:windowlength/2];
% cfg.t_ftimwin = ones(1,length(cfg.foi))*0.2;
% cfg.tapsmofrq = ones(1,length(cfg.foi))*10;
% cfg.output = 'pow';
% spindle_freq_high = ft_freqanalysis(cfg, spindle_eventdata);
%
% % concatenate low and high frequency TFRs
% spindle_freq.powspctrm = cat(3,spindle_freq.powspctrm,spindle_freq_high.powspctrm); % concatenate along dimension 3 (freq)
% spindle_freq.cumtapcnt = cat(2,spindle_freq.cumtapcnt,spindle_freq_high.cumtapcnt); % concatenate along dimension 2 (freq)
% spindle_freq.freq = [spindle_freq.freq spindle_freq_high.freq]; % concatenate vector (freq)
% spindle_freq.cfg2 = spindle_freq_high.cfg; % add second cfg field
%% NNMF
if gcfg.nnmf
% [A,Y,numIter,tElapsed,finalResidual]=nmfrule(X,k);
% nnmf.mat works with matlab2009 but not matlab2012!
cfg = [];
cfg.freqsize = size(spindle_freq.powspctrm); % [trials, channels, freqbins, timebins]
cfg.channel = [1];
cfg.freqwindow = [21:101];
cfg.timewindow = [150:250];
cfg.k = 2; % number of components
input = reshape(spindle_freq.powspctrm(:,cfg.channel,cfg.freqwindow,cfg.timewindow),cfg.freqsize(1),[]); % reshape
[w,h] = nnmf(input,cfg.k,'algorithm','als');
% output = rehape(input,cfg.freqsize(1),length(cfg.freqwindow),length(cfg.timewindow)); % re-reshape
end
spindles_freq.cfg = gcfg;
end % of function