| | 11 | def ccf(x, y, axis=None): |
|---|
| | 12 | """Computes the cross-correlation function of two series `x` and `y`. |
|---|
| | 13 | Note that the computations are performed on anomalies (deviations from |
|---|
| | 14 | average). |
|---|
| | 15 | Returns the values of the cross-correlation at different lags. |
|---|
| | 16 | Lags are given as [0,1,2,...,n,n-1,n-2,...,-2,-1] (not any more) |
|---|
| | 17 | |
|---|
| | 18 | :Parameters: |
|---|
| | 19 | `x` : 1D MaskedArray |
|---|
| | 20 | Time series. |
|---|
| | 21 | `y` : 1D MaskedArray |
|---|
| | 22 | Time series. |
|---|
| | 23 | `axis` : integer *[None]* |
|---|
| | 24 | Axis along which to compute (0 for rows, 1 for cols). |
|---|
| | 25 | If `None`, the array is flattened first. |
|---|
| | 26 | """ |
|---|
| | 27 | assert(x.ndim == y.ndim, "Inconsistent shape !") |
|---|
| | 28 | # assert(x.shape == y.shape, "Inconsistent shape !") |
|---|
| | 29 | if axis is None: |
|---|
| | 30 | if x.ndim > 1: |
|---|
| | 31 | x = x.ravel() |
|---|
| | 32 | y = y.ravel() |
|---|
| | 33 | npad = x.size + y.size |
|---|
| | 34 | xanom = (x - x.mean(axis=None)) |
|---|
| | 35 | yanom = (y - y.mean(axis=None)) |
|---|
| | 36 | Fx = numpy.fft.fft(xanom, npad, ) |
|---|
| | 37 | Fy = numpy.fft.fft(yanom, npad, ) |
|---|
| | 38 | iFxy = numpy.fft.ifft(Fx.conj()*Fy).real |
|---|
| | 39 | varxy = numpy.sqrt(numpy.inner(xanom,xanom) * numpy.inner(yanom,yanom)) |
|---|
| | 40 | else: |
|---|
| | 41 | npad = x.shape[axis] + y.shape[axis] |
|---|
| | 42 | if axis == 1: |
|---|
| | 43 | if x.shape[0] != y.shape[0]: |
|---|
| | 44 | raise ValueError, "Arrays should have the same length!" |
|---|
| | 45 | xanom = (x - x.mean(axis=1)[:,None]) |
|---|
| | 46 | yanom = (y - y.mean(axis=1)[:,None]) |
|---|
| | 47 | varxy = numpy.sqrt((xanom*xanom).sum(1) * (yanom*yanom).sum(1))[:,None] |
|---|
| | 48 | else: |
|---|
| | 49 | if x.shape[1] != y.shape[1]: |
|---|
| | 50 | raise ValueError, "Arrays should have the same width!" |
|---|
| | 51 | xanom = (x - x.mean(axis=0)) |
|---|
| | 52 | yanom = (y - y.mean(axis=0)) |
|---|
| | 53 | varxy = numpy.sqrt((xanom*xanom).sum(0) * (yanom*yanom).sum(0)) |
|---|
| | 54 | Fx = numpy.fft.fft(xanom, npad, axis=axis) |
|---|
| | 55 | Fy = numpy.fft.fft(yanom, npad, axis=axis) |
|---|
| | 56 | iFxy = numpy.fft.ifft(Fx.conj()*Fy,n=npad,axis=axis).real |
|---|
| | 57 | # We juste turn the lags into correct positions: |
|---|
| | 58 | iFxy = numpy.concatenate((iFxy[len(iFxy)/2:len(iFxy)],iFxy[0:len(iFxy)/2])) |
|---|
| | 59 | return iFxy/varxy |
|---|
| 12 | | def record(output, cfilename = 'SpikeTrainPlay.wav', fs=44100, enc = 'pcm26'): |
|---|
| 13 | | """ record the 'sound' produced by a neuron. Takes a spike train as the |
|---|
| 14 | | output. |
|---|
| 15 | | |
|---|
| 16 | | >>> record(my_spike_train) |
|---|
| 17 | | |
|---|
| 18 | | """ |
|---|
| 19 | | |
|---|
| 20 | | |
|---|
| 21 | | # from the spike list |
|---|
| 22 | | simtime_seconds = (output.t_stop - output.t_start)/1000. |
|---|
| 23 | | #time = numpy.linspace(0, simtime_seconds , fs*simtime_seconds) |
|---|
| 24 | | (trace,time) = numpy.histogram(output.spike_times*1000., fs*simtime_seconds) |
|---|
| 25 | | |
|---|
| 26 | | |
|---|
| 27 | | # TODO convolve with proper spike... |
|---|
| 28 | | spike = numpy.ones((fs/1000.,)) # one ms |
|---|
| 29 | | |
|---|
| 30 | | trace = numpy.convolve(trace, spike, mode='same')#/2.0 |
|---|
| 31 | | trace /= numpy.abs(trace).max() * 1.1 |
|---|
| 32 | | |
|---|
| 33 | | from scikits.audiolab import wavwrite |
|---|
| 34 | | wavwrite(trace, cfilename, fs = fs, enc = enc) |
|---|
| 35 | | |
|---|
| 36 | | def play(output): |
|---|
| 37 | | """ |
|---|
| 38 | | plays a spike list to the audio output |
|---|
| 39 | | |
|---|
| 40 | | play(spike_list) where spike_list is a spike_list object |
|---|
| 41 | | |
|---|
| 42 | | see playing_with_simple_single_neuron.py for a sample use |
|---|
| 43 | | |
|---|
| 44 | | >>> play(my_spike_train) |
|---|
| 45 | | |
|---|
| 46 | | TODO: make it possible to play multiple spike trains in stereo |
|---|
| 47 | | """ |
|---|
| 48 | | |
|---|
| 49 | | |
|---|
| 50 | | from tempfile import mkstemp |
|---|
| 51 | | fd, cfilename = mkstemp('SpikeListPlay.wav') |
|---|
| 52 | | try: |
|---|
| 53 | | record(output, cfilename) |
|---|
| 54 | | import pyaudio |
|---|
| 55 | | import wave |
|---|
| 56 | | |
|---|
| 57 | | chunk = 1024 |
|---|
| 58 | | wf = wave.open(cfilename, 'rb') |
|---|
| 59 | | p = pyaudio.PyAudio() |
|---|
| 60 | | |
|---|
| 61 | | # open stream |
|---|
| 62 | | stream = p.open(format = |
|---|
| 63 | | p.get_format_from_width(wf.getsampwidth()), |
|---|
| 64 | | channels = wf.getnchannels(), |
|---|
| 65 | | rate = wf.getframerate(), |
|---|
| 66 | | output = True) |
|---|
| 67 | | |
|---|
| 68 | | # read data |
|---|
| 69 | | data = wf.readframes(chunk) |
|---|
| 70 | | |
|---|
| 71 | | # play stream |
|---|
| 72 | | while data != '': |
|---|
| 73 | | stream.write(data) |
|---|
| 74 | | data = wf.readframes(chunk) |
|---|
| 75 | | |
|---|
| 76 | | stream.close() |
|---|
| 77 | | p.terminate() |
|---|
| 78 | | except: |
|---|
| 79 | | print "Error playing the SpikeTrain " |
|---|
| 80 | | # finally |
|---|
| 81 | | os.remove(cfilename) |
|---|
| 82 | | |
|---|
| 83 | | # Python 2.4 compatibility |
|---|
| 84 | | # finally: |
|---|
| 85 | | os.remove(cfilename) |
|---|
