CrossSpecMatpt.m 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131
  1. function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
  2. %
  3. %
  4. % Multi-taper cross-spectral matrix - another routine, this one allows for multiple trials and channels
  5. % but does not do confidence intervals. Also this routine always averages
  6. % over trials - point process as times
  7. %
  8. % Usage:
  9. %
  10. % [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
  11. % Input:
  12. % Note units have to be consistent. See chronux.m for more information.
  13. % data (as a struct array with dimensions channels x trials) - note
  14. % that times of measurement have to be consistent, we assume all
  15. % times are specified relative to the start time of the trials which
  16. % are taken to be zero.
  17. % win (duration of non-overlapping window)
  18. % trialduration (since it is not possible to infer trial duration
  19. % from spike times, this is an optional argument. If not specified
  20. % the routine uses the minimum and maximum spike time (across all
  21. % channels and trials) as the window of calculation.) -
  22. % optional
  23. % params: structure with fields tapers, pad, Fs, fpass
  24. % - optional
  25. % tapers : precalculated tapers from dpss or in the one of the following
  26. % forms:
  27. % (1) A numeric vector [TW K] where TW is the
  28. % time-bandwidth product and K is the number of
  29. % tapers to be used (less than or equal to
  30. % 2TW-1).
  31. % (2) A numeric vector [W T p] where W is the
  32. % bandwidth, T is the duration of the data and p
  33. % is an integer such that 2TW-p tapers are used. In
  34. % this form there is no default i.e. to specify
  35. % the bandwidth, you have to specify T and p as
  36. % well. Note that the units of W and T have to be
  37. % consistent: if W is in Hz, T must be in seconds
  38. % and vice versa. Note that these units must also
  39. % be consistent with the units of params.Fs: W can
  40. % be in Hz if and only if params.Fs is in Hz.
  41. % The default is to use form 1 with TW=3 and K=5
  42. %
  43. % pad (padding factor for the FFT) - optional (can take values -1,0,1,2...).
  44. % -1 corresponds to no padding, 0 corresponds to padding
  45. % to the next highest power of 2 etc.
  46. % e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
  47. % to 512 points, if pad=1, we pad to 1024 points etc.
  48. % Defaults to 0.
  49. % Fs (sampling frequency) - optional. Default 1.
  50. % fpass (frequency band to be used in the calculation in the form
  51. % [fmin fmax])- optional.
  52. % Default all frequencies between 0 and Fs/2
  53. % Output:
  54. % Sc (cross spectral matrix frequency x channels x channels)
  55. % Cmat Coherence matrix frequency x channels x channels
  56. % Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
  57. % Cvec leading Eigenvector (frequency x channels)
  58. % Cent A different measure of total coherence: GM/AM of SV^2s
  59. % f (frequencies)
  60. d=ndims(data);
  61. if size(d,1)==1; error('Need multiple channels; are you sure your format is channels x trials ?');end;
  62. [C,Ntr]=size(data);
  63. mintime=0;
  64. if nargin < 3; [mintime,maxtime]=minmaxsptimes(data);clear mintime;
  65. else maxtime=T; end;
  66. if nargin < 4; params=[]; end;
  67. [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
  68. clear err trialave params
  69. Nwin=round(Fs*win); % number of samples in window
  70. nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
  71. [f,findx]=getfgrid(Fs,nfft,fpass);
  72. tapers=dpsschk(tapers,Nwin,Fs); % check tapers
  73. twin=linspace(0,win,Nwin); % times of occurrence of "samples" within window - times of evaluation of tapers
  74. Sc=zeros(length(findx),C,C);
  75. tn=mintime:win:maxtime-win;
  76. Nwins=length(tn);
  77. if d==3, % If there are multiple trials
  78. for iwin=1:Nwins,
  79. t=[tn(iwin) tn(iwin)+T];
  80. for i=1:Ntr,
  81. data1=data(:,i);
  82. data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
  83. J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
  84. for k=1:C,
  85. for l=1:C,
  86. spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
  87. Sc(:,k,l)=Sc(:,k,l)+spec;
  88. end
  89. end
  90. end
  91. end
  92. Sc=Sc/(Nwins*Ntr);
  93. end
  94. if d==2, % only one trial
  95. for iwin=1:Nwins,
  96. data1=data(:,i);
  97. data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
  98. J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
  99. for k=1:C,
  100. for l=1:C,
  101. Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
  102. end
  103. end
  104. end
  105. Sc=Sc/Nwins;
  106. end
  107. Cmat=Sc;
  108. Sdiag=zeros(length(findx),C);
  109. for k=1:C,
  110. Sdiag(:,k)=squeeze(Sc(:,k,k));
  111. end
  112. for k=1:C,
  113. for l=1:C,
  114. Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
  115. end
  116. end
  117. Ctot=zeros(length(findx),1); Cent=Ctot;
  118. Cvec=zeros(length(findx),C);
  119. for i=1:length(findx),
  120. [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
  121. % Ctot(i)=s(1)/sum(s); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);
  122. Ctot(i)=s(1)/sum(s); Cent(i)=exp(mean(log(s)))/mean(s);
  123. Cvec(i,:)=transpose(u(:,1));
  124. end