%COMPUTECSD	Fast vectorized CSD computation. 
%	[dataCSD] = computeCSD( data )
%	Uses the following kernel to filter
%	low pass in the spatial domain before
%	computing CSDs.
%	    1
%	 1  4  1
%	    1
%	Uses the following kernel to compute CSDs
%	   -1
%	-1  4 -1
%	   -1
%	Asumes current is zero accross the external 
%	boundaries of the electrode array.


function [dataCSD] = computeCSD( data );

dataFilt = BudapestbandPass(data,1000,5,200);

dataCh = dataFilt(:,2:65);  % create copy of data without the timestamp column


%--- Apply low pass spatial filter ---

dataChSoft = zeros( size( dataCh ) );

for ch = 1:64
   
   nNeighbors = 0;
   chCol = mod( ch - 1, 8 ) + 1;
   
   if( chCol - 1 >= 1 )
      dataChSoft (:, ch) = dataChSoft(:, ch) + dataCh(:, ch - 1);
      nNeighbors = nNeighbors + 1;
   end
   
   if( chCol + 1 <= 8 )
      dataChSoft(:, ch) = dataChSoft(:, ch) + dataCh(:, ch + 1);
      nNeighbors = nNeighbors + 1;
   end
   
   if( ch - 8 >= 1 )
      dataChSoft(:,ch) = dataChSoft(:, ch) + dataCh(:,ch - 8);
      nNeighbors = nNeighbors + 1;
   end
   
   if( ch + 8 <= 64 )
      dataChSoft(:,ch) = dataChSoft(:, ch) + dataCh(:,ch + 8);
      nNeighbors = nNeighbors + 1;
   end
   
   dataChSoft(:, ch) = dataChSoft(:, ch) + 4 * dataCh(:, ch);
   dataChSoft(:, ch) = dataChSoft(:, ch) / (nNeighbors + 4);
   
end


%--- Compute CSD ---

dataChCSD = zeros( size( dataCh ) );

for ch = 1:64
   
   nNeighbors = 0;
   chCol = mod( ch - 1, 8 ) + 1;
   
   if( chCol - 1 >= 1 )
      dataChCSD(:, ch) = dataChCSD(:, ch) - dataChSoft(:, ch - 1);
      nNeighbors = nNeighbors + 1;
   end
   
   if( chCol + 1 <= 8 )
      dataChCSD(:, ch) = dataChCSD(:, ch) - dataChSoft(:, ch + 1);
      nNeighbors = nNeighbors + 1;
   end
   
   if( ch - 8 >= 1 )
      dataChCSD(:,ch) = dataChCSD(:, ch) - dataChSoft(:,ch - 8);
      nNeighbors = nNeighbors + 1;
   end
   
   if( ch + 8 <= 64 )
      dataChCSD(:,ch) = dataChCSD(:, ch) - dataChSoft(:,ch + 8);
      nNeighbors = nNeighbors + 1;
   end
   
   dataChCSD(:, ch) = dataChCSD(:, ch) + nNeighbors * dataChSoft(:, ch);
   
end


dataCSD = [dataFilt(:,1),dataChCSD];  % paste the original time stamp column to the computed CSDs.

t=dataFilt(1:end,1);

   figure
   for n=1:64
   subplot (8,8,n); plot (t,(); axis([0 100 -5 5]); axis off;
   end