A simple function for stimulus estimation/reconstruction in a neural system using a Wiener-Kolmogorov filter.

stim_rec.py

##A simple function for stimulus estimation/reconstruction
##in a neural system using a Wiener-Kolmogorov filter
##by Alexandros Kourkoulas-Chondorizos
##v0.3
##This function takes two 1-by-N arrays as input, the input 
##signal presented to the neuron or neural net and the neural
##response. It also takes two integers nfft and tstep, where 
##nfft is the number of data points used in each block for the
##FFT and tstep is the sampling frequency. nfft must be even 
##and a power of 2 is most efficient. tstep is an integer 
##declaring the time step. It creates a WK filter that is then
##convolved with the neural response in order to produce an 
##estimate of the input that produced it. It returns the input
##signal estimate, the zero-centered input signal and the filter.


from scipy import *
from scipy.signal import *
from matplotlib.pyplot import *

def stim_rec(I,firings,nfft,tstep):
    if I.ndim!=1 or firings.ndim!=1 or not(isinstance(nfft,int))\
    or not(nfft%2==0) or not(isinstance(tstep,int)) or not(tstep>0):
        raise ValueError('I and firings should be 1-by-N arrays, nfft an even integer and tstep a positive integer')

    tstep_s=tstep*0.001  #converts to sec
    Fs=1/tstep_s #in Hz
    stim=I-I.mean() #subtracts the mean stimulus
    spk=firings*Fs
    spk=spk-spk.mean() #subtracts mean firing rate and converts to units of spikes/sec
    win=get_window('bartlett',nfft,False)
    #power spectral density of the neural activity
    Pxx,pxx_freqs=psd(spk,nfft,Fs,window=win,noverlap=nfft/2)
    #cross power spectral density of input and neural activity
    Pxy,pxy_freqs=csd(stim,spk,nfft,Fs,window=win,noverlap=nfft/2)

    Tfs=Pxy/Pxx #Estimates the transfer function
    Tfl=zeros([nfft])
    Tfl[:nfft/2],Tfl[nfft/2:]=Tfs[:nfft/2],flipud(Tfs[1:]).conj()
    
    T=real(ifft(Tfl))*nfft #inverse Fourier transform with negative phase

    h=zeros([nfft+1])
    h[:nfft/2],h[nfft/2:nfft+1]=T[(nfft/2):nfft],T[:nfft/2+1] #unwrap the filter

    Iest=convolve(h*tstep_s,spk,'same') #estimate input
    return Iest, stim, h