function s = SpikeJitter(N,t,L,Lmethod,Gmethod)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function s = SpikeJitter(N,t,L,Lmethod,Gmethod)
%
% DESCRIPTION OF JITTER:
%
% Independently and uniformly jitters the spike times in t over jitter
% windows of length L.  The UNSORTED jittered spike train is stored in 
% s(:,1).  This is repeated N times with the kth jittered spike train 
% stored in s(:,k).  sort(s) gives the sorted jittered spike trains.
%
% DESCRIPTION OF INPUT:
%
% N is the the number of MonteCarlo trials.
% t is a sequence of spike times.
% L is a non-negative scalar (duration of jitter window).
% Lmethod is one of {'centered','fixed'}, 'fixed' is default
%   'centered': the jitter window for t(k) is t(k) + [-L/2,L/2].
%   'fixed': the jitter window for t(k) is floor(t(k)/L)*L + [0,L].
% Gmethod is one of {'real','integer'}, 'integer' is default
%   'real': no changes from above
%   'integer': t and jittered spikes are constrained to the integers
%       L is assumed to be an integer >= 1
%       the jitter positions are
%       'centered': t(k)-floor(L/2) + {0,1,...,L-1}
%       'fixed': floor(t(k)/L)*L + {0,1,...,L-1}
%
% Tested on Matlab 7.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matt Harrison 
% 3/20/2006

% Copyright (C) 2008 Matthew T Harrison
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    See <http://www.gnu.org/licenses/> for a copy of the license.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if numel(N) ~= 1 || N ~= round(N) || N < 0 || isinf(N)
    error('bad N')
end
if numel(L) ~= 1 || L < 0
    error('bad L')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the methods and 
% compute the appropriate intervals.
% More input checking.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 5, Gmethod = 'integer'; end
if nargin < 4, Lmethod = 'fixed'; end

switch lower(Lmethod)

    case 'centered'

        switch lower(Gmethod)

            case 'real'
                Greal = true;
                u = t(:) - L/2;
            case 'integer'
                Greal = false;
                if L ~= round(L) | L < 1, error('bad L'), end
                if any(t ~= round(t)), error('non integer t'), end
                u = round(t(:) - floor(L/2));
            otherwise
                error('unknown Gmethod')
        end
        
    case 'fixed'

        switch lower(Gmethod)
            
            case 'real'
                Greal = true;
                if L > 0
                    u = floor(t(:)/L)*L;
                else
                    u = t(:);
                end
            case 'integer'
                Greal = false;
                if L ~= round(L) | L < 1, error('bad L'), end
                if any(t ~= round(t)), error('non integer t'), end
                u = round(floor(t(:)/L)*L);
            otherwise
                error('unknown Gmethod')
        end
        
    otherwise
        error('unknown Lmethod')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JITTER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m = numel(t);

s = repmat(u,1,N)+L*rand(m,N);
if ~Greal
    s = floor(s);
end
