function x = lyngby_design_events2design(e, TR, L, varargin)

% lyngby_design_events2design - Convert events to design matrix column
%
%      function x = lyngby_design_events2design(e, TR, L, ...
%         'PropertyName', 'PropertyValue')
%
%      Input:  e       Events (in seconds)
%              TR      Sampling time (in seconds)
%              L       Length (number of scans)
%              Property:
%                Decimation     [ {FourierZero} |
%                               FirstOrderInterpolation ] Type of
%                               decimation
%                Oversampling   {100*TR} Number of oversampling
%                               compared to TR. Default is TR*10 that
%                               corresponds to 10ms
%
%      Output: x       Vector corresponding to decimated paradigm
%
%      This function will take a set of events denoted by the time of
%      occurence in eg. seconds are convert them to a vector that is
%      suitable as a design variable or "paradigm" input to a filter. 
%
%      The first element in x, ie. x(1), corresponds to zero
%      seconds. The second element x(2) to the time at TR. x is length L. 
%
%      The input 'e' should be a column vector or a matrix with size
%      Nx2. If 'e' is a vector then the times in 'e' is taken as
%      "Dirac"-events and if 'e' is a matrix then the first column
%      should describe on-times while the second column should
%      describe off-times.
%
%      Example: 
%        e = [ 20 31.56 45.20 52.4 70.1 ]';
%        x = lyngby_design_events2design(e, 2, 50)
%        plot(0:2:98, x, '-x')
%
%      See also LYNGBY.
%
% $Id: lyngby_design_events2design.m,v 1.4 2002/05/19 22:25:30 fnielsen Exp $



    % Check arguments
    if nargin < 3
      error('Not enough arguments')
    end
    if ~any(size(e,2) == [ 1 2 ])
      error('the number of columns in ''e'' should be 1 or 2')
    end


    % Default values
    oversampling = round(100 * TR);
    decimation = 'fourierzero';

    
    % Parse Properties
    n = 1;
    while n <= length(varargin)
      arg = lower(varargin{n});

      if strcmp(arg, 'oversampling')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg)
	  if length(arg) == 1
	    if arg > 0
	      oversampling = arg;
	    else
	      error(sprintf(['''Oversampling'' should be larger than 0. ' ...
		    'It was: %d'], arg))
	    end
	  else
	    error(['The argument with ''Oversampling'' should be ' ...
		  'a scalar.']); 
	  end
	else
	  error('The argument with ''Oversampling'' should be a integer.'); 
	end

      elseif strcmp(arg, 'decimation')
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg)
	  if strcmp(arg, 'fourierzero') | strcmp(arg, 'firstorderinterpolation')
	    decimation = arg;
	  else
	    error(sprintf(['Argument to ''Decimation'' should be' ...
		  ' ''FourierZero'' or ''FirstOrderInterpolation'' ']))
	  end
	else
	  error(['The argument with ''Decimation'' should be ' ...
		'a string.']); 
	end
	      
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    
    
    
    % Length of oversampled sequence
    L_os = L * oversampling;
    
    % Output vector at oversampled frequency 
    x_os = zeros(ceil(L_os),1); 

      
    if size(e,2) == 2
      % Block events
      
      on  = e(:,1);
      off = e(:,2);
      
      for n = 1:length(on)
	if on(n) > off(n)
	  warning(sprintf([ ...
		' ''off'' event before an ''on'' event: ' ...
		'e(%d) = %f - %f' ...
		], n, on(n), off(n)));
	end
	non_os = round(on(n) * oversampling / TR);
	noff_os = round(off(n) * oversampling / TR);
	if non_os > L_os | noff_os > L_os
	  warning( ...
	      sprintf('Event outside window: %f - %f', on(n), off(n)))
	elseif non_os < 0 | noff_os < 0
	  warning( ...
	      sprintf('Ignoring event: %f - %f', on(n), off(n)))
	else
	  x_os( (non_os:noff_os)+1 ) = x_os( (non_os:noff_os)+1 ) + 1; 
	end

      end
      
      normalize = 1;
      
    elseif size(e,2) == 1
      % "Dirac"-events

      % Insert a one for each event
      % x_os(1) is at zero time
      for n = 1:length(e)
	n_os = round(e(n) * oversampling / TR);
	if n_os > L_os
	  warning(sprintf('Event outside window: %f', e(n)))
	elseif n_os < 0
	  warning(sprintf('Ignoring event before zero: %f', e(n)))
	else
	  x_os(n_os+1) = x_os(n_os+1) + 1; 
	end
      end
      
      normalize = oversampling;

    end

    
    if strcmp(decimation, 'fourierzero')

      % Forward Fourier transform
      f = fft(x_os);
      
      % Filter in Fourier domain: Zero high frequency components
      f(1+ceil(L/2):end-floor(L/2)) = 0;
      
      % Inverse Fourier transform with imaginary numerical errors
      % discarded
      x0_os = real(ifft(f));
      
      % Decimate
      x = x0_os(round(1:oversampling:L_os)) * normalize;
      
    elseif strcmp(decimation, 'firstorderinterpolation')
      
      x = sum(reshape(x_os, oversampling, L)) / oversampling * normalize;
      
    else
      error('Internal error')
    end









