%% License Info
%A Model of Feedforward, Global, and Lateral Inhibition in the Locust Visual System.

%This model examines the architecture and function of inhibitory mechanisms in the 
%visual system of locusts, namely those involved in the processing of inputs to a
%key looming-sensitive neuron, the lobula giant movement detector (LGMD).

%Copyright (c) 2026, Erik Olson, Travis Wiens, and Jack Gray

%CITATION:
%When using the model code for scientific publications, cite the following work:
%Olson EGN, Wiens TK, Gray JR. A model of feedforward, global, and lateral inhibition 
%in the locust visual system predicts responses to looming stimuli. Biol Cybern. 2021 
%Jun;115(3):245-265. doi: 10.1007/s00422-021-00876-8. Epub 2021 May 16. PMID: 33997912.

%This program is free software: you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation, either version 3 of the License, or
%(at your option) any later version.

%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.

%You should have received a copy of the GNU General Public License
%along with this program.  If not, see <https://www.gnu.org/licenses/>.

%Contact: erik.olson@usask.ca

%% Code
function [f,hist,t_bins] = histoFiringRate(G,t,bin_width,std_dev)
%NOTE: THIS VERSION IS FOR USE WITH SPIKE PULSES, NOT SYNAPTIC OUTPUTS
%Processes a signal with corresponding time values to detect spikes into a 
%histogram.  This is then converted to a series of firing rate values and 
%associated time values.  Each spike is assumed to have a height of one.  
%Variables:
%G = the original spiketrain data
%t = the corresponding time values [s]
%bin_width = bin width used for histogram generation [s]
%std_dev = standard deviation of the gaussian window in samples - must be
%an integer!

N = size(G,1); %number of points in the input dataset
spikes = []; %time values of spikes are stored in this array

window = 5*std_dev;

bin_count = ceil((t(N)+bin_width)/bin_width);%number of bins for the histogram
t_bins = bin_width*linspace(0,bin_count-1,bin_count)';%time values at bin centers
bin_edges = bin_width*linspace(-0.5,bin_count-0.5,bin_count+1)';%bin edge values for
%matlab's histogram function

for i = 2:N %iterate through the dataset
   if G(i)>0
       spike_count = G(i);
       spikes = [spikes;t(i)*ones(spike_count,1)];
   end
end

%Create the histogram from spike timing data
hist = histcounts(spikes,bin_edges);

%Smooth the histogram
w = gausswin(window);%Create gaussian window for convolution

hist_filt = conv(hist,w,'same')';

weights = conv(ones(1,bin_count),w,'same')';%weightings to normalize the 
%histogram relative to the gaussian

hist_filt = hist_filt./weights;

%Calculate firing rates
f = hist_filt/bin_width;