specerr

 Function to compute lower and upper confidence intervals on the spectrum 
 Usage: Serr=specerr(S,J,err,trialave,numsp)
 Outputs: Serr (Serr(1,...) - lower confidence level, Serr(2,...) upper confidence level)

 Inputs:
 S - spectrum
 J - 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
 numsp    - number of spikes in each channel. specify only when finite
            size correction required (and of course, only for point
            process data)

 Outputs:
 Serr - error estimates. Only for err(1)>=1. If err=[1 p] or [2 p] Serr(...,1) and Serr(...,2)
 contain the lower and upper error bars with the specified method.

0001 function Serr=specerr(S,J,err,trialave,numsp)
0002 % Function to compute lower and upper confidence intervals on the spectrum
0003 % Usage: Serr=specerr(S,J,err,trialave,numsp)
0004 % Outputs: Serr (Serr(1,...) - lower confidence level, Serr(2,...) upper confidence level)
0005 %
0006 % Inputs:
0007 % S - spectrum
0008 % J - tapered fourier transforms
0009 % err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates;
0010 %                   p - p value for error estimates)
0011 % trialave - 0: no averaging over trials/channels
0012 %            1 : perform trial averaging
0013 % numsp    - number of spikes in each channel. specify only when finite
0014 %            size correction required (and of course, only for point
0015 %            process data)
0016 %
0017 % Outputs:
0018 % Serr - error estimates. Only for err(1)>=1. If err=[1 p] or [2 p] Serr(...,1) and Serr(...,2)
0019 % contain the lower and upper error bars with the specified method.
0020 if nargin < 4; error('Need at least 4 input arguments'); end;
0021 if err(1)==0; error('Need err=[1 p] or [2 p] for error bar calculation. Make sure you are not asking for the output of Serr'); end;
0022 [nf,K,C]=size(J);
0023 errchk=err(1);
0024 p=err(2);
0025 pp=1-p/2;
0026 qq=1-pp;
0027 
0028 if trialave
0029    dim=K*C;
0030    C=1;
0031    dof=2*dim;
0032    if nargin==5; dof = fix(1/(1/dof + 1/(2*sum(numsp)))); end
0033    J=reshape(J,nf,dim);
0034 else
0035    dim=K;
0036    dof=2*dim*ones(1,C);
0037    for ch=1:C;
0038      if nargin==5; dof(ch) = fix(1/(1/dof + 1/(2*numsp(ch)))); end 
0039    end;
0040 end;
0041 Serr=zeros(2,nf,C);
0042 if errchk==1;
0043    Qp=chi2inv(pp,dof);
0044    Qq=chi2inv(qq,dof);
0045    Serr(1,:,:)=dof(ones(nf,1),:).*S./Qp(ones(nf,1),:);
0046    Serr(2,:,:)=dof(ones(nf,1),:).*S./Qq(ones(nf,1),:);
0047 elseif errchk==2;
0048    tcrit=tinv(pp,dim-1);
0049    for k=1:dim;
0050        indices=setdiff(1:dim,k);
0051        Jjk=J(:,indices,:); % 1-drop projection
0052        eJjk=squeeze(sum(Jjk.*conj(Jjk),2));
0053        Sjk(k,:,:)=eJjk/(dim-1); % 1-drop spectrum
0054    end;
0055    sigma=sqrt(dim-1)*squeeze(std(log(Sjk),1,1)); if C==1; sigma=sigma'; end; 
0056    conf=repmat(tcrit,nf,C).*sigma;
0057    conf=squeeze(conf); 
0058    Serr(1,:,:)=S.*exp(-conf); Serr(2,:,:)=S.*exp(conf);
0059 end;
0060 Serr=squeeze(Serr);