forked from jakewesterberg/bastos-lab-hackathon-tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathanalysis.m
155 lines (123 loc) · 4.55 KB
/
analysis.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
%Good one: Lucky-Search-05182018-001
cfg = [];
cfg.viewmode = 'vertical';
% cfg.channel = 1:16;
% cfg.channel = 17:32;
cfg.channel = 33:48;
ft_databrowser(cfg, data)
U1 = 1:16; %Area PFC
U2 = 17:32; %Area 7A
U3 = 33:48; %Area V4 Foveal
%interpolate bad channels
% badchan1 = 7;
% badchan2 = 9;
% badchan3 = [];
%
% lfp(:,U1(badchan1),:) = (lfp(:,U1(badchan1)+1,:) + lfp(:,U1(badchan1)-1,:))./2;
% lfp(:,U2(badchan2),:) = (lfp(:,U2(badchan2)+1,:) + lfp(:,U2(badchan2)-1,:))./2;
% lfp(:,U3(badchan3),:) = (lfp(:,U3(badchan3)+1,:) + lfp(:,U3(badchan3)-1,:))./2;
figure;
subplot(1,3,1); imagesc(squeeze(mean(lfp(:,U1,:),3))'); caxis([-0.1 .1])
subplot(1,3,2); imagesc(squeeze(mean(lfp(:,U2,:),3))'); caxis([-0.1 .1])
title('LFP ERP')
subplot(1,3,3); imagesc(squeeze(mean(lfp(:,U3,:),3))'); caxis([-0.1 .1])
data =[];
for t = 1:size(lfp, 3)
data.time{t,1} = -1.5:0.001:3.0;
data.trial{t,1} = zeros(size(electrodeInfo,1), 4501); %num channel x num time
data.trial{t,1} = squeeze(lfp(:,:,t))';
end
data.fsample = 1000;
data.label = {};
%make each channel unique
for c = 1:size(electrodeInfo,1)
data.label{c} = [electrodeInfo.area{c} int2str(electrodeInfo.channel(c))];
end
data.trialinfo = zeros(size(lfp, 3), 1);
data.trialinfo(:,1) = 1:size(lfp, 3);
cfg = [];
cfg.latency = [-1 1];
data = ft_selectdata(cfg, data)
cfg = [];
% cfg.trials = 1:100;
data = ft_selectdata(cfg, data)
cfg = [];
cfg.method = 'mtmfft';
% cfg.tapsmofrq = 5;
cfg.taper = 'hanning';
cfg.output = 'pow';
cfg.keeptrials = 'no';
cfg.trials = 1:100;
cfg.foi = [0:150];
mpow = ft_freqanalysis(cfg, data)
% mpow = ft_freqdescriptives([], pow);
relpow1 = mpow.powspctrm(U1,:) ./ repmat(max(mpow.powspctrm(U1,:)), [16 ,1]);
relpow2 = mpow.powspctrm(U2,:) ./ repmat(max(mpow.powspctrm(U2,:)), [16 ,1]);
relpow3 = mpow.powspctrm(U3,:) ./ repmat(max(mpow.powspctrm(U3,:)), [16,1]);
figure;
subplot(1,3,1); imagesc(relpow1)
subplot(1,3,2); imagesc(relpow2)
title('LFP relative power')
subplot(1,3,3); imagesc(relpow3)
b1 = nearest(mpow.freq, 8); b2 = nearest(mpow.freq, 30);
t1 = nearest(mpow.freq, 1); t2 = nearest(mpow.freq, 4);
g1 = nearest(mpow.freq, 50); g2 = nearest(mpow.freq, 150);
figure;
subplot(1,3,1);
plot(mean(relpow1(:,t1:t2),2), 1:16, 'k'); hold on;
plot(mean(relpow1(:,b1:b2),2), 1:16, 'r'); hold on;
plot(mean(relpow1(:,g1:g2),2), 1:16, 'g'); hold on;
set(gca, 'ydir', 'reverse')
title('PFC')
subplot(1,3,2);
plot(mean(relpow2(:,t1:t2),2), 1:16, 'k'); hold on;
plot(mean(relpow2(:,b1:b2),2), 1:16, 'r'); hold on;
plot(mean(relpow2(:,g1:g2),2), 1:16, 'g'); hold on;
set(gca, 'ydir', 'reverse')
title('7A')
subplot(1,3,3);
plot(mean(relpow3(:,t1:t2),2), 1:16, 'k'); hold on;
plot(mean(relpow3(:,b1:b2),2), 1:16, 'r'); hold on;
plot(mean(relpow3(:,g1:g2),2), 1:16, 'g'); hold on;
set(gca, 'ydir', 'reverse')
title('V4')
epoched_U1 = lfp(:,U1,:);
epoched_U2 = lfp(:,U2,:);
epoched_U3 = lfp(:,U3,:);
NSampl=4501; %number of time samples per trial
NCh=16;
NTrl=size(lfp,3);
sig=0.4; %cortical conductivity
dis=200.*10.^-3; %inter-site distance in micrometers
sep = 2; %separation between channels
CSD1=zeros(NSampl,NCh-sep*2,NTrl);
CSD2=zeros(NSampl,NCh-sep*2,NTrl);
CSD3=zeros(NSampl,NCh-sep*2,NTrl);
for ch = sep+1:NCh-sep
CSD1(:,ch-sep,:) = -sig*(epoched_U1(:,ch-sep,:)-(2.*epoched_U1(:,ch,:))+epoched_U1(:,ch+sep,:))/(((dis).*sep).^2);
CSD2(:,ch-sep,:) = -sig*(epoched_U2(:,ch-sep,:)-(2.*epoched_U2(:,ch,:))+epoched_U2(:,ch+sep,:))/(((dis).*sep).^2);
CSD3(:,ch-sep,:) = -sig*(epoched_U3(:,ch-sep,:)-(2.*epoched_U3(:,ch,:))+epoched_U3(:,ch+sep,:))/(((dis).*sep).^2);
end
mCSD1 = squeeze(mean(CSD1(:,:,:),3)); %1000:2500
mCSD2 = squeeze(mean(CSD2(:,:,:),3));
mCSD3 = squeeze(mean(CSD3(:,:,:),3));
figure;
subplot(1,3,1); imagesc(mCSD1'); caxis([-0.05 .05])
subplot(1,3,2); imagesc(mCSD2'); caxis([-0.025 .025])
title('CSD')
subplot(1,3,3); imagesc(mCSD3'); caxis([-0.15 .15])
figure;
subplot(1,3,1); imagesc(mCSD1'); xlim([1400 2500]); caxis([-0.05 .05])
subplot(1,3,2); imagesc(mCSD2'); xlim([1400 2500]); caxis([-0.025 .025])
title('CSD stimulus')
subplot(1,3,3); imagesc(mCSD3'); xlim([1400 2500]); caxis([-0.15 .15])
figure;
subplot(1,3,1); imagesc(mCSD1'); xlim([3600 4000]); caxis([-0.05 .05])
subplot(1,3,2); imagesc(mCSD2'); xlim([3600 4000]); caxis([-0.025 .025])
title('CSD test on')
subplot(1,3,3); imagesc(mCSD3'); xlim([3600 4000]); caxis([-0.15 .15])
figure;
subplot(1,3,1); imagesc(mCSD1'); xlim([3800 4200]); caxis([-0.05 .05])
subplot(1,3,2); imagesc(mCSD2'); xlim([3800 4200]); caxis([-0.025 .025])
title('CSD response locked')
subplot(1,3,3); imagesc(mCSD3'); xlim([3800 4200]); caxis([-0.15 .15])