Home > chronux_1_50 > 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 (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 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 (can take values -1,0,1,2...).
0028 %                    -1 corresponds to no padding, 0 corresponds to padding
0029 %                    to the next highest power of 2 etc.
0030 %                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
0031 %                       to 512 points, if pad=1, we pad to 1024 points etc.
0032 %                       Defaults to 0.
0033 %           Fs   (sampling frequency) - optional. Default 1.
0034 %           fpass    (frequency band to be used in the calculation in the form
0035 %                                   [fmin fmax])- optional.
0036 %                                   Default all frequencies between 0 and Fs/2
0037 % Output:
0038 %       Sc (cross spectral matrix frequency x channels x channels)
0039 %       Cmat Coherence matrix frequency x channels x channels
0040 %       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
0041 %       Cvec leading Eigenvector (frequency x channels)
0042 %       Cent A different measure of total coherence: GM/AM of SV^2s
0043 %       f (frequencies)
0044 d=ndims(data);
0045 if size(d,1)==1; error('Need multiple channels; are you sure your format is channels x trials ?');end;
0046 [C,Ntr]=size(data);
0047 mintime=0;
0048 if nargin < 3; [mintime,maxtime]=minmaxsptimes(data);clear mintime;
0049 else maxtime=T; end;
0050 if nargin < 4; params=[]; end;
0051 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
0052 clear err trialave params
0053 Nwin=round(Fs*win); % number of samples in window
0054 nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
0055 [f,findx]=getfgrid(Fs,nfft,fpass); 
0056 tapers=dpsschk(tapers,Nwin,Fs); % check tapers
0057 twin=linspace(0,win,Nwin); % times of occurrence of "samples" within window - times of evaluation of tapers
0058 
0059 Sc=zeros(length(findx),C,C);
0060 tn=mintime:win:maxtime-win;
0061 Nwins=length(tn);
0062 if d==3, % If there are multiple trials
0063 for iwin=1:Nwins,
0064     t=[tn(iwin) tn(iwin)+T];
0065     for i=1:Ntr, 
0066         data1=data(:,i);
0067         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
0068         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
0069         for k=1:C,
0070             for l=1:C,
0071                 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 
0072             Sc(:,k,l)=Sc(:,k,l)+spec;
0073             end
0074         end
0075     end
0076 end
0077 Sc=Sc/(Nwins*Ntr);
0078 end
0079 
0080 if d==2, % only one trial
0081 for iwin=1:Nwins,
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             Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
0088             end
0089         end
0090 end
0091 Sc=Sc/Nwins;
0092 end
0093 
0094 Cmat=Sc;
0095 Sdiag=zeros(length(findx),C);
0096 for k=1:C,
0097     Sdiag(:,k)=squeeze(Sc(:,k,k));
0098 end
0099 
0100 for k=1:C,
0101     for l=1:C,
0102         Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
0103     end
0104 end
0105 
0106 Ctot=zeros(length(findx),1); Cent=Ctot;
0107 Cvec=zeros(length(findx),C);
0108 for i=1:length(findx),
0109     [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
0110     Ctot(i)=s(1).^2/sum(s.^2); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);             
0111     Cvec(i,:)=transpose(u(:,1));
0112 
0113 end
0114

Generated on Mon 09-Oct-2006 00:54:52 by m2html © 2003