mtfftpt

 Multi-taper fourier transform for point process given as times

 Usage:
 [J,Msp,Nsp]=mtfftpt (data,tapers,nfft,t,f,findx) - all arguments required
 Input: 
       data        (struct array of times with dimension channels/trials; 
                   also takes in 1d array of spike times as a column vector) 
       tapers      (precalculated tapers from dpss) 
       nfft        (length of padded data) 
       t           (time points at which tapers are calculated)
       f           (frequencies of evaluation)
       findx       (index corresponding to frequencies f) 
 Output:
       J (fft in form frequency index x taper index x channels/trials)
       Msp (number of spikes per sample in each channel)
       Nsp (number of spikes in each channel)

0001 function [J,Msp,Nsp]=mtfftpt(data,tapers,nfft,t,f,findx)
0002 % Multi-taper fourier transform for point process given as times
0003 %
0004 % Usage:
0005 % [J,Msp,Nsp]=mtfftpt (data,tapers,nfft,t,f,findx) - all arguments required
0006 % Input:
0007 %       data        (struct array of times with dimension channels/trials;
0008 %                   also takes in 1d array of spike times as a column vector)
0009 %       tapers      (precalculated tapers from dpss)
0010 %       nfft        (length of padded data)
0011 %       t           (time points at which tapers are calculated)
0012 %       f           (frequencies of evaluation)
0013 %       findx       (index corresponding to frequencies f)
0014 % Output:
0015 %       J (fft in form frequency index x taper index x channels/trials)
0016 %       Msp (number of spikes per sample in each channel)
0017 %       Nsp (number of spikes in each channel)
0018 if nargin < 6; error('Need all input arguments'); end;
0019 if isstruct(data); C=length(data); else C=1; end% number of channels
0020 K=size(tapers,2); % number of tapers
0021 nfreq=length(f); % number of frequencies
0022 if nfreq~=length(findx); error('frequency information (last two arguments) inconsistent'); end;
0023 H=fft(tapers,nfft,1);  % fft of tapers
0024 H=H(findx,:); % restrict fft of tapers to required frequencies
0025 w=2*pi*f; % angular frequencies at which ft is to be evaluated
0026 Nsp=zeros(1,C); Msp=zeros(1,C);
0027 for ch=1:C;
0028   if isstruct(data);
0029      fnames=fieldnames(data);
0030      eval(['dtmp=data(ch).' fnames{1} ';'])
0031      indx=find(dtmp>=min(t)&dtmp<=max(t));
0032      if ~isempty(indx); dtmp=dtmp(indx);
0033      end;
0034   else
0035      dtmp=data;
0036      indx=find(dtmp>=min(t)&dtmp<=max(t));
0037      if ~isempty(indx); dtmp=dtmp(indx);
0038      end;
0039   end;
0040   Nsp(ch)=length(dtmp);
0041   Msp(ch)=Nsp(ch)/length(t);
0042   if Msp(ch)~=0;
0043       data_proj=interp1(t',tapers,dtmp);
0044       exponential=exp(-i*w'*(dtmp-t(1))');
0045       J(:,:,ch)=exponential*data_proj-H*Msp(ch);
0046   else
0047       J(1:nfreq,1:K,ch)=0;
0048   end;
0049 end;