root/trunk/src/analysis.py

#!/usr/bin/env python
# -*- coding: utf8 -*-
"""
analysis.py:

$Id$
"""

import os, numpy


def record(output, cfilename = 'SpikeTrainPlay.wav', fs=44100, enc = 'pcm26'):
    """ record the 'sound' produced by a neuron. Takes a spike train as the
    output.

    >>> record(my_spike_train)

    """


    # from the spike list
    simtime_seconds = (output.t_stop - output.t_start)/1000.
    #time = numpy.linspace(0, simtime_seconds , fs*simtime_seconds)
    (trace,time) = numpy.histogram(output.spike_times*1000., fs*simtime_seconds)


    # TODO convolve with proper spike...
    spike = numpy.ones((fs/1000.,)) # one ms

    trace = numpy.convolve(trace, spike, mode='same')#/2.0
    trace /= numpy.abs(trace).max() * 1.1

    from scikits.audiolab import wavwrite
    wavwrite(trace, cfilename, fs = fs, enc = enc)

def play(output):
    """
    plays a spike list to the audio output

    play(spike_list) where spike_list is a spike_list object

    see playing_with_simple_single_neuron.py for a sample use

    >>> play(my_spike_train)

    TODO: make it possible to play multiple spike trains in stereo
    """


    from tempfile import mkstemp
    fd, cfilename   = mkstemp('SpikeListPlay.wav')
    try:
        record(output, cfilename)
        import pyaudio
        import wave

        chunk = 1024
        wf = wave.open(cfilename, 'rb')
        p = pyaudio.PyAudio()

        # open stream
        stream = p.open(format =
                        p.get_format_from_width(wf.getsampwidth()),
                        channels = wf.getnchannels(),
                        rate = wf.getframerate(),
                        output = True)

        # read data
        data = wf.readframes(chunk)

        # play stream
        while data != '':
            stream.write(data)
            data = wf.readframes(chunk)

        stream.close()
        p.terminate()
    except:
        print "Error playing the SpikeTrain "
        # finally
        os.remove(cfilename)

    # Python 2.4 compatibility
    # finally:
    os.remove(cfilename)

def _dict_max(D):
    """
    For a dict containing numerical values, contain the key for the
    highest value. If there is more than one item with the same highest
    value, return one of them (arbitrary - depends on the order produced
    by the iterator).
    """
    max_val = max(D.values())
    for k in D:
        if D[k] == max_val:
            return k
        
def simple_frequency_spectrum(x):
    """
    Very simple calculation of frequency spectrum with no detrending,
    windowing, etc. Just the first half (positive frequency components) of
    abs(fft(x))
    """
    spec = numpy.absolute(numpy.fft.fft(x))
    spec = spec[:len(x)/2] # take positive frequency components
    spec /= len(x)         # normalize
    spec *= 2.0            # to get amplitudes of sine components, need to multiply by 2
    spec[0] /= 2.0         # except for the dc component
    return spec

class TuningCurve(object):
    """Class to facilitate working with tuning curves."""

    def __init__(self, D=None):
        """
        If `D` is a dict, it is used to give initial values to the tuning curve.
        """
        self._tuning_curves = {}
        self._counts = {}
        if D is not None:
            for k,v in D.items():
                self._tuning_curves[k] = [v]
                self._counts[k] = 1
                self.n = 1
        else:
            self.n = 0

    def add(self, D):
        for k,v  in D.items():
            self._tuning_curves[k].append(v)
            self._counts[k] += 1
        self.n += 1

    def __getitem__(self, i):
        D = {}
        for k,v in self._tuning_curves[k].items():
            D[k] = v[i]
        return D
    
    def __repr__(self):
        return "TuningCurve: %s" % self._tuning_curves

    def stats(self):
        """Return the mean tuning curve with stderrs."""
        mean = {}
        stderr = {}
        n = self.n
        for k in self._tuning_curves.keys():
            arr = numpy.array(self._tuning_curves[k])
            mean[k] = arr.mean()
            stderr[k] = arr.std()*n/(n-1)/numpy.sqrt(n)
        return mean, stderr

    def max(self):
        """Return the key of the max value and the max value."""
        k = _dict_max(self._tuning_curves)
        return k, self._tuning_curves[k]