123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320 |
- % 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;
- pname='data';
- params.Fs=1000; % sampling frequency
- params.fpass=[10 100]; % band of frequencies to be kept
- params. tapers=[3 5]; % taper parameters
- params.pad=2; % pad factor for fft
- params.err=[2 0.05];
- params.trialave=1;
- movingwin=[0.5 0.05];
- segave=1;
- direction=5;
- wintrig=[5*movingwin(1) 5*movingwin(1)];
- winseg=2*movingwin(1);
- %
- % Load data
- %
- eval(['load LIPdata.mat']);
- % %
- % % Create rearranged data blocks for further analysis: We are going to
- % % 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
- % %
- % %
- % indx1=find(targets==5);indx2=find(targets==1); % trials to preferred and antipreferred direction
- % E1=targon(indx1); E2=targon(indx2); % target on times trials to preferred adn anti-preferred directions
- % dlfp1=createdatamatc(dlfp(:,1),E1,Fs,wintrig);dlfp2=createdatamatc(dlfp(:,1),E2,Fs,wintrig); % extract event triggered segments of the first LFP channel
- % dsp1=createdatamatpt(dsp(1),E1,wintrig); dsp2=createdatamatpt(dsp(1),E2,wintrig); % the same for one of the spike trains
- % compute spectrum of the first few seconds of LFP channels 1-2
- NT=round(params.Fs*10*movingwin(1));
- data=dlfp(1:NT,:); data1=data(:,1:2);
- [S,f,Serr]=mtspectrumc(data1,params);
- figure;
- plot(f,10*log10(S),f,10*log10(Serr(1,:)),f,10*log10(Serr(2,:))); xlabel('Frequency Hz'); ylabel('Spectrum');
- %%% pause
- % compute derivative of the spectrum for the same data
- phi=[0 pi/2];
- [dS,f]=mtdspectrumc(data1,phi,params);
- figure;
- plot(f,dS(1,:),f,dS(2,:)); xlabel('frequency Hz'); ylabel('Derivatives'); legend('Time','Frequency');
- %%% pause
- % compute coherency between channels 1-2 and 3-4
- data1=data(:,1:2);data2=data(:,3:4);
- [C,phi,S12,S1,S2,f,confC,phierr,Cerr]=coherencyc(data1,data2,params);
- figure;plot(f,C,f,Cerr(1,:),f,Cerr(2,:));xlabel('frequency'); ylabel('Coherency');
- %%% pause
- % coherency matrix of data1
- [C,phi,S12,f,confC,phierr,Cerr]=cohmatrixc(data1,params);
- % compute spectrogram for 1-2
- [S,t,f,Serr]=mtspecgramc(data1,movingwin,params);
- figure;imagesc(t,f,10*log10(S)'); axis xy; colorbar
- %%% pause
- % compute time-frequency derivative of the spectrogram for 1-2
- phi=[0 pi/2];
- [dS,t,f]=mtdspecgramc(data1,movingwin,phi,params);
- % pause
- % figure;subplot(211);imagesc(t,f,squeeze(dS(1,:,:))'); axis xy; colorbar;
- % subplot(212);imagesc(t,f,squeeze(dS(2,:,:))'); axis xy; colorbar;
- % %%% pause
- % compute coherogram between 1-2 and 3-4
- NT=round(movingwin(1)*Fs);
- [C,phi,S12,S1,S2,t,f,confC,phierr,Cerr]=cohgramc(data1,data2,movingwin,params);
- figure;imagesc(t,f,C'); axis xy; colorbar;
- %%% pause
- % compute segmented spectrum of 1
- NT=10*round(winseg*Fs);
- data1=dlfp(1:NT,1);
- [S,f,varS,C,Serr]=mtspectrumsegc(data1,winseg,params,segave);
- figure; subplot(211);plot(f,10*log(S));
- imagesc(f,f,C); axis xy;colorbar;
- %%% pause
- % compute segmented coherency between 1 and 2
- NT=10*round(winseg*Fs);
- data1=dlfp(1:NT,1); data2=dlfp(1:NT,2);
- [C,phi,S12,S1,S2,f,confC,phierr,Cerr]=coherencysegc(data1,data2,winseg,params);
- figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1)); subplot(313);plot(f,10*log10(S2))
- %%% pause
- % compute spectrum of channel 1 triggered to events E
- E1=targon(find(targets==direction));
- data1=dlfp(:,1);
- [S,f,Serr]=mtspectrumtrigc(data1,E1,wintrig,params);
- figure;plot(f,10*log10(S)'); axis xy; colorbar;
- %%% pause
- % compute spectrogram of channel 1 triggered to events E
- E1=targon(find(targets==direction));
- data1=dlfp(:,1);
- [S,t,f,Serr]=mtspecgramtrigc(data1,E1,wintrig,movingwin,params);
- figure; imagesc(t,f,10*log10(S)'); axis xy; colorbar;
- %%% pause
- %
- % Analysis - point process stored as times
- %
- % dsp contains 2 channels of spikes
- data=extractdatapt(dsp,[20 30]); % extract spikes occurring between 20 and 30 seconds and compute their spectrum
- [S,f,R,Serr]=mtspectrumpt(data,params);
- 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)]);
- %%% pause
- %
- % Compute the derivative of the spectrum
- %
- phi=[0 pi/2];
- [dS,f]=mtdspectrumpt(data,phi,params);
- figure; plot(f,dS);
- %%% pause
- %
- % Compute the derivative of the time-frequency spectrum
- %
- data=extractdatapt(dsp,[20 30]);
- data1=data(1); data2=data(2);fscorr=[];t=[];
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencypt(data1,data2,params);
- figure; plot(f,C);
- %%% pause
- %
- % Compute event triggered average spectrum for one of the directions
- %
- E1=targon(find(targets==direction));
- data=dsp(1);
- [S,f,R,Serr]=mtspectrumtrigpt(data,E1,wintrig,params);
- 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)]);
- %%% pause
- %
- % Compute the matrix of coherencies
- %
- data=extractdatapt(dsp,[20 30]);
- [C,phi,S12,f,zerosp,confC,phierr,Cerr]=cohmatrixpt(data,params,fscorr);
- %
- % Event triggered spectrogram - first way way
- %
- data=createdatamatpt(dsp(1),E1,wintrig);
- [S,t,f,R,Serr]=mtspecgrampt(data,movingwin,params);
- figure;imagesc(t,f,10*log10(S)'); axis xy; colorbar;
- %%% pause
- %
- % Derivative of the time-frequency spectrum
- %
- data=createdatamatpt(dsp(1),E1,wintrig);
- phi=[0 pi/2];
- [dS,t,f]=mtdspecgrampt(data,movingwin,phi,params);
- figure; subplot(211); imagesc(t,f,squeeze(dS(1,:,:))'); axis xy; colorbar;
- subplot(212); imagesc(t,f,squeeze(dS(2,:,:))'); axis xy; colorbar;
- %
- % Coherogram between the two spike trains
- %
- data1=createdatamatpt(dsp(1),E1,wintrig);
- data2=createdatamatpt(dsp(2),E1,wintrig);
- [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr);
- figure;imagesc(t,f,C');axis xy; colorbar
- % %%% pause
- %
- % Event Triggered spectrogram another way
- %
- data=dsp(1);
- [S,t,f,R,Serr]=mtspecgramtrigpt(data,E1,wintrig,movingwin,params);
- imagesc(t,f,10*log10(S)'); axis xy; colorbar
- %
- % Segmented spectrum
- %
- data=extractdatapt(dsp,[20 30]);
- data=data(1);
- [S,f,R,varS]=mtspectrumsegpt(data,winseg,params);
- plot(f,10*log10(S)); line(get(gca,'xlim'),[10*log10(R) 10*log10(R)]);
- %
- % Segmented coherency
- %
- data=extractdatapt(dsp,[20 30]);
- data1=data(1);data2=data(2);
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegpt(data1,data2,winseg,params);
- figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1));subplot(313); plot(f,10*log10(S2))
- %
- % Analysis - hybrid: one continous and one point process stored as times
- %
- offset=1;
- data1=dlfp(20000:30000,1); data2=extractdatapt(dsp,[20 30],offset);data2=data2(1).times;
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencycpt(data1,data2,params);
- figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1));subplot(313); plot(f,10*log10(S2))
- data1=dlfp(20000:30000,1); data2=extractdatapt(dsp,[20 30],offset);data2=data2(1).times;
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegcpt(data1,data2,winseg,params,segave,fscorr);
- figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1));subplot(313); plot(f,10*log10(S2))
- data1=dlfp(20000:30000,1); data2=extractdatapt(dsp,[20 30],offset);data2=data2(1).times;
- [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgramcpt(data1,data2,movingwin,params,fscorr);
- 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
- % Analysis: Binned spike counts
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- [S,f,R,Serr]=mtspectrumpb(data,params);
- 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)]);
- data=dN;
- phi=[0 pi/2];
- [dS,f]=mtdspectrumpb(data,phi,params);
- figure; plot(f,dS);
- data=dN;
- data1=data(:,1); data2=data(:,2);fscorr=[];t=[];
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencypb(data1,data2,params);
- figure; subplot(311); plot(f,C);subplot(312); plot(f,10*log10(S1)); subplot(313); plot(f,10*log10(S2));
- E=targon(find(targets==direction)); E=E(find(E>20 & E<450)); [dN,t]=binspikes(dsp,params.Fs,[20 500]);data=dN(:,1);
- [S,f,R,Serr]=mtspectrumtrigpb(data,E,wintrig,params);
- figure;plot(f,10*log10(S),f,10*log10(Serr(1,:)),f,10*log10(Serr(2,:)));
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- [C,phi,S12,f,zerosp,confC,phierr,Cerr]=cohmatrixpb(data,params,fscorr);
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- [S,t,f,R,Serr]=mtspecgrampb(data,movingwin,params);
- figure;imagesc(t,f,10*log10(S)'); axis xy; colorbar
- clear R Serr Cerr S C phierr dN dS data1 data2
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- phi=[0 pi/2];
- [dS,t,f]=mtdspecgrampb(data,movingwin,phi,params);
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- data1=data(:,1); data2=data(:,2);
- [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgrampb(data1,data2,movingwin,params,fscorr);
- [dN,t]=binspikes(dsp,params.Fs); % extract spikes
- dN=dN(:,1);
- E=E(1:6);
- data=dN;
- [S,t,f,R,Serr]=mtspecgramtrigpb(data,E,wintrig,movingwin,params);
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- data=data(:,1);
- [S,f,R,varS]=mtspectrumsegpb(data,winseg,params,segave,fscorr);
- figure; plot(f,10*log10(S)); line(get(gca,'xlim'),10*log10([R R]));
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data=dN;
- data1=data(:,1);data2=data(:,2);
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegpb(data1,data2,winseg,params);
- %
- % Analysis - hybrid: one continous and one point process stored as counts
- %
- data1=dlfp(20000:30000,:);data1=data1(:,1:2); [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data2=dN; data2=data2(1:end,:);
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencycpb(data1,data2,params);
- data1=data1(:,1); data2=data2(:,1);
- [C,phi,S12,S1,S2,f,zerosp,confC,phierr,Cerr]=coherencysegcpb(data1,data2,winseg,params,segave,fscorr);
- data1=dlfp(20000:30000,:); data1=data1(:,1:2);
- [dN,t]=binspikes(dsp,params.Fs,[20 30]); % extract spikes occurring between 20 and 30 seconds
- data2=dN; data2=data2(1:end,:);
- [C,phi,S12,S1,S2,t,f,zerosp,confC,phierr,Cerr]=cohgramcpb(data1,data2,movingwin,params,fscorr);
|