diff userProgramsTim/runningACF.m @ 38:c2204b18f4a2 tip

End nov big change

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