# 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()
