function out = PlotPSTH(data,range)
% Raster & PSTH Display プログラム
% figureハンドル = PlotPSTH(data, 試行範囲)

out = figure('Position',[0 0 600 600]);

timings=[]; 	   %%% 全試行のスパイクタイミング列

trial=size(range,1);
len=range(1,2)-range(1,1);

for i=1:trial
   train=data(range(i,1)<=data & data<range(i,2))-range(i,1);
			%%% 試行ごとに切り出し、時刻を揃える

   timings=[timings;train];			%%% スパイクタイミング列に追加

   %%%%%%%%%%%%%%%%%% ラスタ表示
   ha=subplot('Position',[0.1 (0.95-0.5*i/trial) 0.8 0.4/trial]); 
			%%% ラスタ表示の位置を指定

   hs=stem(train, ones(length(train)));		%%% X軸からの垂直線
   set(hs,'marker','none');
   set(hs,'color',[0 0 1]);
   axis off;
   set(ha,'xlim',[0 len]);	%%% 範囲：試行１の範囲
   set(ha,'box','off');
   set(ha,'yticklabel',[]);
   set(ha,'color',[0.8 0.8 0.8]);		%%% 背景を同じに
end

%%%%%%%%%%%%%%%%%% PSTH表示
ha=subplot('Position', [0.1 0.1 0.8 0.3]); 	%%% PSTH表示の位置を指定

stairs(0:50:len, [hist(timings, 25:50:(len-25) )/50*1000/trial 0]);
			%%% 50msビンでヒストグラムにして発火頻度に換算
			%%% 階段プロットで表示

set(ha,'xlim',[0 len]);	%%% 範囲：試行１の範囲
set(ha,'box','off');
set(ha,'color',[0.8 0.8 0.8]);
axis on;
xlabel('Time [ms]');
ylabel('Frequency [Hz]');