function [raster,psth,isih] = rasterize( data, range)
% Raster: make a binary matrix from event data
% data: Lx1 timing data file
% range: Mx2: [start end] time of M trials
% raster: DxM raster matrix
% psth: Dx1 PSTH
% isih: isi histogram

trials = size( range, 1);
duration = range(1,2) - range(1,1);	% trial duration

raster = zeros( duration, trials);
isi = [];

for i = 1:trials
	% spike timing
	tsp = data(range(i,1)<=data & data<range(i,2)) - range(i,1);
	raster(1+tsp,i) = 1;
	isi = [isi; tsp(2:end)-tsp(1:end-1)];
end
% imagesc( raster);

% psth
% psth = sum( raster, 2);	% spike counts
psth = mean( raster, 2);	% normalized
% bar( psth);

isih = hist( isi, 0:duration-1);