Skip to content
GitLab
Commit 60313847 authored by Marc van Wanrooij's avatar Marc van Wanrooij
Browse files

no message

parent c3eab578
Hide whitespace changes
Inline Side-by-side
DoingBayesianDataAnalysis/FilConJags.m
View file @ 60313847
function FilConJags
function FilConJags
% FILCONJAGS
%
% Filtration-condensation experiment
......
NIRS/nirs_av_figure5_onesubject.m
View file @ 60313847
......@@ -17,7 +17,7 @@ end
Txt = Txt(2:end,:);
figure(2)
clf
subj = 17;
subj = 1;
sel = nmbr == subj;
fnames = {d(sel).name};
nfiles = numel(fnames);
......@@ -144,7 +144,7 @@ on = on(1:2:end);
stim = nirs.event.stim;
stim = stim(1:2:end);
ustim = unique(stim);
% disp('===================')
for chanIdx = chan
y = S(:,chanIdx);
......@@ -164,8 +164,13 @@ for chanIdx = chan
x = HDR';
y = S(:,chanIdx);
b = glmfit(x,y,'normal'); % A AV V
[b,~,stats] = glmfit(x,y,'normal'); % A AV V
b
stats.t>0
stats.t
stats.p
b(1) = 0;
figure(2)
yfit = glmval(b,x,'identity');
......@@ -183,7 +188,7 @@ x = HDR';
end
% keyboard
subplot(121)
ax = axis;
for ii = 1:numel(on)
......@@ -239,7 +244,7 @@ legend('Visual','Audiovisual','Auditory','Location','SE');
pa_datadir;
print('-depsc','-painters',[mfilename 'HbO2']);
pause
%% HHB
close all
% Load mat files
......@@ -258,7 +263,7 @@ end
Txt = Txt(2:end,:);
figure(2)
clf
sel = nmbr == 17;
sel = nmbr == subj;
fnames = {d(sel).name};
nfiles = numel(fnames);
......@@ -391,9 +396,15 @@ ustim = unique(stim);
for chanIdx = 1
x = HDR';
y = S(:,chanIdx);
b = glmfit(x,y,'normal'); % A AV V
[b,~,stats] = glmfit(x,y,'normal'); % A AV V
offset = b(1);
b(1) = 0;
b
stats.t<0
stats.t
stats.p
figure(2)
yfit = glmval(b,x,'identity');
yfit(isnan(y)) = NaN;
......
NIRS/nirs_av_figure6_grandaverage.m
View file @ 60313847
% function nirs_av_grandaverage
% determine grand average
% signals are obtained from
% clear all
clear all
close all
% Load mat files
cd('/Users/marcw/DATA/NIRS/OXY3_v14112014'); %#ok<*UNRCH> % contains all relevant data files
......@@ -26,159 +26,186 @@ DOXYA = [];
DOXYV = [];
DOXYAV = [];
for ii = 1:(nnmbr-5)
sel = nmbr == unmbr(ii);
fnames = {d(sel).name};
nfiles = numel(fnames);
disp(ii)
disp(fnames)
S = [];
So = [];
T = [];
M = [];
E = [];
for jj = 1:nfiles
load(fnames{jj})
S = [S;nirs.signal]; %#ok<*AGROW>
So = [So;nirs.deepchan]; %#ok<*AGROW>
% S = [S;pa_zscore(nirs.signal)]; %#ok<*AGROW>
% So = [So;pa_zscore(nirs.deepchan)]; %#ok<*AGROW>
%% add timings, continuous
if ~isempty(E)
e = T(end)*nirs.Fs;
t = T(end);
T = [T;nirs.time+t];
E = [E [nirs.event.sample]+e];
else
T = [T;nirs.time];
E = [E [nirs.event.sample]];
end
fname = fnames{jj};
fname = fname(5:end-4);
sel = strcmpi(fname,Txt(:,1));
%% Check for interleaved blocks
modal = Txt(sel,5);
if strncmp(modal,'Random',5);
M = [M {nirs.event.stim}];
elseif strncmp(modal,'Audiovisual',11);
a = {nirs.event.stim};
for kk = 1:size(a,2)
a {kk} = 'AV';
end
M = [M a];
elseif strncmp(modal,'Auditory',8);
a = {nirs.event.stim};
for kk = 1:size(a,2)
a {kk} = 'A';
end
M = [M a];
elseif strncmp(modal,'Visual',6);
a = {nirs.event.stim};
for kk = 1:size(a,2)
a {kk} = 'V';
end
M = [M a];
else
M = [M {nirs.event.stim}];
end
end
S = pa_zscore(S);
% S = pa_zscore(So);
disp('==============')
sel = nmbr == unmbr(ii);
fnames = {d(sel).name};
nfiles = numel(fnames);
disp(ii)
disp(fnames)
S = [];
So = [];
T = [];
M = [];
E = [];
for jj = 1:nfiles
load(fnames{jj})
sci = nirs.sci;
sel = sci>=0.9;
% s = nirs.signal(:,sel);
s = nirs.signal;
% tmp = NaN(size(s(:,~sel)));
%% Block average
clear N
N.event.sample = E';
N.event.stim = M;
N.Fs = nirs.Fs;
N.fsdown = nirs.fsdown;
% s(:,~sel) = tmp;
figure(1)
clf
nchan = size(S,2);
mod = {'A';'V';'AV'};
k = 0;
muA = NaN(nchan,405);
muV = NaN(nchan,405);
muAV = NaN(nchan,405);
for chanIdx = 1:nchan
k = k+1;
for modIdx = 1:3
mu = pa_nirs_blockavg(N,S(:,chanIdx),mod{modIdx});
mu = mu(:,1:405);
x = 1:length(mu);
x = x/10;
whos mu
switch modIdx
case 1
muA(k,:) = nanmean(mu);
case 2
muV(k,:) = nanmean(mu);
case 3
muAV(k,:) = nanmean(mu);
end
figure(1)
plot(x,mu)
hold on
end
S = [S;s]; %#ok<*AGROW>
%% add timings, continuous
if ~isempty(E)
e = T(end)*nirs.Fs;
t = T(end);
T = [T;nirs.time+t];
E = [E [nirs.event.sample]+e];
else
T = [T;nirs.time];
E = [E [nirs.event.sample]];
end
oxy = muA(2:2:end,:);
OXYA = [OXYA;oxy];
oxy = muV(2:2:end,:);
OXYV = [OXYV;oxy];
oxy = muAV(2:2:end,:);
OXYAV = [OXYAV;oxy];
fname = fnames{jj};
fname = fname(5:end-4);
sel = strcmpi(fname,Txt(:,1));
doxy = muA(1:2:end,:);
DOXYA = [DOXYA;doxy];
doxy = muV(1:2:end,:);
DOXYV = [DOXYV;doxy];
doxy = muAV(1:2:end,:);
DOXYAV = [DOXYAV;doxy];
drawnow
%% Check for interleaved blocks
modal = Txt(sel,5);
if strncmp(modal,'Random',5);
M = [M {nirs.event.stim}];
elseif strncmp(modal,'Audiovisual',11);
a = {nirs.event.stim};
for kk = 1:size(a,2)
a {kk} = 'AV';
end
M = [M a];
elseif strncmp(modal,'Auditory',8);
a = {nirs.event.stim};
for kk = 1:size(a,2)
a {kk} = 'A';
end
M = [M a];
elseif strncmp(modal,'Visual',6);
a = {nirs.event.stim};
for kk = 1:size(a,2)
a {kk} = 'V';
end
M = [M a];
else
M = [M {nirs.event.stim}];
end
end
S = pa_zscore(S);
%% Block average
clear N
N.event.sample = E';
N.event.stim = M;
N.Fs = nirs.Fs;
N.fsdown = nirs.fsdown;
figure(1)
clf
nchan = size(S,2);
mod = {'A';'V';'AV'};
k = 0;
muA = NaN(nchan,405);
muV = NaN(nchan,405);
muAV = NaN(nchan,405);
sdA = NaN(nchan,405);
sdV = NaN(nchan,405);
sdAV = NaN(nchan,405);
for chanIdx = 1:nchan
k = k+1;
for modIdx = 1:3
try
mu = pa_nirs_blockavg(N,S(:,chanIdx),mod{modIdx});
mu = mu(:,1:405);
catch
mu = [];
end
x = 1:length(mu);
x = x/10;
switch modIdx
case 1
muA(k,:) = nanmean(mu);
sdA(k,:) = nanstd(mu)/sqrt(size(mu,1));
case 2
muV(k,:) = nanmean(mu);
sdV(k,:) = nanstd(mu)/sqrt(size(mu,1));
case 3
muAV(k,:) = nanmean(mu);
sdAV(k,:) = nanstd(mu)/sqrt(size(mu,1));
end
end
figure(1)
clf
pa_errorpatch(x,muA(k,:),sdA(k,:),'b');
hold on
pa_errorpatch(x,muV(k,:),sdV(k,:),'r');
pa_errorpatch(x,muAV(k,:),sdAV(k,:),'g');
legend('A','V','AV')
axis square
box off
title(['File = ' num2str(ii) ', channel = ' num2str(k)]);
end
oxy = muA(2:2:end,:);
OXYA = [OXYA;oxy];
oxy = muV(2:2:end,:);
OXYV = [OXYV;oxy];
oxy = muAV(2:2:end,:);
OXYAV = [OXYAV;oxy];
doxy = muA(1:2:end,:);
DOXYA = [DOXYA;doxy];
doxy = muV(1:2:end,:);
DOXYV = [DOXYV;doxy];
doxy = muAV(1:2:end,:);
DOXYAV = [DOXYAV;doxy];
drawnow
end
%%
close all
figure(666)
subplot(221)
oxy = OXYV;
mu = nanmean(oxy);
n = size(oxy,1);
sd = std(oxy)./sqrt(n);
sel = all(isnan(oxy'));
n = sum(~sel);
sd = nanstd(oxy)./sqrt(n);
t = 1:length(mu);
t = t/10;
patch([10 10 30 30],[-2 2 2 -2],[.9 .9 .9]);
hold on
pa_errorpatch(t,mu,sd,[0 0 .9])
pa_errorpatch(t,mu,2*sd,[0 0 .9])
hold on
oxy = OXYA;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.8 0 0]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.8 0 0]);
NCHAN(1) = n;
oxy = OXYAV;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[0 .7 0]);
sel = all(isnan(oxy'));
n = sum(~sel);
% n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[0 .7 0]);
NCHAN(2) = n;
oxy = OXYA+OXYV;
mu = nanmean(oxy);
n = size(oxy,1)
sel = all(isnan(oxy'));
n = sum(~sel);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.7 .7 .7]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.7 .7 .7]);
NCHAN(3) = n;
axis square;
......@@ -190,40 +217,45 @@ xlabel('Time re stimulus onset (s)');
ylabel('\DeltaHbO_2');
title('Normal-hearing, left & right AC)');
%%
subplot(223)
oxy = DOXYV;
mu = nanmean(oxy);
n = size(oxy,1);
sel = all(isnan(oxy'));
n = sum(~sel);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
sd = nanstd(s)./sqrt(n);
t = 1:length(mu);
t = t/10;
patch([10 10 30 30],[-2 2 2 -2],[.9 .9 .9]);
hold on
pa_errorpatch(t,mu,sd,[0 0 .9])
pa_errorpatch(t,mu,2*sd,[0 0 .9])
hold on
oxy = DOXYA;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.8 0 0]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.8 0 0]);
NCHAN(4) = n;
oxy = DOXYAV;
mu = nanmean(oxy);
n = size(oxy,1);
sel = all(isnan(oxy'));
n = sum(~sel);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[0 .7 0]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[0 .7 0]);
NCHAN(5) = n;
oxy = DOXYA+DOXYV;
mu = nanmean(oxy);
n = size(oxy,1);
sel = all(isnan(oxy'));
n = sum(~sel);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.7 .7 .7]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.7 .7 .7]);
NCHAN(6) = n;
axis square;
......@@ -234,8 +266,7 @@ ylim([-1.4 0.2])
xlabel('Time re stimulus onset (s)');
ylabel('\DeltaHbR');
% pa_errorpatch(t,mu,sd,'r')
% pa_errorpatch(t,mu,2*sd,'r')
% plot(nanmean(OXYA))
% plot(nanmean(OXYAV))
......@@ -271,7 +302,6 @@ for ii = (nnmbr-4):nnmbr
disp(ii)
disp(fnames)
S = [];
So = [];
T = [];
M = [];
E = [];
......@@ -322,7 +352,6 @@ for ii = (nnmbr-4):nnmbr
end
end
S = pa_zscore(S);
% So = pa_zscore(So);
%% Block average
......@@ -338,10 +367,10 @@ for ii = (nnmbr-4):nnmbr
mod = {'A';'V';'AV'};
k = 0;
muA = NaN(nchan,406);
muV = NaN(nchan,406);
muAV = NaN(nchan,406);
muV = NaN(nchan,406);
muAV = NaN(nchan,406);
for chanIdx = 1:nchan
k = k+1;
for modIdx = 1:3
......@@ -350,18 +379,20 @@ for ii = (nnmbr-4):nnmbr
x = x/10;
switch modIdx
case 1
muA(k,:) = nanmean(mu);
muA(k,:) = nanmean(mu);
case 2
muV(k,:) = nanmean(mu);
muV(k,:) = nanmean(mu);
case 3
muAV(k,:) = nanmean(mu);
muAV(k,:) = nanmean(mu);
end
figure(1)
plot(x,mu)
hold on
end
figure(1)
clf
plot(x,muA(k,:))
hold on
plot(x,muV(k,:))
plot(x,muAV(k,:))
legend('A','V','AV')
end
oxy = muA(2:2:end,:);
OXYA = [OXYA;oxy];
......@@ -376,7 +407,7 @@ for ii = (nnmbr-4):nnmbr
DOXYV = [DOXYV;doxy];
doxy = muAV(1:2:end,:);
DOXYAV = [DOXYAV;doxy];
drawnow
end
end
......@@ -389,26 +420,26 @@ oxy = OXYA;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
sd = nanstd(s)./sqrt(n);
t = 1:length(mu);
t = t/10;
patch([10 10 30 30],[-2 2 2 -2],[.9 .9 .9]);
hold on
pa_errorpatch(t,mu,sd,[0 0 .9])
pa_errorpatch(t,mu,2*sd,[0 0 .9])
hold on
oxy = OXYV;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.8 0 0]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.8 0 0]);
oxy = OXYAV;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[0 .7 0]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[0 .7 0]);
......@@ -416,8 +447,8 @@ oxy = OXYA+OXYV;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.7 .7 .7]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.7 .7 .7]);
axis square;
......@@ -434,34 +465,34 @@ oxy = DOXYA;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
sd = nanstd(s)./sqrt(n);
t = 1:length(mu);
t = t/10;
patch([10 10 30 30],[-2 2 2 -2],[.9 .9 .9]);
hold on
pa_errorpatch(t,mu,sd,[0 0 .9])
pa_errorpatch(t,mu,2*sd,[0 0 .9])
hold on
oxy = DOXYV;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);
pa_errorpatch(t,mu,sd,[.8 0 0]);
sd = nanstd(s)./sqrt(n);
pa_errorpatch(t,mu,2*sd,[.8 0 0]);
oxy = DOXYAV;
mu = nanmean(oxy);
n = size(oxy,1);
s = bsxfun(@minus,oxy,nanmean(oxy,2));
sd = std(s)./sqrt(n);