% Att göra:
% 1) Att läsa in och skapa ett powerspektrum för vart och ett av de
% fyra fallen, samt olika kopplingar, precis som i figur 4a
% 2) Att läsa in och skapa DTF för vart och ett av de fyra fallen,
% pröva några olika frekvenser
% BADA: cue efter 2000 ms, 500 ms lång, 7000 ms långt delay
% CUE PFC: cue efter 2000 ms, 7000 ms lång, inget delay
% CUE PPC: cue efter 2000 ms, 7000 ms lång, inget delay
% INGEN CUE: total tid 7000 ms.
%
% Föreslår att vi skippar de första 500 ms på delayperioden,
% då får vi 4000 ms tid för varje typ
clear
filer{1}=extract_sim_dirs('MYSIMS\DTF\SIM_SERIER\SERIES_testDTF.txt');
filer{2}=extract_sim_dirs('MYSIMS\DTF\SIM_SERIER\SERIES_testDTF_CUE.txt');
titlename{1} = 'Delay';
titlename{2} = 'Cue till PFC';
% plots BUMP EEG.
% dirname: The name of the simulation directory
% smooth = degree of smoothing of histogram. Default = 2
%
% Version 2.0
% Author: Fredrik Edin, 2004
% Address: freedin@nada.kth.se
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% LOAD FILES %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
smooth = 2;
thisdir = pwd;
close all
LFP_ppc = cell(length(filer),1);
LFP_pfc = cell(length(filer),1);
for i = 1:length(filer)
for j = 1:length(filer{i})
cd( filer{i}{j} )
% load data files
load('Params.txt')
version = Params(1);
load('Q.txt')
nmod = Params(2);
if version == 3
tStart = Params(4);
tmp = [ Params(1:6) ; 100 ; 100 ];
for k = 1:nmod
tmp = [ tmp ; Params(3+4*k) ; 1000 ; Params(5+4*k) ];
tmp = [ tmp ; Params(4+4*k) ; 1000 ; Params(6+4*k) ];
end
Params = tmp;
elseif version == 4
tStart = Params(4);
Params = [ Params(1:6) ; 100 ; 100 ; Params(7:end) ];
elseif version == 5
tStart = Params(4);
end
tStop = Params(5);
% Name of the simulation
filename = strcat( filer{i}{j}, '/', 'LFP_prox_dist.txt' );
if exist(filename(1:end-4),'file')
movefile(filename(1:end-4), filename);
end
tmp=load( filename );
% plot the LFPgrams of the cell
filename = strcat( filer{i}{j}, '/', 'Q.txt' );
load( filename )
if length(Q) > 0 %& Q(1,5) ~= 0
% Normalize data to 0 mean and variance 1
tmp2 = tmp(Q(1,1)-tStart+1000:Q(1,1)-tStart+5000,2);
tmp2 = (tmp2 - mean(tmp2)) / std(tmp2);
tmp3 = tmp(Q(1,1)-tStart+1000:Q(1,1)-tStart+5000,3);
tmp3 = (tmp3 - mean(tmp3)) / std(tmp3);
% Store as LFPs
LFP_ppc{i} = [LFP_ppc{i} tmp2]; % OBS, har ar ordningen ratt igen eftersom jag anvander nya simulatorn
LFP_pfc{i} = [LFP_pfc{i} tmp3]; % OBS, har ar ordningen ratt igen eftersom jag anvander nya simulatorn
end
cd( thisdir )
end
LFP_ppc{i} = LFP_ppc{i} - ones(size(LFP_ppc{i},1),1)*mean(LFP_ppc{i});
LFP_pfc{i} = LFP_pfc{i} - ones(size(LFP_ppc{i},1),1)*mean(LFP_pfc{i});
end
dt = tmp(2,1)-tmp(1,1);
fS = 1000/dt;
fN = fS/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% BANDPASS %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fl = 0; % Lower cut-off in Hz
fu = 55; % Higher cut-off in Hz
Wp = [fl fu]/fN;
fwinL = 0.005; % The width of the difference between lower pass- and stopbands
fwinU = 0.05; % The width of the difference between upper pass- and stopbands
Ws = [];
Rp = 1;
Rs = 60;
% Create a bandpass filter
if fu == 100
Wp(2) = [];
Ws = min(Wp-fwinL,0.5*Wp);
[n,Wn] = ellipord(Ws, Wp, Rp, Rs);
[b,a] = ellip(n,1,60,Wn,'high');
else
Ws = min(Wp(end)+fwinU,Wp(end)+0.25*(1-Wp(end)));
if fl > 0
Ws = [max(Wp(1)-fwinL,0.5*Wp(1)) Ws];
else
Wp(1) = [];
end
[n,Wn] = ellipord(Wp, Ws, Rp, Rs);
[b,a] = ellip(n,1,60,Wn);
end
leg = [];
wa = [];
Pp = [];
Pf = [];
for i = 1:length(LFP_pfc)
Pp{i} = []; % Power spectrum of the
Pf{i} = [];
for j = 1:size(LFP_pfc{i},2)
filt_Lp{i}(:,j) = filter(b,a,reshape(LFP_ppc{i}(:,j),[],1));
filt_Lf{i}(:,j) = filter(b,a,reshape(LFP_pfc{i}(:,j),[],1));
[Pp{i}(:,j), Ppc{i}(:,2*j-1:2*j), f] = psd(filt_Lp{i}(:,j),1000,fS);
[Pf{i}(:,j), Pfc{i}(:,2*j-1:2*j), f] = psd(filt_Lf{i}(:,j),1000,fS);
Pp_o{i}(:,j) = psd(reshape(LFP_ppc{i}(:,j),[],1),1000,fS);
Pf_o{i}(:,j) = psd(reshape(LFP_pfc{i}(:,j),[],1),1000,fS);
end
end
cd(thisdir)
dt = tmp(2,1)-tmp(1,1);
% t = 0:dt:dt*(size(filt_Lp{1},1)-1);
f = 0:1:500;
color = 'rgb';
for i = 1:length(LFP_pfc)
meanPf_o(:,i) = [mean(Pf_o{i}')]';
meanPp_o(:,i) = [mean(Pp_o{i}')]';
end
figure(1)
clf
subplot(1,2,1)
color = 'bgrkm';
for i = 1:length(Pf_o)
plot(f,log10(Pf_o{i})*10,'color',color(i))
hold on
end
hold on
xlim([0 60])
ylim([-20 20])
legend({'0','100'},2)
title('PFC')
subplot(1,2,2)
for i = 1:length(Pp_o)
plot(f,log10(Pp_o{i})*10,'color',color(i))
hold on
end
hold on
xlim([0 60])
ylim([-20 20])
legend({'0','100'},3)
title('PPC')
cd(thisdir)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% DTF %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thisdir = pwd;
pmin = 1;
pmax = 100;
orsel = 'fpe';
for i = 1:length(filer)
A0 = [];
A25 = [];
for j = 1:1
% Do ARfit for prox_dist
% A0 = [A0 ; filt_Lp{i}(:,j) filt_Lf{i}(:,j)];
% A25 = [A25 ; filt_Lp{i}(:,j+1) filt_Lf{i}(:,j+1)];
% A50 = [A50 ; filt_Lp{i}(:,j+2) filt_Lf{i}(:,j+2)];
A0 = [A0 ; LFP_ppc{i}(:,j) LFP_pfc{i}(:,j)];
A25 = [A25 ; LFP_ppc{i}(:,j+1) LFP_pfc{i}(:,j+1)];
end
[w0{i},AR0{i},C0{i},SBC0,FPE0,th0]=arfit(A0,pmin,pmax,orsel,'zero');
order0(i) = size(AR0{i},2)/2;
[w25{i},AR25{i},C25{i},SBC25,FPE25,th25]=arfit(A25,pmin,pmax,orsel,'zero');
order25(i) = size(AR25{i},2)/2;
% [w50{i},AR50{i},C50{i},SBC50,FPE50,th50]=arfit(A50,pmin,pmax,orsel,'zero');
% order50(i) = size(AR50{i},2)/2;
% results of ARFIT
disp(['Order ' titlename{i} ' 0%:' int2str(order0(i))])
disp(['Order ' titlename{i} ' 100%:' int2str(order25(i))])
% disp(['Order ' titlename{i} ' 50%:' int2str(order50(i))])
end
% Do directed transfer function
f = 0:1:50;
dt = 0.001;
for i = 1:length(filer)
for j = 1:length(f)
[DTF0{i}{j},H0]=DirTransFunc(AR0{i},f(j),dt);
wp0{i}(j) = DTF0{i}{j}(2,1) - DTF0{i}{j}(1,2);
[DTF25{i}{j},H25]=DirTransFunc(AR25{i},f(j),dt);
wp25{i}(j) = DTF25{i}{j}(2,1) - DTF25{i}{j}(1,2);
% [DTF50{i}{j},H50]=DirTransFunc(AR50{i},f(j),dt);
% wp50{i}(j) = DTF50{i}{j}(2,1) - DTF50{i}{j}(1,2);
end
end
style = '--';
tit = [];
for i = 1:length(filer)
figure(2)
plot(f,wp0{i},style(1:i),f,wp25{i},style(1:i))%,f,wp50{i})
hold on
tit = [ tit titlename{i} ' (' style(1:i) ') ' ];
title(tit)
legend('0','100')%,'50')
ylim([-1 1])
line(xlim, [0 0], 'Color', 'k', 'Linestyle', '--')
yl = ylim;
xl = xlim;
text(xl(1) + 0.9*diff(xl), yl(1) + 0.25*diff(yl), 'P->F < F->P','HorizontalAlignment', 'right')
text(xl(1) + 0.9*diff(xl), yl(1) + 0.20*diff(yl), 'DTF > 0 means netflow goes in direction P->F','HorizontalAlignment', 'right')
disp(titlename{i})
disp(['0: ' num2str(sum(wp0{i})) ' 100: ' num2str(sum(wp25{i}))])
end
for i = 1:length(filer)
disp([titlename{i} ' 40Hz'])
disp(['0: ' num2str(wp0{i}(41)) ' 100: ' num2str(wp25{i}(41))])
end
for i = 1:length(filer)
disp([titlename{i} ' Gamma'])
disp(['0: ' num2str(sum(wp0{i}(31:end))) ' 100: ' num2str(sum(wp25{i}(31:end)))])
end
