# input.py
# Generate cochlear input spikes to use auditory thalamocortical model
import numpy as np
rng = np.random.RandomState()
def cochlearInputSpikes(freqRange=[4800, 5200], #[125, 20000], #[9000, 11000],
numCenterFreqs=4, #100,
numCells=4*100, #10000, # should be ~100 * numCFs
duration=1000,
toneFreq=5000,
loudnessDBs=50,
plotRaster=False):
import matplotlib.pyplot as plt
import scipy.signal as dsp
import cochlea
fs = 100e3
# Make sound
t = np.arange(0, duration/1000.0, 1/fs)
s = dsp.chirp(t, toneFreq-100, t[-1], toneFreq+100)
s = cochlea.set_dbspl(s, loudnessDBs)
pad = np.zeros(int(10e-3 * fs))
sound = np.concatenate( (s, pad) )
# Run model
anf = cochlea.run_zilany2014(
sound,
fs,
anf_num=(numCells/numCenterFreqs, 0, 0), # the desired number of auditory nerve fibers per frequency channel (CF)
cf=(freqRange[0], freqRange[1], numCenterFreqs), # the center frequency(s) of the simulated auditory nerve fibers
seed=0,
powerlaw='approximate',
species='human',
)
# Accumulate spike trains
#anf_acc = th.accumulate(anf, keep=['cf', 'duration'])
#anf_acc.sort_values('cf', ascending=False, inplace=True)
# Plot auditory nerve response
if plotRaster:
pass
# generate list of spk times
spkTimes = [list(anf.iloc[i]['spikes']*1000.) for i in range(numCells)]
return spkTimes
# def cochlearInputSpikesBrianHears(freqRange=[9000,11000], # orig: [20, 20000],
# numCells=200, # orig: [3000]
# duration=1000,
# toneFreq=10000,
# plotRaster=False):
# from brian2 import *
# from brian2hears import *
# print(' Generating cochlear-like auditory input spikes using Brian Hears ...')
# cfmin, cfmax, cfN = freqRange[0]*Hz, freqRange[1]*Hz, numCells
# cf = erbspace(cfmin, cfmax, cfN)
# sound1 = tone(toneFreq*Hz, duration*ms)
# sound2 = whitenoise(duration*ms)
# sound = sound1+sound2
# sound = sound.ramp()
# gfb = Gammatone(sound, cf)
# ihc = FunctionFilterbank(gfb, lambda x: 3 * clip(x, 0, Inf)**(1.0 / 3.0))
# # Leaky integrate-and-fire model with noise and refractoriness
# eqs = '''
# dv/dt = (I-v)/(1*ms)+0.2*xi*(2/(1*ms))**.5 : 1 (unless refractory)
# I : 1
# '''
# G = FilterbankGroup(ihc, 'I', eqs, reset='v=0', threshold='v>1', refractory=5*ms, method='euler')
# # Run, and raster plot of the spikes
# M = SpikeMonitor(G)
# run(sound.duration)
# if plotRaster:
# plot(M.t / ms, M.i, '.')
# plt.show()
# # generate list of spk times
# spkts = list(M.t)
# spkids = list(M.i)
# spkTimes = [[] for i in range(numCells)]
# for spkt, spkid in zip(spkts, spkids):
# spkTimes[spkid].append(float(spkt)*1000.)
# return spkTimes
def poisson_generator(rate, t_start=0.0, t_stop=1000.0, seed=None):
"""
Returns a SpikeTrain whose spikes are a realization of a Poisson process
with the given rate (Hz) and stopping time t_stop (milliseconds).
Note: t_start is always 0.0, thus all realizations are as if
they spiked at t=0.0, though this spike is not included in the SpikeList.
Inputs:
-------
rate - the rate of the discharge (in Hz)
t_start - the beginning of the SpikeTrain (in ms)
t_stop - the end of the SpikeTrain (in ms)
array - if True, a np array of sorted spikes is returned,
rather than a SpikeTrain object.
Examples:
--------
>> gen.poisson_generator(50, 0, 1000)
>> gen.poisson_generator(20, 5000, 10000, array=True)
See also:
--------
inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
"""
rng = np.random.RandomState(seed)
#number = int((t_stop-t_start)/1000.0*2.0*rate)
# less wasteful than double length method above
n = (t_stop-t_start)/1000.0*rate
number = np.ceil(n+3*np.sqrt(n))
if number<100:
number = min(5+np.ceil(2*n),100)
if number > 0:
isi = rng.exponential(1.0/rate, int(number))*1000.0
if number > 1:
spikes = np.add.accumulate(isi)
else:
spikes = isi
else:
spikes = np.array([])
spikes+=t_start
i = np.searchsorted(spikes, t_stop)
extra_spikes = []
if i==len(spikes):
# ISI buf overrun
t_last = spikes[-1] + rng.exponential(1.0/rate, 1)[0]*1000.0
while (t_last<t_stop):
extra_spikes.append(t_last)
t_last += rng.exponential(1.0/rate, 1)[0]*1000.0
spikes = np.concatenate((spikes,extra_spikes))
else:
spikes = np.resize(spikes,(i,))
return spikes
def inh_poisson_generator(rate, t, t_stop, seed=None):
"""
Returns a SpikeTrain whose spikes are a realization of an inhomogeneous
poisson process (dynamic rate). The implementation uses the thinning
method, as presented in the references.
Inputs:
-------
rate - an array of the rates (Hz) where rate[i] is active on interval
[t[i],t[i+1]]
t - an array specifying the time bins (in milliseconds) at which to
specify the rate
t_stop - length of time to simulate process (in ms)
array - if True, a np array of sorted spikes is returned,
rather than a SpikeList object.
Note:
-----
t_start=t[0]
References:
-----------
Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier
Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories
Neural Comput. 2007 19: 2958-3010.
Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
Examples:
--------
>> time = arange(0,1000)
>> stgen.inh_poisson_generator(time,sin(time), 1000)
See also:
--------
poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
"""
rng = np.random.RandomState(seed)
if np.shape(t)!=np.shape(rate):
raise ValueError('shape mismatch: t,rate must be of the same shape')
# get max rate and generate poisson process to be thinned
rmax = np.max(rate)
ps = poisson_generator(rate=rmax, t_start=t[0], t_stop=t_stop, seed=seed)
# return empty if no spikes
if len(ps) == 0:
np.array([])
# gen uniform rand on 0,1 for each spike
rn = np.array(rng.uniform(0, 1, len(ps)))
# instantaneous rate for each spike
idx = np.searchsorted(t, ps) - 1
#spike_rate = rate[idx]
spike_rate = np.array([rate[i] for i in idx])
# thin and return spikes
spike_train = ps[rn<spike_rate/rmax]
return list(spike_train)
# main
#spk=cochlearInputSpikes(plotRaster=1)