Home > chronux > spectral_analysis > pointtimes > CrossSpecMatpt.m

# CrossSpecMatpt

## 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:
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...).
to the next highest power of 2 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:
• extractdatapt Extract segements of spike times between t(1) and t(2)
• minmaxsptimes Find the minimum and maximum of the spike times in each channel
• mtfftpt Multi-taper fourier transform for point process given as times
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...).
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;
0068 clear err trialave params
0069 Nwin=round(Fs*win); % number of samples in window
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     Ctot(i)=s(1)/sum(s); Cent(i)=exp(mean(log(s)))/mean(s);
0128     Cvec(i,:)=transpose(u(:,1));
0129
0130 end
0131```

Generated on Fri 28-Sep-2012 12:34:30 by m2html © 2005