import os
import numpy as np
import scipy as sp
import scipy.stats as stats
import scipy.signal as signal
import matplotlib.pyplot as plt
from sympy import *
import math
%matplotlib inline
# define an Impulse Response Function:
def pupil_IRF(timepoints, s=1.0/(10**26), n=10.1, tmax=0.93):
""" pupil_IRF defines the IRF (response to a transient input) of the pupil.
Parameters
----------
t: IRF is defined with respect to 't'
s: scaling factor
n: sets the width
tmax: sets the time to peak
IRF_len: function is evaluated for t = [0:IRF_len]
Returns
-------
y: IRF evaluated for t = [0:IRF_len]
yprime: IRF first derivative evaluated for t = [0:IRF_len]
"""
# in sympy:
t = Symbol('t')
y = ( (s) * (t**n) * (math.e**((-n*t)/tmax)) )
yprime = y.diff(t)
# lambdify:
y = lambdify(t, y, "numpy")
yprime = lambdify(t, yprime, "numpy")
# evaluate:
y = y(timepoints)
yprime = yprime(timepoints)
return (y, yprime)
# create the IRF:
sample_rate = 10.0
IRF_len = 3.0 # in seconds
timepoints = np.linspace(0,IRF_len,IRF_len*sample_rate)
IRF, IRF_prime = pupil_IRF(timepoints=timepoints)
IRF = IRF / IRF.std()
IRF_prime = IRF_prime / IRF_prime.std()
# plot the IRF:
fig = plt.figure()
plt.plot(timepoints, IRF, color='r')
# plt.plot(timepoints, IRF_prime, color='g')
plt.legend(['IRF'])
plt.title('Impulse Response Function')
plt.xlabel('time (s)')
plt.ylabel('a.u.')
<matplotlib.text.Text at 0x110a08e50>
# input:
duration = 35 # in seconds
transient_times = np.array([5, 10, 15, 20, 25, 30]) # in seconds
ramp_durs = np.array([2.4, 2.2, 2.0, 1.8, 1.6, 1.4]) # in seconds
ramp_times = transient_times - ramp_durs
ramp_heights = 1 / (transient_times-ramp_times) / 2.0
input_signal_transient = np.zeros(duration*sample_rate)
input_signal_ramp = np.zeros(duration*sample_rate)
for i in range(len(transient_times)):
input_signal_transient[transient_times[i]*sample_rate] = 1
input_signal_ramp[ramp_times[i]*sample_rate:transient_times[i]*sample_rate] = np.linspace(0,ramp_heights[i], ramp_durs[i]*sample_rate)
input_signal = input_signal_transient + input_signal_ramp
# plot input
timepoints = np.linspace(0,duration,duration*sample_rate)
fig = plt.figure()
plt.plot(timepoints, input_signal, 'r')
plt.legend(['input signal'], loc=2)
plt.title('input signal')
plt.xlabel('time (s)')
plt.ylabel('a.u.')
# convolve inputs with IRF:
convolved_signal = (sp.convolve(input_signal, IRF, 'full'))[:-(IRF.shape[0]-1)]
# plot simulated convolved signal without noise:
fig = plt.figure()
plt.plot(timepoints, input_signal, 'r')
plt.plot(timepoints, convolved_signal)
plt.legend(['input signal', 'simulated convolved signal'], loc=2)
plt.title('simulated convolved signal without noise')
plt.xlabel('time (s)')
plt.ylabel('a.u.')
# let's add some noise:
convolved_signal_noise = convolved_signal + np.random.normal(0,0.5,len(convolved_signal))
# plot simulated convolved signal with noise:
fig = plt.figure()
plt.plot(timepoints, input_signal, 'r')
plt.plot(timepoints, convolved_signal_noise)
plt.legend(['input signal', 'simulated convolved signal'], loc=2)
plt.title('simulated convolved signal with noise')
plt.xlabel('time (s)')
plt.ylabel('a.u.')
<matplotlib.text.Text at 0x110b18710>
Let's epoch the data and compute mean response:
# times for epoching:
epoch_times = [-3, 3] # in seconds
# mean input:
epochs = np.vstack([input_signal[i+(epoch_times[0]*sample_rate):i+(epoch_times[1]*sample_rate)] for i in transient_times*sample_rate])
mean_input = np.mean(epochs, axis=0)
# mean response:
epochs = np.vstack([convolved_signal_noise[i+(epoch_times[0]*sample_rate):i+(epoch_times[1]*sample_rate)] for i in transient_times*sample_rate])
mean_response = np.mean(epochs, axis=0)
# plot mean response versus IRF:
timepoints = np.linspace(epoch_times[0],epoch_times[1],(epoch_times[1]-epoch_times[0])*sample_rate)
fig = plt.figure()
plt.plot(timepoints, mean_input, color='r')
plt.plot(timepoints, mean_response, color='b')
plt.axvline(0, color='k')
plt.legend(['mean response'], loc=2)
plt.title('mean response')
plt.xlabel('time from response (s)')
plt.ylabel('a.u.')
<matplotlib.text.Text at 0x110b51190>
This is uninformative! Due to the sluggishness of our impulse response function, it is hard to see what exactly is driving this mean response. Something seems to be going on already before t=0, but what exactly?
Let's do fft deconvolution!
Given that fft(IRF) has low power at high frequencies, we should filter all the high frequency noise from our measured time series. Let's first determine a sensible cut-off frequency:
freqs = np.fft.fftfreq(int(IRF_len*sample_rate), 1.0/sample_rate)
max_freq = np.ceil(max(freqs))
print max_freq
fig = plt.figure()
plt.plot(freqs, abs(np.fft.fft(IRF)))
plt.axvline(max_freq, color='r')
plt.xlabel('frequencies (hz)')
plt.ylabel('power')
5.0
<matplotlib.text.Text at 0x110afb5d0>
hp = max_freq
N = 2 # Filter order
Wn = hp / (convolved_signal.shape[0] / sample_rate) # Cutoff frequency
B, A = signal.butter(N, Wn, output='ba')
convolved_signal_filtered = signal.filtfilt(B, A, convolved_signal_noise)
# cutoff_freq = 5.0
# FS = 10.0 # sampling rate
# FC = cutoff_freq / (cutoff_freq*sample_rate) # cutoff frequency at 0.05 Hz
# N = 30 # number of filter taps
# a = 1 # filter denominator
# b = signal.firwin(N, cutoff=FC, window='hamming') # filter numerator
# M = FS*60 # number of samples (60 seconds)
# n = np.arange(M) # time index
# convolved_signal_filtered = signal.lfilter(b, a, convolved_signal_noise) # filtered output
# plt.plot(convolved_signal_noise)
plt.figure()
plt.plot(convolved_signal_filtered, 'r')
plt.plot(convolved_signal)
plt.xlabel('time (s)')
plt.ylabel('a.u.')
plt.legend(['filtered signal, {} hz'.format(hp), 'true convolved signal'])
<matplotlib.legend.Legend at 0x110b83250>
convolved_signal_fft = np.fft.fft(convolved_signal_filtered)
IRF_fft = np.fft.fft(IRF, convolved_signal.shape[0])
D = np.fft.ifft(convolved_signal_fft / IRF_fft)
plt.figure()
timepoints = np.linspace(0,duration,duration*sample_rate)
plt.plot(timepoints, input_signal, 'r')
plt.plot(timepoints, D, 'k')
plt.legend(['true input signal', 'fft deconvolved filtered signal'], loc=2)
plt.title('simulated convolved signal with noise')
plt.xlabel('time (s)')
plt.ylabel('a.u.')
<matplotlib.text.Text at 0x110d94650>
convolved_signal_fft = np.fft.fft(convolved_signal)
IRF_fft = np.fft.fft(IRF, convolved_signal.shape[0])
D = np.fft.ifft(convolved_signal_fft / IRF_fft)
plt.figure()
timepoints = np.linspace(0,duration,duration*sample_rate)
plt.plot(timepoints, input_signal, 'r')
plt.plot(timepoints, D, 'k')
plt.legend(['true input signal', 'fft deconvolved convolved signal'], loc=2)
plt.title('simulated convolved signal with noise')
plt.xlabel('time (s)')
plt.ylabel('a.u.')
<matplotlib.text.Text at 0x110df0f90>