% 2019-01-27

clearvars;
close all;

KSE = 2;
            % 1: Single spike train
            % 2: Multiple spike trains
            
SPKTRNKSE = 3;
            % 1: Uniform spike train distribution
            % 2: Poisson spike train distribution
            % 3: Gaussian spike train distribution


% Time definitions

Tmax = 1000;
dt = 0.1;
t = 0:dt:Tmax;

            
% White (Gaussian) noise

eta = randn(1,length(t));

if KSE == 1

    % Spike-train generation

    if SPKTRNKSE == 1

        % Uniform distribution

        Freq = 30;              
        ISI = 1000/Freq;        
        tspk = dt:ISI:Tmax;
        Spk = zeros(1,length(t));
        jspk = floor(tspk/dt);
        Nspk = length(tspk);
        Spk(jspk) = 1;


     elseif SPKTRNKSE == 2


        % Poisson Distributed Spike Train:
        % Generation of a Poisson distributed input spike train with a
        % refractory period equal to the minimum of the ISIs from the uniform
        % distributions: ISI(1) (refractoriness: rejection of ISIin<ISI(1))
        % Rate: Freqin

        Freq = 30;   
        ISIm = 1;

        [tspk] = PoissionDistributedSpikeTrains(Freq,ISIm,t,Tmax);

        ISI = diff(tspk);
        jspk = floor(tspk/dt);
        Nspk = length(tspk);
        Spk = zeros(1,length(t));
        Spk(jspk) = 1;


    elseif SPKTRNKSE == 3

        % Gaussian distribution

        Freqmean = 30;     
        Freqvar = 10;
        Freqstd = sqrt(Freqvar);
        ISImean = 1000/Freqmean;      
        tspk = zeros(1);
        ISI = zeros(1);

        for j=2:floor(Tmax/dt)
            ISI(j-1) = Freqmean+Freqstd*randn;
            if tspk(j-1)+ISI(j-1) < Tmax
                tspk(j) = tspk(j-1)+ISI(j-1);
            else
                break
            end                
        end

        Spk = zeros(1,length(t));
        jspk = floor(tspk/dt);
        jspk(1) = 1;
        Nspk = length(tspk);
        Spk(jspk) = 1;

    end

    figure
    hold on
    plot(t,Spk,'-b','linewidth',2);
    axis([0 Tmax 0 1.2]);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');

    figure
    hold on
    hist(ISI,40);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');
    ylabel('ISI');
    title('ISI Histogram');

    % Coefficient of variation

    CV = std(ISI)/mean(ISI);

elseif KSE == 2
    
    Ntr = 40;
    
    if SPKTRNKSE == 1
        
        % Uniform distribution (with uniform noise)

        alpha = 4;
        beta = -2;

        Freq = 30;              
        ISIbaseline = 1000/Freq;  
        tspkbaseline = dt:ISIbaseline:Tmax;
        tspk = tspkbaseline'+alpha*rand(length(tspkbaseline),Ntr)-1;
        for k=1:Ntr
            if tspk(1,k)<0
                tspk(1,k)= dt;
                jspk(1,k) = 1;
            end
        end
        ISI = diff(tspk);   
        jspk = floor(tspk/dt);
        for k=1:Ntr
            if jspk(1,k) == 0
                jspk(1,k) = 1;
            end
        end
        Nspk = zeros(1,Ntr);
        for k=1:Ntr
            Nspk(k) = length(tspk(:,k));
        end
        
        
       
        
        ISIglobal = ISI(:,1);
        for k=2:Ntr
            ISIglobal = [ISIglobal; ISI(:,k)];
        end
            
    elseif SPKTRNKSE == 2


        % Poisson Distributed Spike Train:
        % Generation of a Poisson distributed input spike train with a
        % refractory period equal to the minimum of the ISIs from the uniform
        % distributions: ISI(1) (refractoriness: rejection of ISIin<ISI(1))
        % Rate: Freqin

        Freq = 40;   
        ISIm = 1;

        tspk = zeros(1,Ntr);
        Nspk = zeros(1,Ntr);
        for k=1:Ntr
            aux = PoissionDistributedSpikeTrains(Freq,ISIm,t,Tmax);
            Nspk(k) = length(aux);
            tspk(1:Nspk(k),k) = aux';
        end
        
        ISI = diff(tspk);
        jspk = floor(tspk/dt);
                      
        ISIglobal = ISI(1:Nspk(1)-1,1);
        for k=2:Ntr
            ISIglobal = [ISIglobal; ISI(1:Nspk(k)-1,k)];
        end
        
    elseif SPKTRNKSE == 3
        
        % Gaussian distribution

        Freqmean = 30;     
        Freqvar = 10;
        Freqstd = sqrt(Freqvar);
        ISImean = 1000/Freqmean;      
        tspk = zeros(1);
        ISI = zeros(1);
       
        

        for k=1:Ntr
            j = 1;
            ISI(j,k) = Freqmean+Freqstd*randn;            
            while sum(ISI(:,k))<Tmax
                j = j+1;
                ISI(j,k) = Freqmean+Freqstd*randn;                
            end
            Nspk(k) = j;
            tspk(1,k) = dt;
            for j=2:Nspk(k)
                tspk(j,k) = tspk(j-1,k)+ISI(j-1,k);
            end   
        end

        ISI = diff(tspk);
        jspk = floor(tspk/dt);
                      
        ISIglobal = ISI(1:Nspk(1)-1,1);
        for k=2:Ntr
            ISIglobal = [ISIglobal; ISI(1:Nspk(k)-1,k)];
        end
        
        
    end
    
    
    
    % Coefficient of variation
    
    CV = std(ISIglobal)/mean(ISIglobal);
    
    % Fano factor
    
    FanoF = var(Nspk)/mean(Nspk);
    
    
    % ISI histogram (all spike trains)
    
    figure
    hold on
    hist(ISIglobal,40);
    set(gca,'fontsize',24);
    xlabel('ISI');    
    title('ISI Histogram');
    
    
    % Spike-count rate histogram: number of spikes per train normalized
    % by the interval time length (Tmax)
    
    figure
    hold on
    hist(Nspk/Tmax*1000,40);
    axis([0 100 0 10]);
    set(gca,'fontsize',24);
    xlabel('Spike-count rate  [Hz]');    
    title('Histogram: # of spikes per train');
    
    % Raster plot
    
    figure
    hold on
    for k=1:Ntr
        for j=1:Nspk(k)
            plot([tspk(j,k),tspk(j,k)],[k-0.9 k-0.1],'-b','linewidth',2);
        end
    end
    axis([0 Tmax 0 Ntr]);
    set(gca,'fontsize',24);
    xlabel('t  [ms]');   
    ylabel('trial #');
    title('Raster Plot');
    
    % Computation of the time-dependent Spike-count firing rage
    
    Nbin = 30;
    Deltabin = Tmax/Nbin;
    Spkbin = zeros(max(Nspk),Ntr);
    
    for k=1:Ntr
        Spkbin(1:Nspk(k),k) = ceil(tspk(1:Nspk(k),k)/Deltabin);        
    end
    Spkbinvec = Spkbin(1:Nspk(1),1);
    for k=2:Ntr
        Spkbinvec = [Spkbinvec;Spkbin(1:Nspk(k),k)];
    end
    
    Nspkbin = zeros(1,Nbin);
    for j=1:length(Spkbinvec)
        Nspkbin(Spkbinvec(j)) =  Nspkbin(Spkbinvec(j))+1;
    end
    Nspkbin = Nspkbin/Ntr/Deltabin*1000;
   
    
    figure
    hold on
    plot(Nspkbin,'ob','linewidth',2);
    plot(Nspkbin,'-b','linewidth',1);
    set(gca,'fontsize',24);
    xlabel('bin #'); 
    ylabel('spike-count rate  [Hz]');
    
    % Smoothing (convolution) for a single spike train
    
    Spk = zeros(Ntr,length(t));
    for k=1:Ntr
        for j=1:Nspk(k)  
            Spk(k,jspk(j,k)) = 1;
        end
    end
    
    % half-size of kernel
    sz = 100; 
    % Gaussian kernel
    gauss_ker = normpdf(-sz:sz,0, sz/2 ); 
    % Square kernel
    sq_ker = ones(1,sz*2+1);
    
    k = 1;
    SpkGauss = conv(gauss_ker,Spk(k,:))/sum(gauss_ker); 
    SpkSq = conv(sq_ker,Spk(k,:))/sum(sq_ker);
    
    figure
    hold on
    plot(Spk(k,:)-1,'ob','linewidth',1);
    plot(SpkGauss(sz+1:end-sz),'-r','linewidth',2);
    plot(SpkSq(sz+1:end-sz),'-g','linewidth',2);
    set(gca,'fontsize',24);
    axis([0 length(t) -0.001 0.05]);
    
    % Autocorrelation
    
    k = 1;
    XC = xcorr(Spk(k,:),Spk(k,:));
    XC_unbiased = xcorr(Spk(k,:),Spk(k,:),'unbiased');   
    
    figure
    hold on
    plot(XC,'-b','linewidth',2);
    plot(XC_unbiased,'-r','linewidth',2);
    set(gca,'fontsize',24);
    title('Autocorrelation')
    
    % Cross-correlation
    
    k1 = 1;
    k2 = 2;
    CC = xcorr(Spk(k1,:),Spk(k2,:));
    CC_unbiased = xcorr(Spk(k1,:),Spk(k2,:),'unbiased');   
    
    figure
    hold on
    plot(CC,'-b','linewidth',2);
    plot(CC_unbiased,'-r','linewidth',2);
    set(gca,'fontsize',24);
    title('Cross-correlation');
    
    
end



