rmeddis@38: function sampledacf = runningACF(inputSignalMatrix, sfreq, params)
rmeddis@38: % runningACF computes within-channel, running autocorrelations (acfs)
rmeddis@38: %  and samples it at given time steps
rmeddis@38: %
rmeddis@38: % INPUT
rmeddis@38: %  	inputSignalMatrix: a matrix (channel x time) of AN activity
rmeddis@38: %  	sfreq: sampling frequency
rmeddis@38: %   
rmeddis@38: % params: 	a list of parmeters to guide the computation:
rmeddis@38: %   params.lags: an array of lags to be computed (seconds)
rmeddis@38: %   params.acfTau: time constant (sec) of the running acf calculations
rmeddis@38: %     
rmeddis@38: %
rmeddis@38: %
rmeddis@38: % OUTPUT
rmeddis@38: %  	sampledacf:  	(time x lags)
rmeddis@38: 
rmeddis@38: %%
rmeddis@38: boundaryValue=[];
rmeddis@38: [nChannels inputLength]= size(inputSignalMatrix);
rmeddis@38: % list of BFs must be repeated is many fiber types used
rmeddis@38: %BFlist=method.nonlinCF;
rmeddis@38: %nfibertypes=nChannels/length(BFlist);
rmeddis@38: %BFlist=repmat(BFlist',2,1)';
rmeddis@38: 
rmeddis@38: dt=1/sfreq;
rmeddis@38: 
rmeddis@38: samplelength=0.003; %sample length in ms
rmeddis@38: gridpoints=round([samplelength:samplelength:inputLength*dt]/dt);
rmeddis@38: gridpoints=[1 gridpoints];
rmeddis@38: 
rmeddis@38: lags=params.lags;
rmeddis@38: 
rmeddis@38: nLags=length(lags);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: % Establish matrix of lag time constants 
rmeddis@38: %   these are all the same if tau<1
rmeddis@38: 
rmeddis@38: if params.acfTau<1          % all taus are the same
rmeddis@38:     acfTaus=repmat(params.acfTau, 1, nLags);
rmeddis@38:     acfDecayFactors=ones(size(lags)).*exp(-dt/params.acfTau);
rmeddis@38: else                     
rmeddis@38:     error('acfTau must be <1');
rmeddis@38: end
rmeddis@38: % make acfDecayFactors into a (channels x lags) matrix for speedy computation
rmeddis@38: acfDecayFactors=repmat(acfDecayFactors,nChannels, 1);
rmeddis@38: 
rmeddis@38: % P function lowpass filter decay (only one value needed)
rmeddis@38: %pDecayFactor=exp(-dt/params.lambda);
rmeddis@38: 
rmeddis@38: % ACF
rmeddis@38: % lagPointers is a list of pointers relative to 'time now'
rmeddis@38: lagPointers=round(lags/dt);
rmeddis@38: lagWeights=ones(nChannels,nLags);
rmeddis@38: if max(lagPointers)+1>inputLength
rmeddis@38:     error([' filteredSACF: not enough signal to evaluate ACF. Max(lag)= ' num2str(max(lags))])
rmeddis@38: end
rmeddis@38: 
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     acf=zeros(nChannels,nLags);
rmeddis@38:     sampledacf = zeros(length(gridpoints),nChannels,nLags);
rmeddis@38:     % create a runup buffer of signal
rmeddis@38:     buffer= zeros(nChannels, max(lagPointers));
rmeddis@38: 
rmeddis@38: inputSignalMatrix=[buffer inputSignalMatrix];
rmeddis@38: [nChannels inputLength]= size(inputSignalMatrix);
rmeddis@38: 
rmeddis@38: timeCounter=0; biggestSACF=0;
rmeddis@38: gridcounter =1;
rmeddis@38: for timePointer= max(lagPointers)+1:inputLength
rmeddis@38:     % acf is a continuously changing channels x lags matrix
rmeddis@38:     %   Only the current value is stored
rmeddis@38:     % sacf is the vertical summary of acf ( a vector) and all values are kept and returned
rmeddis@38:     % P is the smoothed version of sacf and all values are kept and returned
rmeddis@38:     % lagWeights emphasise some BF/lag combinations and ignore others
rmeddis@38:     % NB time now begins at the longest lag.
rmeddis@38:     % E.g. if max(lags) is .04 then this is when the ACf will begin.
rmeddis@38:     %            AN                                       AN delayed                             weights               filtering
rmeddis@38:     
rmeddis@38:     % This is the ACF calculation
rmeddis@38:     timeCounter=timeCounter+1;
rmeddis@38:     acf = (repmat(inputSignalMatrix(:, timePointer), 1, nLags) .* inputSignalMatrix(:, timePointer-lagPointers)).*lagWeights *dt + acf.* acfDecayFactors;
rmeddis@38:     
rmeddis@38:     %just sample on certain samplepoints
rmeddis@38:     if ~isempty(find(timeCounter==gridpoints))
rmeddis@38:         sampledacf(gridcounter,:,:)=acf;
rmeddis@38:         gridcounter = gridcounter+1;
rmeddis@38:     end
rmeddis@38:     
rmeddis@38:     end
rmeddis@38: end
rmeddis@38: 
rmeddis@38: 
rmeddis@38: 
rmeddis@38: