Cortical model with reinforcement learning drives realistic virtual arm (Dura-Bernal et al 2015)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:183014
We developed a 3-layer sensorimotor cortical network of consisting of 704 spiking model-neurons, including excitatory, fast-spiking and low-threshold spiking interneurons. Neurons were interconnected with AMPA/NMDA, and GABAA synapses. We trained our model using spike-timing-dependent reinforcement learning to control a virtual musculoskeletal human arm, with realistic anatomical and biomechanical properties, to reach a target. Virtual arm position was used to simultaneously control a robot arm via a network interface.
References:
1 . Dura-Bernal S, Zhou X, Neymotin SA, Przekwas A, Francis JT, Lytton WW (2015) Cortical Spiking Network Interfaced with Virtual Musculoskeletal Arm and Robotic Arm. Front Neurorobot 9:13 [PubMed]
2 . Dura-Bernal S, Li K, Neymotin SA, Francis JT, Principe JC, Lytton WW (2016) Restoring behavior via inverse neurocontroller in a lesioned cortical spiking model driving a virtual arm. Front. Neurosci. Neuroprosthetics 10:28
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Neocortex M1 pyramidal pyramidal tract L5B cell; Neocortex M1 pyramidal intratelencephalic L2-5 cell; Neocortex M1 interneuron basket PV cell; Neocortex fast spiking (FS) interneuron; Neostriatum fast spiking interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; Python (web link to model);
Model Concept(s): Synaptic Plasticity; Learning; Reinforcement Learning; STDP; Reward-modulated STDP; Sensory processing; Motor control;
Implementer(s): Neymotin, Sam [samn at neurosim.downstate.edu]; Dura, Salvador [ salvadordura at gmail.com];
Search NeuronDB for information about:  Neocortex M1 pyramidal intratelencephalic L2-5 cell; Neocortex M1 pyramidal pyramidal tract L5B cell; Neocortex M1 interneuron basket PV cell; GabaA; AMPA; NMDA; Gaba; Glutamate;
/
arm2dms_modeldb
mod
msarm
stimdata
README.html
analyse_funcs.py
analysis.py
armGraphs.py
arminterface_pipe.py
basestdp.hoc
bicolormap.py
boxes.hoc *
bpf.h *
col.hoc
colors.hoc *
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
filtutils.hoc *
grvec.hoc
hinton.hoc *
hocinterface.py
infot.hoc *
init.hoc
intfsw.hoc *
labels.hoc
load.hoc
load.py
local.hoc *
main.hoc
main_demo.hoc
main_neurostim.hoc
misc.h *
misc.py *
msarm.hoc
network.hoc
neuroplot.py *
neurostim.hoc
nload.hoc
nqs.hoc *
nqsnet.hoc *
nrnoc.hoc
params.hoc
perturb.hoc
python.hoc
pywrap.hoc *
run.hoc
runbatch_neurostim.py
runsim_neurostim
samutils.hoc *
saveoutput.hoc
saveoutput2.hoc
setup.hoc *
sim.hoc
sim.py
sim_demo.py
simctrl.hoc *
stats.hoc *
stim.hoc
syncode.hoc *
units.hoc *
vector.py
xgetargs.hoc *
                            
# neuroplot.py -- Python script with routines for making pylab plots of neural
#   data.

# Script with pylab routines for plotting and processing neural data.
# 
# Last update: 10/31/11 (georgec)

# Import pylab so the functions and classes don't have to.
import pylab

#
# Data file I/O functions
#

# Load numpy arrays from a file. 
def ldnarrs(fstem):
   fname = 'savedpydata/%s.txt' % fstem
   return pylab.loadtxt(fname)

# Save numpy arrays to a file.
def svnarrs(fstem, X):
   fname = 'savedpydata/%s.txt' % fstem
   pylab.savetxt(fname, X)

#
# Data plotting functions
#

# Plot numpy arrays from a file.
# Note that tvec needs to be the first saved vector.  The other vectors
# will plot as separate traces. 
def plotnarrs(fstem, numtraces=1):
   vecsstr = 'tvec,vec'
   for ii in range(numtraces-1):
      vecsstr += ',vec%d' % (ii+2)
   exec('%s = ldnarrs("%s")' % (vecsstr, fstem))
   
   # Plot first plot.
   pylab.plot(tvec, vec)
   
   # Plot the remaining plots.
   for ii in range(numtraces-1):
      exec('pylab.plot(tvec,vec%d)' % (ii+2))

def plot_spike_times(tvec):
   xmarg = 10   # x display margin
   ymarg = 0.1  # y display margin
   if (tvec == []):
      return
   avec = pylab.ones(len(tvec))
   pylab.stem(tvec, avec)
#   pylab.xlim(-xmarg, sim_duration + xmarg)
   pylab.ylim(-ymarg, 1 + ymarg)
   pylab.title('Spike Times')
   pylab.xlabel('Time')

def plot_spike_raster(tvec, vec):
   xmarg = 10  # x display margin
   ymarg = 10  # y display margin
   marks = pylab.scatter(tvec, vec, marker='+')
   pylab.setp(marks, color='b')
#   pylab.xlim(-xmarg + simview_min_t, simview_max_t + xmarg)
   pylab.ylim(-ymarg, max(vec) + ymarg)
   pylab.title('Spike Events')
   pylab.xlabel('Time')
   pylab.ylabel('Unit ID #')

def plot_power_spect(x_fvec,x_pvec,logx=True,logy=False,\
   putupbands=True,putuplabels=True):
   freqs = pylab.copy(x_fvec)
   if ((not logx) and (not logy)):
      pylab.plot(freqs,x_pvec)
      if putuplabels:
         pylab.ylabel('Power')
   elif (logx and (not logy)):
      if (freqs[0] == 0.0):
         freqs[0] = 0.1
      pylab.semilogx(freqs,x_pvec)
      if putuplabels:
         pylab.ylabel('Power')
   elif ((not logx) and logy):
      pylab.semilogy(freqs,x_pvec)
      if putuplabels:
         pylab.ylabel('Log Power')
   else:
      if (freqs[0] == 0.0):
         freqs[0] = 0.1
      pylab.loglog(freqs,x_pvec)
      if putuplabels:
         pylab.ylabel('Log Power')
   if putuplabels:
      pylab.title('Spectrum')
   pylab.xlim(0,freqs[-1])
   if putuplabels:
      pylab.xlabel('Frequency (Hz)')
   if (putupbands):
      put_up_band_bounds(logx,logy)

def put_up_band_bounds(logx=True,logy=False):
   line = pylab.axvline(1, color='k')
   ax = line.get_axes()
   ymax = ax.get_ybound()[-1]
   pylab.axvline(4, color='k')
   pylab.axvline(7.5, color='k')
   pylab.axvline(14, color='k')
   pylab.axvline(30, color='k')
   pylab.axvline(70, color='k')
   if logx:
      if (not logy):
         pylab.text(1.75,ymax * 0.75,'delta')
         pylab.text(4.5,ymax * 0.75,'theta')
         pylab.text(8,ymax * 0.75,'alpha')
         pylab.text(17,ymax * 0.75,'beta')
         pylab.text(33,ymax * 0.75,'gamma')
      else:
         pylab.text(1.75,ymax * 0.10,'delta')
         pylab.text(4.5,ymax * 0.05,'theta')
         pylab.text(8,ymax * 0.10,'alpha')
         pylab.text(17,ymax * 0.05,'beta')
         pylab.text(33,ymax * 0.10,'gamma')
   else:
      if (not logy):
         pylab.text(1,ymax * 0.75,'delta')
         pylab.text(2,ymax * 0.7,'theta')
         pylab.text(8,ymax * 0.75,'alpha')
         pylab.text(20,ymax * 0.7,'beta')
         pylab.text(45,ymax * 0.75,'gamma')
      else:
         pylab.text(1,ymax * 0.10,'delta')
         pylab.text(2,ymax * 0.05,'theta')
         pylab.text(8,ymax * 0.10,'alpha')
         pylab.text(20,ymax * 0.05,'beta')
         pylab.text(45,ymax * 0.10,'gamma')

def plot_specgram(freqs, ts, Pxx, logz=False):
   if logz:
      pylab.imshow(pylab.log10(Pxx), origin='lower', aspect='auto', \
         extent=(0,ts[-1]*1000.0+self.ts[0]*1000.0,0,freqs[-1]))
      pylab.title('Log Spectrogram')
   else:
      pylab.imshow(Pxx, origin='lower', aspect='auto', \
         extent=(0,ts[-1]*1000.0+self.ts[0]*1000.0,0,freqs[-1]))
      pylab.title('Spectrogram')
   pylab.xlabel('Time (ms)')
   pylab.ylabel('Frequency (Hz)')
   pylab.spectral()
   pylab.colorbar()

def plot_band_specgram(bandnames, bts, Ptxx):
   pylab.imshow(Ptxx, origin='lower', aspect='auto', \
      extent=(0,bts[-1]*1000.0+self.bts[0]*1000.0,0,len(bandnames)), \
      interpolation='nearest')
   pylab.title('Band Spectrogram')
   pylab.xlabel('Time (ms)')
   pylab.ylabel('Spectral Band')
   pylab.yticks(pylab.arange(len(bandnames))+0.5, bandnames)
   pylab.spectral()
   pylab.colorbar()

def plot_specgram_spectdist(freqs,ts,Pxx,show='quartiles',logx=True,logy=False,putupbands=True):
   import scipy.stats

   # Get the mean, median, range, and standard deviation scores.
   mean_vec = Pxx.mean(axis=1)      # average over time
   std_vec = Pxx.std(axis=1)        # stdev over time
   upperstd_vec = mean_vec + std_vec
   lowerstd_vec = mean_vec - std_vec
   upperrng_vec = Pxx.max(axis=1)
   lowerrng_vec = Pxx.min(axis=1)

   # Calculate the 25th, 50th (median) and 75th percentile scores.
   p25_vec = pylab.zeros(len(mean_vec))
   median_vec = pylab.zeros(len(mean_vec))
   p75_vec = pylab.zeros(len(mean_vec))
   for ii in range(len(freqs)):
      p25_vec[ii] = scipy.stats.scoreatpercentile(Pxx[ii,:], 25)
      median_vec[ii] = scipy.stats.scoreatpercentile(Pxx[ii,:], 50)
      p75_vec[ii] = scipy.stats.scoreatpercentile(Pxx[ii,:], 75)

   if ((not logx) and (not logy)):
      if (show in ['meanstd','meanmed']):
         pylab.plot(freqs,mean_vec)
      if (show in ['meanmed','quartiles']):
         pylab.plot(freqs,median_vec)
      if (show in ['meanstd']):
         pylab.plot(freqs,upperstd_vec,'g')
         pylab.plot(freqs,lowerstd_vec,'g')
      if (show in ['quartiles']):
         pylab.plot(freqs,p75_vec,'g')
         pylab.plot(freqs,p25_vec,'g') 
      if (show in ['meanstd','meanmed','quartiles']):
         pylab.plot(freqs,upperrng_vec,'r')
         pylab.plot(freqs,lowerrng_vec,'r')
      pylab.ylabel('Power')
   elif (logx and (not logy)):
      if (freqs[0] == 0.0):
         freqs[0] = 0.1
      if (show in ['meanstd','meanmed']):
         pylab.semilogx(freqs,mean_vec)
      if (show in ['meanmed','quartiles']):
         pylab.semilogx(freqs,median_vec)
      if (show in ['meanstd']):
         pylab.semilogx(freqs,upperstd_vec,'g')
         pylab.semilogx(freqs,lowerstd_vec,'g')
      if (show in ['quartiles']):
         pylab.semilogx(freqs,p75_vec,'g')
         pylab.semilogx(freqs,p25_vec,'g') 
      if (show in ['meanstd','meanmed','quartiles']):
         pylab.semilogx(freqs,upperrng_vec,'r')
         pylab.semilogx(freqs,lowerrng_vec,'r')
      pylab.ylabel('Power')
   elif ((not logx) and logy):
      if (show in ['meanstd','meanmed']):
         pylab.semilogy(freqs,mean_vec)
      if (show in ['meanmed','quartiles']):
         pylab.semilogy(freqs,median_vec)
      if (show in ['meanstd']):
         pylab.semilogy(freqs,upperstd_vec,'g')
         pylab.semilogy(freqs,lowerstd_vec,'g')
      if (show in ['quartiles']):
         pylab.semilogy(freqs,p75_vec,'g')
         pylab.semilogy(freqs,p25_vec,'g') 
      if (show in ['meanstd','meanmed','quartiles']):
         pylab.semilogy(freqs,upperrng_vec,'r')
         pylab.semilogy(freqs,lowerrng_vec,'r')
      pylab.ylabel('Log Power')
   else:
      if (freqs[0] == 0.0):
         freqs[0] = 0.1
      if (show in ['meanstd','meanmed']):
         pylab.loglog(freqs,mean_vec)
      if (show in ['meanmed','quartiles']):
         pylab.loglog(freqs,median_vec)
      if (show in ['meanstd']):
         pylab.loglog(freqs,upperstd_vec,'g')
         pylab.loglog(freqs,lowerstd_vec,'g')
      if (show in ['quartiles']):
         pylab.loglog(freqs,p75_vec,'g')
         pylab.loglog(freqs,p25_vec,'g') 
      if (show in ['meanstd','meanmed','quartiles']):
         pylab.loglog(freqs,upperrng_vec,'r')
         pylab.loglog(freqs,lowerrng_vec,'r')
      pylab.ylabel('Log Power')
   pylab.title('Spectrogram Spectrum Distribution')
   pylab.xlim(0,freqs[-1])
   pylab.xlabel('Frequency (Hz)')
   if (putupbands):
      put_up_band_bounds(logx,logy)

#
# Data processing / analysis functions
#

def get_vec_subset(tvec, vec, mintvec=-1, maxtvec=-1, minvec=-1, maxvec=-1):
   if (mintvec == -1):
      mintvec = tvec.min()
   if (maxtvec == -1):
      maxtvec = tvec.max() + 0.0005  # allow inclusion of the max
   if (minvec == -1):
      minvec = vec.min()
   if (maxvec == -1):
      maxvec = vec.max() + 0.0005    # allow inclusion of the max
   filt_tvec = []
   filt_vec = []
   for ii in range(len(tvec)):
      # If the tvec and vec values are in the window, then include them.  Note
      # that the min boundaries include the minimum value, but the max bounds
      # do not include the max value, only values up to the max.
      if ((tvec[ii] >= mintvec) and (tvec[ii] < maxtvec) and \
         (vec[ii] >= minvec) and (vec[ii] < maxvec)):
         filt_tvec.append(tvec[ii])
         filt_vec.append(vec[ii])
   filt_tvec = pylab.array(filt_tvec)
   filt_vec = pylab.array(filt_vec)
   return filt_tvec, filt_vec

# Downsample a time sequence.
#   x_tvec - time stamps for sequence (units assumed to be in ms)
#   x_vec - amplitudes for sequence
#   newfs - new sample frequency (in Hz, must evenly divide the old fs)  
def downsample (x_tvec, x_vec, newfs):
   fsratio = (1000.0 / (x_tvec[1] - x_tvec[0])) / float(newfs)
   
   # If the fs ratio is invalid, give an error.
   if ((float(int(fsratio)) != fsratio) or (fsratio < 1)):
      print 'Error: The ratio of oldfs and newfs must be a positive integer'
      return None
      
   # Integerize the ratio.
   fsratio = int(fsratio)
   
   # Do the downsampling by averaging x_vec blockwise by blocks of size fsratio.
   npts = len(x_vec) / fsratio
   tvec = pylab.zeros(npts)
   vec = pylab.zeros(npts)
   for ii in range(npts):
      tvec[ii] = x_tvec[ii * fsratio]
      vec[ii] = x_vec[ii * fsratio:((ii+1) * fsratio)].mean()
      
   return tvec, vec

# Use convolution with a Hanning window to smooth a vector.
def smooth(x, win_size=11):
   y = pylab.zeros(x.size)
   smwin = pylab.hanning(win_size)
   smwin /= sum(smwin)  # scale the window so the total area is 1.0
   filtout = pylab.convolve(x,smwin,mode='valid')
   y[0:filtout.size] = filtout
   return y

def get_ave_fire_freqs(spks_tvec,spks_vec,num_units,fire_dur=0.0):
   if (fire_dur == 0.0):
      fire_dur = spks_tvec.max() - spks_tvec.min()
   fire_counts = pylab.zeros(num_units)
   for ii in range(len(spks_tvec)):
      fire_counts[spks_vec[ii]] += 1
   return fire_counts * 1000.0 / fire_dur

def get_unit_fire_bins(spks_tvec,spks_vec,num_units,fire_dur=0.0,bin_dur=100.0):
   if (fire_dur == 0.0):
      fire_dur = spks_tvec.max() - spks_tvec.min()
   num_bins = int(fire_dur / bin_dur)
   fire_bin_counts = pylab.zeros((num_units,num_bins))
   for ii in range(num_units):
      tvec, vec = get_vec_subset(spks_tvec, spks_vec, minvec=ii, maxvec=ii+1)
      n, bins = pylab.histogram(tvec,bins=num_bins,range=(0.0,fire_dur))
      fire_bin_counts[ii,:] = n[:]
   return bins,fire_bin_counts

# Given full spectrogram information, return a matrix where columns are all of 
# the time slices and rows are frequency band and the values are the sum, within 
# the time slice, of all of the power in that band.
def get_band_waterfall(freqs,ts,Pxx,normlz=False):
   bandnames = ['<Delta','Delta','Theta','Alpha','Beta','Gamma','>Gamma']
   bts = pylab.copy(ts)
   Ptxx = pylab.zeros((len(bandnames),len(ts)))
   for tt in range(len(ts)):
      for ff in range(len(freqs)):
         if (freqs[ff] < 1.0):
            Ptxx[0,tt] += Pxx[ff,tt]  # sub-delta
         elif ((freqs[ff] >= 1.0) and (freqs[ff] < 4.0)):
            Ptxx[1,tt] += Pxx[ff,tt]  # delta
         elif ((freqs[ff] >= 4.0) and (freqs[ff] < 7.5)):
            Ptxx[2,tt] += Pxx[ff,tt]  # theta
         elif ((freqs[ff] >= 7.5) and (freqs[ff] < 14.0)):
            Ptxx[3,tt] += Pxx[ff,tt]  # alpha
         elif ((freqs[ff] >= 14.0) and (freqs[ff] < 30.0)):
            Ptxx[4,tt] += Pxx[ff,tt]  # beta
         elif ((freqs[ff] >= 30.0) and (freqs[ff] < 70.0)):
            Ptxx[5,tt] += Pxx[ff,tt]  # gamma
         else:
            Ptxx[6,tt] += Pxx[ff,tt]  # super-gamma
      if (normlz):
         Ptxx[:,tt] /= sum(Ptxx[:,tt])
   return bandnames,bts,Ptxx

#
# Data plotting / browsing classes
#

class SpecgramBrowser:
   """
   Spectrogram browser for looking at spectra time-slices in a waterfall plot.

   Methods:
     constructor(xxx):

   Attributes:
     count (static): number of instances

   Usage:
   >>> sb = SpecgramBrowser(freqs,ts,Pxx)
   """

   count = 0

   def __init__(self, freqs, ts, Pxx, Pxx2=None, logz=False, logx=True, logy=False, \
      putupbands=True, nointerp=False):
      # Set up browser parameters.
      self.freqs = freqs
      self.ts = ts
      self.Pxx = Pxx
      self.Pxx2 = Pxx2
      self.logz = logz
      self.logx = logx
      self.logy = logy
      self.ymaxlock = False
      self.putupbands = putupbands
      self.nointerp = nointerp
      self.specgram_fig = pylab.figure()  # handle to specgram_browser spectrogram figure
      if (self.Pxx2 == None):
         self.specgram_fig2 = 0   # handle to specgram_browser spectrogram 2 figure
      else:
         self.specgram_fig2 = pylab.figure()  # handle to specgram_browser spectrogram figure
      self.spect_fig = pylab.figure()     # handle to specgram_browser spectrum figure
      self.tindex = 0           # index to the time slice to be shown in the spectrum figure/s
      self.specgram_mark = 0    # handle to spectrogram time slice marker
      self.specgram_mark2 = 0   # handle to spectrogram 2 time slice marker

      # Increment the instance count.
      SpecgramBrowser.count += 1 

      # Draw the initial spectrogram/s.
      self._draw_specgram()

      # Draw the initial spectrum/a.
      self._draw_spect()

      # Attach an onclick event to grab mouse clicks for spectrogram window.
      cid = self.specgram_fig.canvas.mpl_connect('button_press_event', self._onclick)

      # Attach an onkey event to grab key presses for spectrogram window.
      cid2 = self.specgram_fig.canvas.mpl_connect('key_press_event', self._onkey)

      if (self.Pxx2 != None):
         # Attach an onclick event to grab mouse clicks for spectrogram window.
         cid3 = self.specgram_fig2.canvas.mpl_connect('button_press_event', self._onclick)

         # Attach an onkey event to grab key presses for spectrogram window.
         cid4 = self.specgram_fig2.canvas.mpl_connect('key_press_event', self._onkey)

      # Attach an onkey event to grab key presses for spectrum window.
      cid5 = self.spect_fig.canvas.mpl_connect('key_press_event', self._onkey2)

   def _draw_specgram(self):
      def figdrawer(self, Pxx):
         pylab.clf()
         if self.logz:
            if self.nointerp:
               pylab.imshow(pylab.log10(Pxx), origin='lower', aspect='auto', \
                  extent=(0,self.ts[-1]*1000.0+self.ts[0]*1000.0,0,self.freqs[-1]), \
                  interpolation='nearest')
            else:
               pylab.imshow(pylab.log10(Pxx), origin='lower', aspect='auto', \
                  extent=(0,self.ts[-1]*1000.0+self.ts[0]*1000.0,0,self.freqs[-1]))
            pylab.title('Log Spectrogram')
         else:
            if self.nointerp:
               pylab.imshow(Pxx, origin='lower', aspect='auto', \
                  extent=(0,self.ts[-1]*1000.0+self.ts[0]*1000.0,0,self.freqs[-1]), \
                  interpolation='nearest')
            else:
               pylab.imshow(Pxx, origin='lower', aspect='auto', \
                  extent=(0,self.ts[-1]*1000.0+self.ts[0]*1000.0,0,self.freqs[-1]))
            pylab.title('Spectrogram')
         pylab.xlabel('Time (ms)')
         pylab.ylabel('Frequency (Hz)')
         pylab.spectral()
         pylab.colorbar()

      pylab.figure(self.specgram_fig.number)
      figdrawer(self, self.Pxx)
      self.specgram_mark = pylab.axvline(self.ts[self.tindex] * 1000.0,color='r')
      if (self.Pxx2 != None):
         pylab.figure(self.specgram_fig2.number)
         figdrawer(self, self.Pxx2)
         self.specgram_mark2 = pylab.axvline(self.ts[self.tindex] * 1000.0,color='r')

   def _move_specgram_mark(self):
      # Reset the line marker according to the new time slice index.  
      self.specgram_mark.set_xdata(pylab.array([self.ts[self.tindex] * 1000.0,self.ts[self.tindex] * 1000.0]))

      # Redraw the spectrogram canvas.
      self.specgram_fig.canvas.draw()

      # If we have another spectrogram...
      if (self.Pxx2 != None):
         # Reset the line marker according to the new time slice index.  
         self.specgram_mark2.set_xdata(pylab.array([self.ts[self.tindex] * \
            1000.0,self.ts[self.tindex] * 1000.0]))

         # Redraw the spectrogram canvas.
         self.specgram_fig2.canvas.draw()

   def _draw_spect(self):
      pylab.figure(self.spect_fig.number)
      pylab.clf()
      freqs2 = pylab.copy(self.freqs)
      x_pvec = pylab.copy(self.Pxx[:,self.tindex])
      if (self.Pxx2 != None):
         x_pvec2 = pylab.copy(self.Pxx2[:,self.tindex])
      if ((not self.logx) and (not self.logy)):
         pylab.plot(freqs2,x_pvec)
         if (self.Pxx2 != None):
            pylab.plot(freqs2,x_pvec2,'r')
         pylab.ylabel('Power')
      elif (self.logx and (not self.logy)):
         if (freqs2[0] == 0.0):
            freqs2[0] = 0.1
         pylab.semilogx(freqs2,x_pvec)
         if (self.Pxx2 != None):
            pylab.semilogx(freqs2,x_pvec2,'r')
         pylab.ylabel('Power')
      elif ((not self.logx) and self.logy):
         pylab.semilogy(freqs2,x_pvec)
         if (self.Pxx2 != None):
            pylab.semilogy(freqs2,x_pvec2,'r')
         pylab.ylabel('Log Power')
      else:
         if (freqs2[0] == 0.0):
            freqs2[0] = 0.1
         pylab.loglog(freqs2,x_pvec)
         if (self.Pxx2 != None):
            pylab.loglog(freqs2,x_pvec2,'r')
         pylab.ylabel('Log Power')
      if (self.Pxx2 == None):
         pylab.title('Spectrum at %.2f s' % self.ts[self.tindex])
      else:
         pylab.title('Spectra at %.2f s' % self.ts[self.tindex])
      pylab.xlim(0,freqs2[-1])
      if (self.ymaxlock):
         maxlock = self.Pxx.max()
         if (self.Pxx2 != None):
            maxlock = max(maxlock, self.Pxx2.max())
         maxlock *= 1.0
         pylab.ylim(0,maxlock)
      pylab.xlabel('Frequency (Hz)')
      if (self.putupbands):
         put_up_band_bounds(self.logx,self.logy)

   def _onclick(self,event):
      print 'Selected...'
#      print 'button=%d, x=%d, y=%d, xdata=%f, ydata=%f' % \
#         (event.button, event.x, event.y, event.xdata, event.ydata)

      # Create freq bin upper boundaries and array of indices where value 
      # falls in bin.
      fbinrad = self.freqs[1] / 2.0
      fbinubs = self.freqs + fbinrad
      fis = pylab.find(fbinubs >= event.ydata)

      # Create time bin upper boundaries and array of indices where value
      # falls in bin.
      tbinrad = (self.ts[1] - self.ts[0]) / 2.0
      tbinubs = (self.ts + tbinrad) * 1000.0
      tis = pylab.find(tbinubs >= event.xdata)

      # Show the bin and value.
      print 'freq bin=%f, time bin=%f, value=%f' % (self.freqs[fis[0]], self.ts[tis[0]], self.Pxx[fis[0]][tis[0]])

      # Switch the spectrum view to the corresponding time slice.
      self.tindex = tis[0]
      self._draw_spect()

      # Redraw the spectrogram (moving the marker).
      self._draw_specgram()

   def _onkey(self,event):
      # Note: for reasons I don't understand, Tkinter callbacks trap the 
      # following keystrokes, so I cannot use them:
      #   's'
      #   'l'   
#      print 'you pressed', event.key, event.xdata, event.ydata
      if (event.key == 'left'):
         if (self.tindex > 0):
            self.tindex -= 1
         self._draw_spect()
         self._draw_specgram()
      elif (event.key == 'right'):
         if (self.tindex < len(self.ts)-1):
            self.tindex += 1
         self._draw_spect()
         self._move_specgram_mark()
      elif (event.key == 'i'):
         self.nointerp = not self.nointerp
         self._draw_specgram()

   def _onkey2(self,event):
      # Note: for reasons I don't understand, Tkinter callbacks trap the 
      # following keystrokes, so I cannot use them:
      #   's'
      #   'l'  
#      print 'you pressed', event.key, event.xdata, event.ydata
      if (event.key == 'left'):
         if (self.tindex > 0):
            self.tindex -= 1
         self._draw_spect()
         self._move_specgram_mark()
      elif (event.key == 'right'):
         if (self.tindex < len(self.ts)-1):
            self.tindex += 1
         self._draw_spect()
         self._move_specgram_mark()
      elif (event.key == 'x'):
         self.logx = not self.logx
         self._draw_spect()
      elif (event.key == 'y'):
         self.logy = not self.logy
         self._draw_spect()
      elif (event.key == 'b'):
         self.putupbands = not self.putupbands
         self._draw_spect()
      elif (event.key == 'm'):
         self.ymaxlock = not self.ymaxlock
         self._draw_spect()

Loading data, please wait...