diff --git a/spike_trains/code/binned_rate.m b/spike_trains/code/binned_rate.m
new file mode 100644
index 0000000..cf3dcf0
--- /dev/null
+++ b/spike_trains/code/binned_rate.m
@@ -0,0 +1,16 @@
+function [time, rate] = binned_rate(spike_times, bin_width, dt, t_max)
+% Calculates the firing rate with the binning method. The hist funciton is
+% used to count the number of spikes in each bin. 
+% Function takes the spike times, the width of the bins in seconds as well
+% as the temporal resolution and the tiral duration.
+%
+% Returns two vectors containing the time and the rate.
+
+time = 0:dt:t_max-dt;
+bins = 0:bin_width:t_max;
+rate = zeros(size(time));
+
+h = hist(spike_times, bins) ./bin_width;
+for i = 2:length(bins)
+    rate(round(bins(i-1)/dt)+1:round(bins(i)/dt)) = h(i);
+end
diff --git a/spike_trains/code/convolution_rate.m b/spike_trains/code/convolution_rate.m
new file mode 100644
index 0000000..8f61b11
--- /dev/null
+++ b/spike_trains/code/convolution_rate.m
@@ -0,0 +1,22 @@
+function [time, rate] = convolution_rate(spike_times, sigma, dt, t_max)
+% Calculates the firing rate with the convolution method. Function takes
+% the spike times, the width of the Gaussian kernel (sigma, i.e. the
+% standard deviation), the temporal resolution, and the the duration of 
+% the trial. 
+%
+% Returns two vectors containing the time and the rate.
+
+time = 0:dt:t_max-dt;
+rate = zeros(size(time));
+spike_indices = round(spike_times/dt);
+rate(spike_indices) = 1;
+kernel = gauss_kernel(sigma, dt);
+
+rate = conv(rate, kernel, 'same');
+end
+
+
+function y = gauss_kernel(s, step)
+    x = -4*s:step:4*s;
+    y = exp(-0.5*(x./s).^2)/sqrt(2*pi)/s;
+end
diff --git a/spike_trains/code/instantaneous_rate.m b/spike_trains/code/instantaneous_rate.m
new file mode 100644
index 0000000..18b5fed
--- /dev/null
+++ b/spike_trains/code/instantaneous_rate.m
@@ -0,0 +1,22 @@
+function [time, rate] = instantaneous_rate(spike_times, dt, t_max)
+% Funttion takes the spike times of a single trial, the temporal
+% rersolution, dt, and the duration of the trial in seconds. 
+% Function calculates the firing rate as the inverse of the interspike
+% interval.
+%
+% Returns the vector representing time and a vector containing the rate.
+%
+
+
+time = 0:dt:t_max-dt;
+rate = zeros(size(time));
+
+isis = diff([0 spike_times]);
+inst_rate = 1./isis;
+spike_indices = [1 round(spike_times ./ dt)];
+
+for i = 2:length(spike_indices)
+    rate(spike_indices(i-1):spike_indices(i)) = inst_rate(i-1);
+end
+    
+
diff --git a/spike_trains/code/psths.m b/spike_trains/code/psths.m
new file mode 100644
index 0000000..c64110c
--- /dev/null
+++ b/spike_trains/code/psths.m
@@ -0,0 +1,86 @@
+%% Instantaneous rate
+clear 
+clc
+close all
+
+load 'lifoustim.mat'
+t_max = 30;
+rates = zeros(length(spikes), 30/dt);
+
+for i = 1:length(spikes)
+    [t, rates(i,:)] = instantaneous_rate(spikes{i}, dt, t_max);
+end
+
+
+f = figure();
+set(f, 'paperunits', 'centimeters')
+set(f, 'papersize', [10, 10])
+set(f, 'paperposition', [0 0 10 10])
+hold on
+plot(t, rates(1,:), 'displayname', 'trial 1')
+plot(t, mean(rates,1), 'displayname', 'average rate')
+
+xlabel('time [s]', 'fontsize', 10)
+ylabel('firing rate [Hz]','fontsize', 10)
+legend('show')
+box('off')
+
+saveas(gcf, 'instantaneous_rate', 'pdf')
+
+%% Binning Method
+clear 
+clc
+close all
+
+load 'lifoustim.mat'
+t_max = 30;
+rates = zeros(length(spikes), 30/dt);
+
+for i = 1:length(spikes)
+    [t, rates(i,:)] = binned_rate(spikes{i}, 0.05, dt, t_max);
+end
+
+
+f = figure();
+set(f, 'paperunits', 'centimeters')
+set(f, 'papersize', [10, 10])
+set(f, 'paperposition', [0 0 10 10])
+hold on
+plot(t, rates(1,:), 'displayname', 'trial 1')
+plot(t, mean(rates,1), 'displayname', 'average rate')
+
+xlabel('time [s]', 'fontsize', 10)
+ylabel('firing rate [Hz]','fontsize', 10)
+legend('show')
+box('off')
+
+saveas(gcf, 'binned_rate', 'pdf')
+
+
+%% Convolution Method
+clear 
+clc
+close all
+
+load 'lifoustim.mat'
+t_max = 30;
+rates = zeros(length(spikes), 30/dt);
+
+for i = 1:length(spikes)
+    [t, rates(i,:)] = convolution_rate(spikes{i}, 0.05, dt, t_max);
+end
+
+
+f = figure();
+set(f, 'paperunits', 'centimeters', 'papersize', [10, 10], 'paperposition', [0 0 10 10])
+hold on
+plot(t, rates(1,:), 'displayname', 'trial 1')
+plot(t, mean(rates,1), 'displayname', 'average rate')
+
+xlabel('time [s]', 'fontsize', 10)
+ylabel('firing rate [Hz]','fontsize', 10)
+legend('show')
+box('off')
+
+saveas(gcf, 'convolved_rate', 'pdf')
+