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 one of the following
                    forms: 
                    (1) A numeric vector [TW K] where TW is the
                        time-bandwidth product and K is the number of
                        tapers to be used (less than or equal to
                        2TW-1). 
                    (2) A numeric vector [W T p] where W is the
                        bandwidth, T is the duration of the data and p 
                        is an integer such that 2TW-p tapers are used. In
                        this form there is no default i.e. to specify
                        the bandwidth, you have to specify T and p as
                        well. Note that the units of W and T have to be
                        consistent: if W is in Hz, T must be in seconds
                        and vice versa. Note that these units must also
                        be consistent with the units of params.Fs: W can
                        be in Hz if and only if params.Fs is in Hz.
                        The default is to use form 1 with TW=3 and K=5

            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...). 
                    -1 corresponds to no padding, 0 corresponds to padding
                    to the next highest power of 2 etc.
                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
                       to 512 points, if pad=1, we pad to 1024 points etc.
                       Defaults to 0.
           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 one of the following
0026 %                    forms:
0027 %                    (1) A numeric vector [TW K] where TW is the
0028 %                        time-bandwidth product and K is the number of
0029 %                        tapers to be used (less than or equal to
0030 %                        2TW-1).
0031 %                    (2) A numeric vector [W T p] where W is the
0032 %                        bandwidth, T is the duration of the data and p
0033 %                        is an integer such that 2TW-p tapers are used. In
0034 %                        this form there is no default i.e. to specify
0035 %                        the bandwidth, you have to specify T and p as
0036 %                        well. Note that the units of W and T have to be
0037 %                        consistent: if W is in Hz, T must be in seconds
0038 %                        and vice versa. Note that these units must also
0039 %                        be consistent with the units of params.Fs: W can
0040 %                        be in Hz if and only if params.Fs is in Hz.
0041 %                        The default is to use form 1 with TW=3 and K=5
0042 %
0043 %            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
0044 %                    -1 corresponds to no padding, 0 corresponds to padding
0045 %                    to the next highest power of 2 etc.
0046 %                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
0047 %                       to 512 points, if pad=1, we pad to 1024 points etc.
0048 %                       Defaults to 0.
0049 %           Fs   (sampling frequency) - optional. Default 1.
0050 %           fpass    (frequency band to be used in the calculation in the form
0051 %                                   [fmin fmax])- optional.
0052 %                                   Default all frequencies between 0 and Fs/2
0053 % Output:
0054 %       Sc (cross spectral matrix frequency x channels x channels)
0055 %       Cmat Coherence matrix frequency x channels x channels
0056 %       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
0057 %       Cvec leading Eigenvector (frequency x channels)
0058 %       Cent A different measure of total coherence: GM/AM of SV^2s
0059 %       f (frequencies)
0060 d=ndims(data);
0061 if size(d,1)==1; error('Need multiple channels; are you sure your format is channels x trials ?');end;
0062 [C,Ntr]=size(data);
0063 mintime=0;
0064 if nargin < 3; [mintime,maxtime]=minmaxsptimes(data);clear mintime;
0065 else maxtime=T; end;
0066 if nargin < 4; params=[]; end;
0067 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
0068 clear err trialave params
0069 Nwin=round(Fs*win); % number of samples in window
0070 nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
0071 [f,findx]=getfgrid(Fs,nfft,fpass); 
0072 tapers=dpsschk(tapers,Nwin,Fs); % check tapers
0073 twin=linspace(0,win,Nwin); % times of occurrence of "samples" within window - times of evaluation of tapers
0074 
0075 Sc=zeros(length(findx),C,C);
0076 tn=mintime:win:maxtime-win;
0077 Nwins=length(tn);
0078 if d==3, % If there are multiple trials
0079 for iwin=1:Nwins,
0080     t=[tn(iwin) tn(iwin)+T];
0081     for i=1:Ntr, 
0082         data1=data(:,i);
0083         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
0084         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
0085         for k=1:C,
0086             for l=1:C,
0087                 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 
0088             Sc(:,k,l)=Sc(:,k,l)+spec;
0089             end
0090         end
0091     end
0092 end
0093 Sc=Sc/(Nwins*Ntr);
0094 end
0095 
0096 if d==2, % only one trial
0097 for iwin=1:Nwins,
0098         data1=data(:,i);
0099         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
0100         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
0101         for k=1:C,
0102             for l=1:C,
0103             Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
0104             end
0105         end
0106 end
0107 Sc=Sc/Nwins;
0108 end
0109 
0110 Cmat=Sc;
0111 Sdiag=zeros(length(findx),C);
0112 for k=1:C,
0113     Sdiag(:,k)=squeeze(Sc(:,k,k));
0114 end
0115 
0116 for k=1:C,
0117     for l=1:C,
0118         Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
0119     end
0120 end
0121 
0122 Ctot=zeros(length(findx),1); Cent=Ctot;
0123 Cvec=zeros(length(findx),C);
0124 for i=1:length(findx),
0125     [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
0126     Ctot(i)=s(1)/sum(s); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);             
0127     Cvec(i,:)=transpose(u(:,1));
0128 
0129 end
0130

Generated on Fri 12-Aug-2011 11:36:15 by m2html © 2005