From 2a4674180643db0992620290d0b1bf73bb7ec561 Mon Sep 17 00:00:00 2001
From: Jan Grewe <jan.grewe@g-node.org>
Date: Wed, 15 Oct 2014 18:42:55 +0200
Subject: [PATCH] nearly done

---
 programming/exercises/psths.m | 56 +++++++++++++++++++++++++++++++----
 1 file changed, 50 insertions(+), 6 deletions(-)

diff --git a/programming/exercises/psths.m b/programming/exercises/psths.m
index 7d76c86..f6b89e9 100644
--- a/programming/exercises/psths.m
+++ b/programming/exercises/psths.m
@@ -1,16 +1,60 @@
 load('ampullary.mat')
+sample_rate = 20000; % Hz
+max_time = 0;
+for i = 1:size(times,2)
+    max_time = max([max_time, max(times{i})]);
+end
 fig = figure();
-
+set(gcf,'Color', 'white')
 %% create PSTH on the basis of the interspike intervals
-fig.sub
+subplot(3,1,1)
+hold on
 % 1. get the interspike intervals for each trial
 for i = 1:size(times,2)
-   isi = diff(times{i});
-    
+    t = times{i};
+    isi = diff(t);
+    plot(t(2:end), 1./isi)
 end
 
-%% create PSTH using the binning method
+xlabel('time [s]')
+ylabel('firing rate [Hz]')
+box('off')
+title('instanataneous firing rate')
 
+%% create PSTH using the binning method
+subplot(3,1,2)
+box('off')
+bin_width = 0.02; % s
+edges = 0:bin_width:max_time;
+firing_rate = [];
+for i = 1:size(times,2)
+    t = times{i};
+    [n, t] = hist(t, edges);
+    if isempty(firing_rate)
+        firing_rate = n / bin_width;
+    else
+        firing_rate = firing_rate + (n / bin_width / size(times,2));
+    end
+end
+plot(t,firing_rate)
+xlabel('time [s]')
+ylabel('firing rate [Hz]')
+title('binning method')
 
+%% create PSTH using the kernel-convolution method
+subplot(3,1,3)
+binary_spikes = zeros(size(times,2), round(max_time*sample_rate));
+resps = zeros(size(binary_spikes));
+window = hann(bin_width/4*sample_rate,'symmetric');
+window = window/sum(window);
 
-%% create PSTH using the kernel-convolution method
\ No newline at end of file
+for i = 1:size(times,2)
+    t = times{i};
+    temp = round(t*sample_rate);
+    if temp(1) <= 0 
+      temp(1) = 1;
+    end
+    binary_spikes(i, temp) = 1;
+    resps(i,:) = conv(binary_spikes(i,:), window, 'same');
+end
+plot((0:1/sample_rate:max_time), mean(resps,2))