% raster_demo.m  
% 
% script demonstrating basic plotting data format for frequency tuning
% curve data collected from a single A1 neuron
%
% what this is meant to teach:  basics of handling data and a few tips for
% generating pretty and interpretable plots
%
% recommended usage:  read code, copy and paste non-commented commands into
% the Matlab command window one line at a time.
%
% SVD 2014-09
%

% start with clean slate
clear
close all

% load the data
load raster_data

% list the relevant data
whos

% spike_times - cell array of spike times (in seconds) during each
% presentation of each tone.  size(spike_times)= 10 x 20 -- 10 reps of 20
% different tones
[trial_count,stim_count]=size(spike_times)

% list the frequency of each tone
tone_frequency

% notice that length of frequency list matches size(spike_times,2):
length(tone_frequency)

% when did the tone start during each trial?
StimOnsetTime

% when did the tone stop?
StimOffsetTime

% how long was each trial?
TrialDuration

% what happened on the seventh repetition of the sizth tone?  set indices:
trial_index=7;
stim_index=6;

% when were there spikes on this trial?
spike_times{trial_index,stim_index}

% how many spikes before tone onset?
sum(spike_times{trial_index,stim_index} <= StimOnsetTime)

% how many spikes after tone onset and before offset?
sum(spike_times{trial_index,stim_index} > StimOnsetTime & spike_times{trial_index,stim_index} <= StimOffsetTime)

% let's make a raster plot for all presentations of this stimulus:
figure;
subplot(3,1,1);
for trial_index=1:trial_count,
    for spike_index=1:length(spike_times{trial_index,stim_index}),
       plot(spike_times{trial_index,stim_index}(spike_index),trial_index,'k.');
       hold on
       
       % make the axes fit nicely
       axis([0 TrialDuration 0 trial_count+1]);
       pause(0.05);
    end
end
% delineate the tone onset/offset times
plot([StimOnsetTime StimOnsetTime],[0 trial_count+1],'k--');
plot([StimOffsetTime StimOffsetTime],[0 trial_count+1],'k--');


% add some labels
xlabel('time (sec)');
ylabel('trial number');
title_string=sprintf('Response raster for %d Hz',tone_frequency(stim_index));
title(title_string);

% now plot the PSTH response for this tone:
subplot(3,1,2);

% how many seconds per bin - 10 ms
psth_bin_size=0.01;

% initialize the psth as a vector of zeros
psth=zeros(ceil(TrialDuration./psth_bin_size),1);

% figure out center of each time bin in the PSTH
time_axis=(psth_bin_size./2):psth_bin_size:TrialDuration;

% find each spike time and increment apprpriate bin:
for trial_index=1:trial_count,
    spike_count_this_trial=length(spike_times{trial_index,stim_index});
    for spike_index=1:spike_count_this_trial,
        psth_bin=ceil(spike_times{trial_index,stim_index}(spike_index)./psth_bin_size);
        psth(psth_bin)=psth(psth_bin)+1./trial_count;
    end
    plot(time_axis,psth);
    title(sprintf('PSTH after %d trials',trial_index));
    pause(0.1);
end

% delineate the tone onset/offset times
axis_height=max(psth)*1.1;
hold on
plot([StimOnsetTime StimOnsetTime],[0 axis_height],'k--');
plot([StimOffsetTime StimOffsetTime],[0 axis_height],'k--');

axis([0 TrialDuration 0 axis_height]);

xlabel('time (sec)');
ylabel('PSTH (spikes/sec)');

% you can also plot the PSTH using bars
subplot(3,1,3);
bar(time_axis,psth,'k');

% delineate the tone onset/offset times
hold on
plot([StimOnsetTime StimOnsetTime],[0 axis_height],'k--');
plot([StimOffsetTime StimOffsetTime],[0 axis_height],'k--');
axis([0 TrialDuration 0 axis_height]);

xlabel('time (sec)');
ylabel('PSTH (spikes/sec)');
title('PSTH plotted using bars');
