


Function to compute lower and upper confidence intervals on the coherency given the tapered fourier transforms,
errchk, trialave.
Usage: [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2)
Inputs:
C - coherence
J1,J2 - tapered fourier transforms
err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates;
p - p value for error estimates)
trialave - 0: no averaging over trials/channels
1 : perform trial averaging
numsp1 - number of spikes for data1. supply only if finite size corrections are required
numsp2 - number of spikes for data2. supply only if finite size corrections are required
Outputs:
confC - confidence level for C - only for err(1)>=1
phistd - theoretical or jackknife standard deviation for phi for err(1)=1 and err(1)=2 respectively
returns zero if coherence is 1
Cerr - Jacknife error bars for C - only for err(1)=2

0001 function [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2) 0002 % Function to compute lower and upper confidence intervals on the coherency given the tapered fourier transforms, 0003 % errchk, trialave. 0004 % Usage: [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2) 0005 % Inputs: 0006 % C - coherence 0007 % J1,J2 - tapered fourier transforms 0008 % err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates; 0009 % p - p value for error estimates) 0010 % trialave - 0: no averaging over trials/channels 0011 % 1 : perform trial averaging 0012 % numsp1 - number of spikes for data1. supply only if finite size corrections are required 0013 % numsp2 - number of spikes for data2. supply only if finite size corrections are required 0014 % 0015 % Outputs: 0016 % confC - confidence level for C - only for err(1)>=1 0017 % phistd - theoretical or jackknife standard deviation for phi for err(1)=1 and err(1)=2 respectively 0018 % returns zero if coherence is 1 0019 % Cerr - Jacknife error bars for C - only for err(1)=2 0020 0021 if nargin < 5; error('Need at least 5 input arguments'); end; 0022 if err(1)==0; error('Need err=[1 p] or [2 p] for error bar calculation'); end; 0023 if nargout==4 && err(1)==1; error('Cerr contains Jackknife errors: only computed when err(1) is 2'); end; 0024 [nf,K,Ch]=size(J1); 0025 errchk=err(1); 0026 p=err(2); 0027 pp=1-p/2; 0028 % 0029 % Find the number of degrees of freedom 0030 % 0031 if trialave; 0032 dim=K*Ch; 0033 dof=2*dim; 0034 dof1=dof; 0035 dof2=dof; 0036 Ch=1; 0037 if nargin>=6 && ~isempty(numsp1) 0038 totspikes1=sum(numsp1); 0039 dof1=fix(2*totspikes1*dof/(2*totspikes1+dof)); 0040 end 0041 if nargin==7 && ~isempty(numsp2); 0042 totspikes2=sum(numsp2); 0043 dof2=fix(2*totspikes2*dof/(2*totspikes2+dof)); 0044 end; 0045 dof=min(dof1,dof2); 0046 J1=reshape(J1,nf,dim); 0047 J2=reshape(J2,nf,dim); 0048 else 0049 dim=K; 0050 dof=2*dim; 0051 dof1=dof; 0052 dof2=dof; 0053 for ch=1:Ch; 0054 if nargin>=6 && ~isempty(numsp1); 0055 totspikes1=numsp1(ch); 0056 dof1=fix(2*totspikes1*dof/(2*totspikes1+dof)); 0057 end; 0058 if nargin==7 && ~isempty(numsp2); 0059 totspikes2=numsp2(ch); 0060 dof2=fix(2*totspikes2*dof/(2*totspikes2+dof)); 0061 end; 0062 dof(ch)=min(dof1,dof2); 0063 end; 0064 end; 0065 % 0066 % variance of the phase 0067 % 0068 % 0069 % Old code is the next few lines - new code is in the if statement below 0070 % beginning line 87 0071 % 0072 % if isempty(find((C-1).^2 < 10^-16)); 0073 % phierr = sqrt((2./dof(ones(nf,1),:)).*(1./(C.^2) - 1)); 0074 % else 0075 % phierr = zeros(nf,Ch); 0076 % end 0077 0078 % 0079 % theoretical, asymptotic confidence level 0080 % 0081 if dof <= 2 0082 confC = 1; 0083 else 0084 df = 1./((dof/2)-1); 0085 confC = sqrt(1 - p.^df); 0086 end; 0087 % 0088 % Phase standard deviation (theoretical and jackknife) and jackknife 0089 % confidence intervals for C 0090 % 0091 if errchk==1; 0092 totnum=nf*Ch; 0093 phistd=zeros(totnum,1); 0094 CC=reshape(C,[totnum,1]); 0095 indx=find(abs(CC-1)>=1.e-16); 0096 dof=repmat(dof,[nf,1]); 0097 dof=reshape(dof,[totnum 1]); 0098 phistd(indx)= sqrt((2./dof(indx).*(1./(C(indx).^2) - 1))); 0099 phistd=reshape(phistd,[nf Ch]); 0100 elseif errchk==2; 0101 tcrit=tinv(pp,dof-1); 0102 for k=1:dim; 0103 indxk=setdiff(1:dim,k); 0104 J1k=J1(:,indxk,:); 0105 J2k=J2(:,indxk,:); 0106 eJ1k=squeeze(sum(J1k.*conj(J1k),2)); 0107 eJ2k=squeeze(sum(J2k.*conj(J2k),2)); 0108 eJ12k=squeeze(sum(conj(J1k).*J2k,2)); 0109 Cxyk=eJ12k./sqrt(eJ1k.*eJ2k); 0110 absCxyk=abs(Cxyk); 0111 atanhCxyk(k,:,:)=sqrt(2*dim-2)*atanh(absCxyk); 0112 phasefactorxyk(k,:,:)=Cxyk./absCxyk; 0113 % indxk=setdiff(1:dim,k); 0114 % J1jk=J1(:,indxk,:); 0115 % J2jk=J2(:,indxk,:); 0116 % eJ1jk=squeeze(sum(J1jk.*conj(J1jk),2)); 0117 % eJ2jk=squeeze(sum(J2jk.*conj(J2jk),2)); 0118 % eJ12jk=squeeze(sum(conj(J1jk).*J2jk,2)); 0119 % atanhCxyjk(k,:,:)=sqrt(2*dim-2)*atanh(abs(eJ12jk)./sqrt(eJ1jk.*eJ2jk)); 0120 end; 0121 atanhC=sqrt(2*dim-2)*atanh(C); 0122 sigma12=sqrt(dim-1)*squeeze(std(atanhCxyk,1,1)); 0123 % sigma12=sqrt(dim-1)*squeeze(std(atanhCxyjk,1,1)); 0124 if Ch==1; sigma12=sigma12'; end; 0125 Cu=atanhC+tcrit(ones(nf,1),:).*sigma12; 0126 Cl=atanhC-tcrit(ones(nf,1),:).*sigma12; 0127 Cerr(1,:,:) = tanh(Cl/sqrt(2*dim-2)); 0128 Cerr(2,:,:) = tanh(Cu/sqrt(2*dim-2)); 0129 phistd=(2*dim-2)*(1-abs(squeeze(mean(phasefactorxyk)))); 0130 if trialave; phistd=phistd'; end; 0131 end 0132 % ncrit=norminv(pp); 0133 % phierr=zeros([2 size(phistd)]); 0134 % phierr(1,:,:)=phi-ncrit*phistd; phierr(2,:,:)=phi+ncrit*phistd;