testscript

 function testscript(pname,direction,movingwin,segave,params,fscorr)

 This script runs a sequence of analysis steps using the test
 data contained in data. The data consists of a single tetrode
 recording from macaque area LIP during a memory saccade experiment
 a la Pesaran et al (2002). The data is already separated into
 spikes and LFPs. LFPs are contained in variable dlfp, spikes from two
 neurons are in in a struct array dsp, and event information is in
 the following set of variables:

 trialtimes - start times of trials
 fixon - fixation light comes on
 fixacq - fixation acquired
 targon - target light on
 targoff - target light off
 fixoff - fixation off
 saccade - saccade

 Note that spikes and event times are in seconds and the sampling
 frequency for the LFP in this experiment was 1kHz.

 the script takes the following input argument - 
 pname - path name on your computer where the data file LIPdata is stored.
 direction - target direction to be analysed (0-7)

 The remaining parameters control various computations and are discussed
 in chronux.m - type Help chronux.m for more information.
 
 if nargin < 4;
     error('Need 6 input parameters - see help');
 end;
 if nargin < 5 | isempty(params);
    [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
 end;
 if nargin < 6 | isempty(fscorr);
     fscorr=1;
 end;

0001 % function testscript(pname,direction,movingwin,segave,params,fscorr)
0002 %
0003 % This script runs a sequence of analysis steps using the test
0004 % data contained in data. The data consists of a single tetrode
0005 % recording from macaque area LIP during a memory saccade experiment
0006 % a la Pesaran et al (2002). The data is already separated into
0007 % spikes and LFPs. LFPs are contained in variable dlfp, spikes from two
0008 % neurons are in in a struct array dsp, and event information is in
0009 % the following set of variables:
0010 %
0011 % trialtimes - start times of trials
0012 % fixon - fixation light comes on
0013 % fixacq - fixation acquired
0014 % targon - target light on
0015 % targoff - target light off
0016 % fixoff - fixation off
0017 % saccade - saccade
0018 %
0019 % Note that spikes and event times are in seconds and the sampling
0020 % frequency for the LFP in this experiment was 1kHz.
0021 %
0022 % the script takes the following input argument -
0023 % pname - path name on your computer where the data file LIPdata is stored.
0024 % direction - target direction to be analysed (0-7)
0025 %
0026 % The remaining parameters control various computations and are discussed
0027 % in chronux.m - type Help chronux.m for more information.
0028 %
0029 % if nargin < 4;
0030 %     error('Need 6 input parameters - see help');
0031 % end;
0032 % if nargin < 5 | isempty(params);
0033 %    [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
0034 % end;
0035 % if nargin < 6 | isempty(fscorr);
0036 %     fscorr=1;
0037 % end;
0038 pname='data';
0039 params.Fs=1000; % sampling frequency
0040 params.fpass=[10 100]; % band of frequencies to be kept
0041 params. tapers=[3 5]; % taper parameters
0042 params.pad=2; % pad factor for fft
0043 params.err=[2 0.05];
0044 params.trialave=1;
0045 movingwin=[0.5 0.05];
0046 segave=1;
0047 direction=5;
0048 
0049 wintrig=[5*movingwin(1) 5*movingwin(1)];
0050 winseg=2*movingwin(1);
0051 %
0052 % Load data
0053 %
0054 eval(['load LIPdata.mat']);
0055 % %
0056 % % Create rearranged data blocks for further analysis: We are going to
0057 % % extract segments of data centered on the target off times from the first channel of LFP data and the from one of the two spike trains
0058 % %
0059 % %
0060 % indx1=find(targets==5);indx2=find(targets==1); % trials to preferred and antipreferred direction
0061 % E1=targon(indx1); E2=targon(indx2); % target on times trials to preferred adn anti-preferred directions
0062 % dlfp1=createdatamatc(dlfp(:,1),E1,Fs,wintrig);dlfp2=createdatamatc(dlfp(:,1),E2,Fs,wintrig); % extract event triggered segments of the first LFP channel
0063 % dsp1=createdatamatpt(dsp(1),E1,wintrig); dsp2=createdatamatpt(dsp(1),E2,wintrig); % the same for one of the spike trains
0064 
0065 % compute spectrum of the first few seconds of LFP channels 1-2
0066 NT=round(params.Fs*10*movingwin(1));
0067 data=dlfp(1:NT,:); data1=data(:,1:2);
0068 [S,f,Serr]=mtspectrumc(data1,params);
0069 figure;
0070 plot(f,10*log10(S),f,10*log10(Serr(1,:)),f,10*log10(Serr(2,:))); xlabel('Frequency Hz'); ylabel('Spectrum');
0071 %%% pause
0072 
0073 % compute derivative of the spectrum for the same data
0074 phi=[0 pi/2];
0075 [dS,f]=mtdspectrumc(data1,phi,params);
0076 figure;
0077 plot(f,dS(1,:),f,dS(2,:)); xlabel('frequency Hz'); ylabel('Derivatives'); legend('Time','Frequency');
0078 %%% pause
0079 
0080 % compute coherency between  channels 1-2 and  3-4
0081 data1=data(:,1:2);data2=data(:,3:4);
0082 [C,phi,S12,S1,S2,f,confC,phierr,Cerr]=coherencyc(data1,data2,params); 
0083 figure;plot(f,C,f,Cerr(1,:),f,Cerr(2,:));xlabel('frequency'); ylabel('Coherency');
0084 %%% pause
0085 
0086 % coherency matrix of data1
0087 [C,phi,S12,f,confC,phierr,Cerr]=cohmatrixc(data1,params);
0088 
0089 % compute spectrogram for 1-2
0090 [S,t,f,Serr]=mtspecgramc(data1,movingwin,params);
0091 figure;imagesc(t,f,10*log10(S)'); axis xy; colorbar
0092 %%% pause
0093 
0094 % compute time-frequency derivative of the spectrogram for 1-2
0095 phi=[0 pi/2];
0096 [dS,t,f]=mtdspecgramc(data1,movingwin,phi,params);
0097 % pause
0098 % figure;subplot(211);imagesc(t,f,squeeze(dS(1,:,:))'); axis xy; colorbar;
0099 % subplot(212);imagesc(t,f,squeeze(dS(2,:,:))'); axis xy; colorbar;
0100 % %%% pause
0101 
0102 % compute coherogram between 1-2 and 3-4
0103 NT=round(movingwin(1)*Fs);
0104 [C,phi,S12,S1,S2,t,f,confC,phierr,Cerr]=cohgramc(data1,data2,movingwin,params);
0105 figure;imagesc(t,f,C'); axis xy; colorbar;
0106 %%% pause
0107 
0108 % compute segmented spectrum of 1
0109 NT=10*round(winseg*Fs);
0110 data1=dlfp(1:NT,1);
0111 [S,f,varS,C,Serr]=mtspectrumsegc(data1,winseg,params,segave);
0112 figure; subplot(211);plot(f,10*log(S));
0113 imagesc(f,f,C); axis xy;colorbar;
0114 %%% pause
0115 
0116 % compute segmented coherency between 1 and 2
0117 NT=10*round(winseg*Fs);
0118 data1=dlfp(1:NT,1); data2=dlfp(1:NT,2);
0119 [C,phi,S12,S1,S2,f,confC,phierr,Cerr]=coherencysegc(data1,data2,winseg,params);
0120 figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1)); subplot(313);plot(f,10*log10(S2))
0121 %%% pause
0122 
0123 % compute spectrum of channel 1 triggered to events E
0124 E1=targon(find(targets==direction)); 
0125 data1=dlfp(:,1); 
0126 [S,f,Serr]=mtspectrumtrigc(data1,E1,wintrig,params);
0127 figure;plot(f,10*log10(S)'); axis xy; colorbar;
0128 %%% pause
0129 
0130 % compute spectrogram of channel 1 triggered to events E
0131 E1=targon(find(targets==direction)); 
0132 data1=dlfp(:,1); 
0133 [S,t,f,Serr]=mtspecgramtrigc(data1,E1,wintrig,movingwin,params);
0134 figure; imagesc(t,f,10*log10(S)'); axis xy; colorbar;
0135 %%% pause
0136 
0137 %
0138 % Analysis - point process stored as times
0139 %
0140 
0141 % dsp contains 2 channels of spikes
0142 
0143 data=extractdatapt(dsp,[20 30]); % extract spikes occurring between 20 and 30 seconds and compute their spectrum
0144 [S,f,R,Serr]=mtspectrumpt(data,params);
0145 figure; plot(f,10*log10(S),f,10*log10(Serr(1,:)),f,10*log10(Serr(2,:)));line(get(gca,'xlim'),[10*log10(R) 10*log10(R)]);
0146 %%% pause
0147 
0148 %
0149 % Compute the derivative of the spectrum
0150 %
0151 phi=[0 pi/2];
0152 [dS,f]=mtdspectrumpt(data,phi,params);
0153 figure; plot(f,dS);
0154 %%% pause
0155 
0156 %
0157 % Compute the derivative of the time-frequency spectrum
0158 %
0159 data=extractdatapt(dsp,[20 30]);
0160 data1=data(1); data2=data(2);fscorr=[];t=[];
0161 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencypt(data1,data2,params);
0162 figure; plot(f,C);
0163 %%% pause
0164 
0165 %
0166 % Compute event triggered average spectrum for one of the directions
0167 %
0168 E1=targon(find(targets==direction));
0169 data=dsp(1); 
0170 [S,f,R,Serr]=mtspectrumtrigpt(data,E1,wintrig,params);
0171 figure;plot(f,10*log10(S),f,10*log10(Serr(1,:)),f,10*log10(Serr(2,:))); line(get(gca,'xlim'),[10*log10(R) 10*log10(R)]); 
0172 %%% pause
0173 
0174 %
0175 % Compute the matrix of coherencies
0176 %
0177 data=extractdatapt(dsp,[20 30]);
0178 [C,phi,S12,f,zerosp,confC,phierr,Cerr]=cohmatrixpt(data,params,fscorr);
0179 
0180 %
0181 % Event triggered spectrogram - first way way
0182 %
0183 data=createdatamatpt(dsp(1),E1,wintrig);
0184 [S,t,f,R,Serr]=mtspecgrampt(data,movingwin,params);
0185 figure;imagesc(t,f,10*log10(S)'); axis xy; colorbar;
0186 %%% pause
0187 
0188 %
0189 % Derivative of the time-frequency spectrum
0190 %
0191 data=createdatamatpt(dsp(1),E1,wintrig);
0192 phi=[0 pi/2];
0193 [dS,t,f]=mtdspecgrampt(data,movingwin,phi,params);
0194 figure; subplot(211); imagesc(t,f,squeeze(dS(1,:,:))'); axis xy; colorbar;
0195 subplot(212); imagesc(t,f,squeeze(dS(2,:,:))'); axis xy; colorbar;
0196 
0197 %
0198 % Coherogram between the two spike trains
0199 %
0200 data1=createdatamatpt(dsp(1),E1,wintrig);
0201 data2=createdatamatpt(dsp(2),E1,wintrig);
0202 [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr);
0203 figure;imagesc(t,f,C');axis xy; colorbar
0204 % %%% pause
0205 
0206 %
0207 % Event Triggered spectrogram another way
0208 %
0209 data=dsp(1);
0210 [S,t,f,R,Serr]=mtspecgramtrigpt(data,E1,wintrig,movingwin,params);
0211 imagesc(t,f,10*log10(S)'); axis xy; colorbar
0212 %
0213 % Segmented spectrum
0214 %
0215 data=extractdatapt(dsp,[20 30]);
0216 data=data(1);
0217 [S,f,R,varS]=mtspectrumsegpt(data,winseg,params);
0218 plot(f,10*log10(S)); line(get(gca,'xlim'),[10*log10(R) 10*log10(R)]);
0219 %
0220 % Segmented coherency
0221 %
0222 data=extractdatapt(dsp,[20 30]);
0223 data1=data(1);data2=data(2);
0224 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegpt(data1,data2,winseg,params);
0225 figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1));subplot(313); plot(f,10*log10(S2))
0226 %
0227 % Analysis - hybrid: one continous and one point process stored as times
0228 %
0229 offset=1;
0230 data1=dlfp(20000:30000,1); data2=extractdatapt(dsp,[20 30],offset);data2=data2(1).times;
0231 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencycpt(data1,data2,params);
0232 figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1));subplot(313); plot(f,10*log10(S2))
0233 
0234 
0235 data1=dlfp(20000:30000,1); data2=extractdatapt(dsp,[20 30],offset);data2=data2(1).times;
0236 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegcpt(data1,data2,winseg,params,segave,fscorr);
0237 figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1));subplot(313); plot(f,10*log10(S2))
0238 
0239 
0240 data1=dlfp(20000:30000,1); data2=extractdatapt(dsp,[20 30],offset);data2=data2(1).times;
0241 [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgramcpt(data1,data2,movingwin,params,fscorr);
0242 figure; subplot(311); imagesc(t,f,C');axis xy; colorbar; subplot(312);imagesc(t,f,10*log10(S1)');axis xy; colorbar; subplot(313); imagesc(t,f,10*log10(S2)');axis xy; colorbar
0243 
0244 
0245 %  Analysis: Binned spike counts
0246 
0247 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0248 
0249 data=dN;
0250 [S,f,R,Serr]=mtspectrumpb(data,params);
0251 plot(f,10*log10(S),f,10*log10(Serr(1,:)), f,10*log10(Serr(2,:))); %line(get(gca,'xlim'),[10*log10(R) 10*log10(R)]);
0252 
0253 data=dN;
0254 phi=[0 pi/2];
0255 [dS,f]=mtdspectrumpb(data,phi,params);
0256 figure; plot(f,dS); 
0257 
0258 data=dN;
0259 data1=data(:,1); data2=data(:,2);fscorr=[];t=[];
0260 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencypb(data1,data2,params);
0261 figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1)); subplot(313); plot(f,10*log10(S2));
0262 
0263 E=targon(find(targets==direction)); E=E(find(E>20 & E<450)); [dN,t]=binspikes(dsp,params.Fs,[20 500]);data=dN(:,1);
0264 [S,f,R,Serr]=mtspectrumtrigpb(data,E,wintrig,params);
0265 figure;plot(f,10*log10(S),f,10*log10(Serr(1,:)),f,10*log10(Serr(2,:)));
0266 
0267 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0268 data=dN;
0269 [C,phi,S12,f,zerosp,confC,phierr,Cerr]=cohmatrixpb(data,params,fscorr);
0270 
0271 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0272 data=dN;
0273 [S,t,f,R,Serr]=mtspecgrampb(data,movingwin,params);
0274 figure;imagesc(t,f,10*log10(S)'); axis xy; colorbar
0275 
0276 clear R Serr Cerr S C phierr dN dS data1 data2
0277 
0278 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0279 data=dN;
0280 phi=[0 pi/2];
0281 [dS,t,f]=mtdspecgrampb(data,movingwin,phi,params);
0282 
0283 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0284 data=dN;
0285 data1=data(:,1); data2=data(:,2);
0286 [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgrampb(data1,data2,movingwin,params,fscorr);
0287 
0288 [dN,t]=binspikes(dsp,params.Fs); % extract spikes
0289 dN=dN(:,1);
0290 E=E(1:6);
0291 data=dN;
0292 [S,t,f,R,Serr]=mtspecgramtrigpb(data,E,wintrig,movingwin,params);
0293 
0294 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0295 data=dN;
0296 data=data(:,1);
0297 [S,f,R,varS]=mtspectrumsegpb(data,winseg,params,segave,fscorr);
0298 figure; plot(f,10*log10(S)); line(get(gca,'xlim'),10*log10([R R])); 
0299 
0300 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0301 data=dN;
0302 data1=data(:,1);data2=data(:,2);
0303 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegpb(data1,data2,winseg,params);
0304 
0305 %
0306 % Analysis - hybrid: one continous and one point process stored as counts
0307 %
0308 data1=dlfp(20000:30000,:);data1=data1(:,1:2); [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0309 data2=dN; data2=data2(1:end,:);
0310 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencycpb(data1,data2,params);
0311 
0312 data1=data1(:,1); data2=data2(:,1);
0313 [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegcpb(data1,data2,winseg,params,segave,fscorr);
0314 
0315 
0316 data1=dlfp(20000:30000,:); data1=data1(:,1:2);
0317 [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
0318 data2=dN; data2=data2(1:end,:);
0319 [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgramcpb(data1,data2,movingwin,params,fscorr);
0320