Home > chronux_1_15 > pointtimes > CrossSpecMatpt.m

CrossSpecMatpt

PURPOSE ^

SYNOPSIS ^

function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)

DESCRIPTION ^


 Multi-taper cross-spectral matrix - another routine, this one allows for multiple trials and channels 
 but does not do confidence intervals. Also this routine always averages
 over trials - point process as times

 Usage:

 [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
 Input: 
 Note units have to be consistent. See chronux.m for more information.
       data (as a struct array with dimensions channels x trials) - note
       that times of measurement have to be consistent, we assume all
       times are specified relative to the start time of the trials which
       are taken to be zero.
       win  (duration of non-overlapping window)
       trialduration (since it is not possible to infer trial duration
               from spike times, this is an optional argument. If not specified
               the routine uses the minimum and maximum spike time (across all
               channels and trials) as the window of calculation.) -
               optional
       params: structure with fields tapers, pad, Fs, fpass
       - optional
           tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional. If not 
                                                 specified, use [NW K]=[3 5]
            pad            (padding factor for the FFT) - optional. Defaults to 0.  
                       e.g. For N = 500, if PAD = 0, we pad the FFT 
                       to 512 points; if PAD = 2, we pad the FFT
                       to 2048 points, etc.
           Fs   (sampling frequency) - optional. Default 1.
           fpass    (frequency band to be used in the calculation in the form
                                   [fmin fmax])- optional. 
                                   Default all frequencies between 0 and Fs/2
 Output:
       Sc (cross spectral matrix frequency x channels x channels)
       Cmat Coherence matrix frequency x channels x channels
       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
       Cvec leading Eigenvector (frequency x channels)
       Cent A different measure of total coherence: GM/AM of SV^2s
       f (frequencies)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
0002 %
0003 %
0004 % Multi-taper cross-spectral matrix - another routine, this one allows for multiple trials and channels
0005 % but does not do confidence intervals. Also this routine always averages
0006 % over trials - point process as times
0007 %
0008 % Usage:
0009 %
0010 % [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
0011 % Input:
0012 % Note units have to be consistent. See chronux.m for more information.
0013 %       data (as a struct array with dimensions channels x trials) - note
0014 %       that times of measurement have to be consistent, we assume all
0015 %       times are specified relative to the start time of the trials which
0016 %       are taken to be zero.
0017 %       win  (duration of non-overlapping window)
0018 %       trialduration (since it is not possible to infer trial duration
0019 %               from spike times, this is an optional argument. If not specified
0020 %               the routine uses the minimum and maximum spike time (across all
0021 %               channels and trials) as the window of calculation.) -
0022 %               optional
0023 %       params: structure with fields tapers, pad, Fs, fpass
0024 %       - optional
0025 %           tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional. If not
0026 %                                                 specified, use [NW K]=[3 5]
0027 %            pad            (padding factor for the FFT) - optional. Defaults to 0.
0028 %                       e.g. For N = 500, if PAD = 0, we pad the FFT
0029 %                       to 512 points; if PAD = 2, we pad the FFT
0030 %                       to 2048 points, etc.
0031 %           Fs   (sampling frequency) - optional. Default 1.
0032 %           fpass    (frequency band to be used in the calculation in the form
0033 %                                   [fmin fmax])- optional.
0034 %                                   Default all frequencies between 0 and Fs/2
0035 % Output:
0036 %       Sc (cross spectral matrix frequency x channels x channels)
0037 %       Cmat Coherence matrix frequency x channels x channels
0038 %       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
0039 %       Cvec leading Eigenvector (frequency x channels)
0040 %       Cent A different measure of total coherence: GM/AM of SV^2s
0041 %       f (frequencies)
0042 d=ndims(data);
0043 if size(d,1)==1; error('Need multiple channels; are you sure your format is channels x trials ?');end;
0044 [C,Ntr]=size(data);
0045 mintime=0;
0046 if nargin < 3; [mintime,maxtime]=minmaxsptimes(data);clear mintime;
0047 else maxtime=T; end;
0048 if nargin < 4; params=[]; end;
0049 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
0050 clear err trialave params
0051 Nwin=round(Fs*win); % number of samples in window
0052 nfft=2^(nextpow2(Nwin)+pad);
0053 [f,findx]=getfgrid(Fs,nfft,fpass); 
0054 tapers=dpsschk(tapers,Nwin,Fs); % check tapers
0055 twin=linspace(0,win,Nwin); % times of occurrence of "samples" within window - times of evaluation of tapers
0056 
0057 Sc=zeros(length(findx),C,C);
0058 tn=mintime:win:maxtime-win;
0059 Nwins=length(tn);
0060 if d==3, % If there are multiple trials
0061 for iwin=1:Nwins,
0062     t=[tn(iwin) tn(iwin)+T];
0063     for i=1:Ntr, 
0064         data1=data(:,i);
0065         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
0066         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
0067         for k=1:C,
0068             for l=1:C,
0069                 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 
0070             Sc(:,k,l)=Sc(:,k,l)+spec;
0071             end
0072         end
0073     end
0074 end
0075 Sc=Sc/(Nwins*Ntr);
0076 end
0077 
0078 if d==2, % only one trial
0079 for iwin=1:Nwins,
0080         data1=data(:,i);
0081         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
0082         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
0083         for k=1:C,
0084             for l=1:C,
0085             Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
0086             end
0087         end
0088 end
0089 Sc=Sc/Nwins;
0090 end
0091 
0092 Cmat=Sc;
0093 Sdiag=zeros(length(findx),C);
0094 for k=1:C,
0095     Sdiag(:,k)=squeeze(Sc(:,k,k));
0096 end
0097 
0098 for k=1:C,
0099     for l=1:C,
0100         Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
0101     end
0102 end
0103 
0104 Ctot=zeros(length(findx),1); Cent=Ctot;
0105 Cvec=zeros(length(findx),C);
0106 for i=1:length(findx),
0107     [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
0108     Ctot(i)=s(1).^2/sum(s.^2); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);             
0109     Cvec(i,:)=transpose(u(:,1));
0110 
0111 end
0112

Generated on Tue 15-Aug-2006 22:51:57 by m2html © 2003