| | 610 | |
| | 611 | def psth(self, events, time_bin=2, t_min=50, t_max=50, display = False, kwargs={}): |
| | 612 | """ |
| | 613 | Return the psth of the spike times contained in the SpikeTrain according to selected events, |
| | 614 | on a time window t_spikes - tmin, t_spikes + tmax |
| | 615 | |
| | 616 | Inputs: |
| | 617 | events - Can be a SpikeTrain object (and events will be the spikes) or just a list |
| | 618 | of times |
| | 619 | t_min - Time (>0) to average the signal before an event, in ms (default 0) |
| | 620 | t_max - Time (>0) to average the signal after an event, in ms (default 100) |
| | 621 | display - if True, a new figure is created. Could also be a subplot. |
| | 622 | kwargs - dictionary contening extra parameters that will be sent to the plot |
| | 623 | function |
| | 624 | |
| | 625 | Examples: |
| | 626 | >> spk.psth(spktrain, t_min = 50, t_max = 150) |
| | 627 | >> spk.psth(spktrain, ) |
| | 628 | >> spk.psth(range(0,1000,10), display=True) |
| | 629 | |
| | 630 | See also |
| | 631 | SpikeTrain.spike_histogram |
| | 632 | """ |
| | 633 | |
| | 634 | if isinstance(events, SpikeTrain): |
| | 635 | events = events.spike_times |
| | 636 | assert (t_min >= 0) and (t_max >= 0), "t_min and t_max should be greater than 0" |
| | 637 | assert len(events) > 0, "events should not be empty and should contained at least one element" |
| | 638 | |
| | 639 | spk_hist = self.time_histogram(time_bin) |
| | 640 | subplot = get_display(display) |
| | 641 | count = 0 |
| | 642 | t_min_l = numpy.floor(t_min/time_bin) |
| | 643 | t_max_l = numpy.floor(t_max/time_bin) |
| | 644 | result = numpy.zeros((t_min_l+t_max_l), numpy.float32) |
| | 645 | t_start = numpy.floor(self.t_start/time_bin) |
| | 646 | t_stop = numpy.floor(self.t_stop/time_bin) |
| | 647 | for ev in events: |
| | 648 | ev = numpy.floor(ev/time_bin) |
| | 649 | if ((ev - t_min_l )> t_start) and (ev + t_max_l ) < t_stop: |
| | 650 | count += 1 |
| | 651 | result += spk_hist[(ev-t_min_l):ev+t_max_l] |
| | 652 | result /= count |
| | 653 | |
| | 654 | if not subplot or not HAVE_PYLAB: |
| | 655 | return result |
| | 656 | else: |
| | 657 | xlabel = "Time (ms)" |
| | 658 | ylabel = "PSTH" |
| | 659 | time = numpy.linspace(-t_min, t_max, (t_min+t_max)/time_bin) |
| | 660 | set_labels(subplot, xlabel, ylabel) |
| | 661 | subplot.plot(time, result, c='k', **kwargs) |
| | 662 | xmin, xmax, ymin, ymax = subplot.axis() |
| | 663 | subplot.plot([0,0],[ymin, ymax], c='r') |
| | 664 | set_axis_limits(subplot, -t_min, t_max, ymin, ymax) |
| | 665 | pylab.draw() |
| | 666 | return result |
