function [y,x_axis] = LIF_firingrate(times,bin_size,time_seconds,dt,interval,varargin)

% LIF_FIRINGRATE instantaneous firing rate generation from spike count
%
%   LIF_FIRINGRATE(A,WINSIZE,T,DT,I) where A is array of event times (in seconds) for the spike train,
%   WINSIZE is the required window-size (in seconds), T is the 
%   simulation period in seconds, DT is the time-step of the simulation, and I is the number of time-steps between window
%   movements (that is, so 1 = dt, 2 = 2*dt etc)
%   Generates a firing rate time-series of the trains specified by A using a moving window average.
%
%   LIF_FIRINGRATE(...,W) specifies the window type to be used, 'alpha' or 'gaussian'. 
%   By default, a Gaussian window is used.
%
%   [Y,X] = LIF_FIRINGRATE(...) returns the Y time-series (in spikes-per-second) generated by the windowing, and X the 
%   corresponding x-axis
%
%   REFERENCE: Dayan, P & Abbott, L. F. (2001) Theoretical Neuroscience. Cambridge, MA: MIT Press [Gaussian + alpha]
%              Nawrot, M., Aersten, A. & Rotter, S. (1999) "Single-trial estimation of neuronal firing rates: From
%              single-neuron spike trains to population activity" Journal of Neuroscience Methods, 94, 81-92. [original source?]
%
%   Mark Humphries 22/4/04

% co-efficient for alpha window
a = 1/bin_size;

% variables and storage
x_axis = dt:dt*interval:time_seconds;     % create x-axis
y = zeros(1,length(x_axis));
counter = 0;

if nargin == 6
    win_type = varargin{1};
else
    win_type = 'gaussian';
end

if ~strcmp(win_type,'gaussian') & ~strcmp(win_type,'alpha')
    error('Incorrect window type specified')        
end

% check there's something to process
if length(times) < 2
    % nothing to process
    return
end

% flip array if necessary
[r1 c1] = size(times);
if r1 > c1 times = times'; end

% % matrix method
% x_mat = repmat(x_axis',1,length(times));
% t_mat = repmat(times,length(x_axis),1);
% diff = x_mat - t_mat;
% 
% % free up memory
% clear x_mat t_mat
% 
% switch win_type
% case 'gaussian'
%     window = 1/(bin_size*sqrt(2*pi)) .* exp(-(diff.^2)./(2*bin_size^2));
% case 'alpha'
%     window = a^2 .* diff .* exp(-a.*diff);
%     window(window < 0) = 0;            % rectify
% end
% y = sum(window');

% for-loop method - kept for comprehension (and for low memory overhead)
% pre-calculate terms for efficiency
g1 = 1/(bin_size*sqrt(2*pi));
g2 = (2*bin_size^2);
a1 = a^2;

for loop = dt:dt*interval:time_seconds
    counter = counter + 1;
    x = loop - times;
    switch win_type
    case 'gaussian'
        window = g1 * exp(-(x.^2)/g2); % note: Gaussian in Dayan & Abbott is slightly wrong (or unclear)
    case 'alpha'
        window = a1 .* x .* exp(-a.*x);
        window(window < 0) = 0;             % rectify
    end
    y(counter) = sum(window);
end