-
Notifications
You must be signed in to change notification settings - Fork 0
/
spect_analysisV2_pipeline.m
168 lines (132 loc) · 5.86 KB
/
spect_analysisV2_pipeline.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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
% this code finds the epoc specific time independent spectrograms using
% pmtm()
%%
nwues=length(AllEEGData_complete_reshaped);
goodelectrodes=[1:126];
currwue=1:6;
freqcutoffout=100;freqcutoffin=1;
%Analyzing complete spectrogram(not temporal) for each of the sub epocs of every trial:
%which frequencies get modulated max during the sub epocs will be analyzed.
%based on these freqs then will the spectral purturbation be analyzed.
events_onoffsetkin_data(2).sessno(2).data(5,11)=events_onoffsetkin_data(2).sessno(2).data(5,11)-1; % for some reason only the 8th traial in 2 wue and 2nd sesion is one shorter in Mara_reshaped
for i=currwue
nsess=length(AllEEGData_complete_reshaped(i).sessno);
max_subepoc1=0;
for j=1:nsess
ntrials=length(AllEEGData_complete_reshaped(i).sessno(j).trial);
max_subepoc=diff(events_onoffsetkin_data(i).sessno(j).data);
max_subepoc=max(max(max_subepoc));
if max_subepoc>max_subepoc1
max_subepoc1=max_subepoc;
end
end
numchans=size(AllEEGData_complete(3).sessno(1).data,1);
specOutmat=NaN*ones(numchans,max_subepoc1,4,ntrials,nsess);
datOutmat=zeros(size(specOutmat));
for j=1:nsess
ntrials=length(AllEEGData_complete_reshaped(i).sessno(j).trial);
for k=1:ntrials
trialen=size(AllEEGData_complete_reshaped(i).sessno(j).trial(k).data,2);
for m=1:4; %looping over sub epocs
startime=events_onoffsetkin_data(i).sessno(j).data(m,k);
endtime=events_onoffsetkin_data(i).sessno(j).data(m+1,k);
for l=1:numchans
dat=AllEEGData_complete(i).sessno(j).data(l,startime:endtime);
specOut=pmtm(dat);
querypnts=linspace(1,length(dat),max_subepoc1); %interpolating to same length as max subepoc
interpdSpec=interp1(specOut,querypnts,'PCHIP');
specOutmat(l,:,m,k,j)=interpdSpec; %specOut: chan x freq x epoc x trial x session
end
end
end
end
spectdata_SUBEPOC(i).data=specOutmat;
end
[pval sigVal]=signifreq(spectdata_SUBEPOC);
%% NOrmalizing the spectdata subepoc for every trial:
spectdata_SUBEPOC_norm=spectdata_SUBEPOC;
for wueno=currwue
sz=size(spectdata_SUBEPOC(wueno).data);
loopsz=prod(sz(3:5));
freqlen=size(spectdata_SUBEPOC(wueno).data,2);
freqlen=[1:freqlen].*500/freqlen;
goodidx=find(freqlen<freqcutoffout);
for i=1:sz(3)
for j=1:sz(4)
for k=1:sz(5)
mn=nanmean(spectdata_SUBEPOC(wueno).data(:,goodidx,i,j,k),2);
stde=nanstd(spectdata_SUBEPOC(wueno).data(:,goodidx,i,j,k),0,2);
mn=repmat(mn,1,length(goodidx));
stde=repmat(stde,1,length(goodidx));
spectdata_SUBEPOC_norm(wueno).data(:,goodidx,i,j,k)=...
(spectdata_SUBEPOC(wueno).data(:,goodidx,i,j,k)-mn)./stde;
end
end
end
end
clearvars mn stde loopsz sz
[pval_norm sigVal_norm]=signifreq(spectdata_SUBEPOC_norm);
%% conversion to dB scale for every trial:
spectdata_SUBEPOC_dB=spectdata_SUBEPOC;
for wueno=currwue
sz=size(spectdata_SUBEPOC(wueno).data);
loopsz=prod(sz(3:5));
freqlen=size(spectdata_SUBEPOC(wueno).data,2);
freqlen=[1:freqlen].*500/freqlen;
goodidx=find(freqlen<freqcutoffout);
for i=1:sz(3)
for j=1:sz(4)
for k=1:sz(5)
mx=max(spectdata_SUBEPOC(wueno).data(:,goodidx,i,j,k),[],2);
mx=repmat(mx,1,length(goodidx));
spectdata_SUBEPOC_dB(wueno).data(:,goodidx,i,j,k)=...
10*log10((spectdata_SUBEPOC(wueno).data(:,goodidx,i,j,k))./mx);
end
end
end
end
clearvars mx loopsz sz
[pval_dB sigVal_dB]=signifreq(spectdata_SUBEPOC_dB);
%% plotting:
imgstoreloc1='C:\Users\PROVOST\Google Drive ([email protected])\scripts\sasicaout\sasica\wicaOn\pipeout\';
imgstoreloc_norm='C:\Users\PROVOST\Google Drive ([email protected])\scripts\sasicaout\sasica\wicaOn\pipeout\norm\';
imgstoreloc_dB='C:\Users\PROVOST\Google Drive ([email protected])\scripts\sasicaout\sasica\wicaOn\pipeout\dB\';
imgst={imgstoreloc1, imgstoreloc_norm,imgstoreloc_dB };
for i=1:3
imgstoreloc=imgst{1,i};
%---------
if i==1
spectdata_SUBEPOC1=spectdata_SUBEPOC;sigval=sigVal;tag='';
elseif i==2
spectdata_SUBEPOC1=spectdata_SUBEPOC_norm;sigval=sigVal_norm;tag='normalized,';
elseif i==3
spectdata_SUBEPOC1=spectdata_SUBEPOC_dB;sigval=sigVal_dB;tag='dB,';
end
for wueno=currwue
freqlen=size(spectdata_SUBEPOC1(wueno).data,2);
freqlen=[1:freqlen].*500/freqlen;
goodidx=find(freqlen<freqcutoffout);
if ~mkdir([imgstoreloc sprintf('wueno%d',wueno)]);
mkdir([imgstoreloc sprintf('wueno%d',wueno)]);
end
for chnonew=1:(size(spectdata_SUBEPOC1(wueno).data,1))
figure();
col=hsv(4);
%goodelectrodes=[11:36 53:58 76:77 122:123 125:126];
%totsesstrials=(size(spectdata_SUBEPOC(wueno).data,4))*(size(spectdata_SUBEPOC(wueno).data,5));
for t=1:4;
plot(freqlen(goodidx),nanmean(nanmean(spectdata_SUBEPOC1(wueno).data(chnonew,(goodidx),t,:,:),4),5),...
'Color',col(t,:)); hold on;
end;
ax=gca;yvl=ax.YTick;yvl=yvl(2);
sigdat=sigval(wueno).data(chnonew,:);
plot(freqlen(goodidx),yvl*sigdat(goodidx),...
'-*');
title(sprintf('wue no. %d, chno. %d (%s avg across trial/sess)',wueno,goodelectrodes(chnonew),tag));
legend({'pre','reach','grasp','reverse'});grid on
xlabel(['Hz']);ylabel(['PSD avg across sessions']);
print([imgstoreloc sprintf('wueno%d',wueno) filesep...
sprintf('chno_%d ',goodelectrodes(chnonew))],'-djpeg');close all;
end
end
end