% Computational Neuroscience II: Foundations of Neural Coding, SS 2003
% Solutions for Assignment 4 (15.05.03) 
% (c) written by Christian Leibold and Richard Kempter

% --- forget about everything in the past
clear;

path(path,"./octave_functions");% where we will find some 
% additional octave repository functions
% that do not come with the main package. Alternatively, you download them from
% http://octave.sourceforge.net/  (and install them properly on your system)


%%%%%%%% Excercise 9:

% --- define grid
DT=.05;
t=0:DT:20;


% ---- Define two sets of data and plot them

tau = 1; f = 0.25;
v1=t.*exp(-t/tau).*sin(2*pi*f*t); % sin modulated by an alpha function

figure(1)
title('Function 1: (tau,f) = (1,0.25)')
xlabel('x-axis')
ylabel('y-axis')
plot(t,v1);

%

tau = 1.5; f = 0.33;
v2=t.*exp(-t/tau).*sin(2*pi*f*t);

figure(2)
title('Function 2: (tau,f) = (1.5,0.33)')
xlabel('x-axis')
ylabel('y-axis')
plot(t,v2);

%----

%---- produce vector of time shifts 
n= length(v1);
s=-(n-1)*DT:DT:(n-1)*DT;


%--- invoke the function 'correlate', and plot the result
Q=correlate(v1,v2);

figure(3);
subplot(3,1,1);
title('Crosscorrelation -- correlate')
xlabel('Difference');
ylabel('Correlation');
plot(s,Q);


%-- now the easy way
Qx = xcorr(v1,v2);

subplot(3,1,2);
title('Crosscorrelation -- xcorr')
xlabel('Difference');
ylabel('Correlation');
plot(s,Qx);

% let's have a look at the noise vectors

load noise1.dat;
load noise2.dat;
maxlag = 100;
[Qnoise,lags] = xcorr(noise1,noise2,maxlag,'unbiased'); % try also the
                                                       % 'coeff' option

subplot(3,1,3);
title('Crosscorrelation of noise vectors -- xcorr')
xlabel('Difference');
ylabel('Correlation');
plot(lags,Qnoise);


%%%%% Exercise 10
%--- and now again, the spike trains

load 'spiketrain1.dat';
load 'spiketrain2.dat';

binsize = 5; % in milliseconds 

number_of_bins1 = max(spiketrain1)/binsize;
number_of_bins2 = max(spiketrain2)/binsize;

h1=hist(spiketrain1,number_of_bins1); 
h2=hist(spiketrain2,number_of_bins2);

[Q11,l11]=xcorr(h1,h1,500/binsize);
[Q12,l12]=xcorr(h1,h2,3000/binsize);
[Q22,l22]=xcorr(h2,h2,2000/binsize);

figure(4)
subplot(3,1,1)
title('s1%s1');
xlabel('time difference [ms]')
plot(binsize*l11, Q11);
subplot(3,1,2)
title('s1%s2');
xlabel('time difference [ms]')
plot(binsize*l12,Q12);
subplot(3,1,3)
title('s2%s2');
xlabel('time difference [ms]')
ylabel('correlation')
plot(binsize*l22, Q22);


% ------- For specialsts: Correlation with  fft
figure(5)

% f12=ifft((fft(v2,2*n-1)'.'.*fft(v1,2*n-1))); 
f12=ifft((conj(fft(v2,2*n-1)) .* fft(v1,2*n-1))); % eqivalent to the
                                                  % previous line! 

xlabel('x-difference')
title('Crosscorrelation -- fft')
plot(s, real(ifftshift(f12)));