


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