Notebook Fourier transform deconvolution -------------------------------

In [4]:

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
# import mne.filter

Let's do inline plotting:

In [5]:

%matplotlib inline

Let's first create and plot the true pupil response to a transient input. This is called an "Inpulse Response Function" (IRF).

In [6]:

# 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 
    
    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
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.')

Out[6]:

<matplotlib.text.Text at 0x112a30cd0>

Let's simulate convolved timeseries data based this IRF.

In [7]:

# input:
nr_events = 20
transient_times = np.cumsum(np.repeat(6, nr_events)) # in seconds
ramp_durs = np.repeat(1.5, nr_events) # in seconds
ramp_times = transient_times - ramp_durs
ramp_heights = 1 / ramp_durs / 2.0
duration = transient_times[-1] + 10 # in seconds
input_signal_transient = np.zeros(duration*sample_rate)
input_signal_ramp = np.zeros(duration*sample_rate)
for i in range(len(transient_times)):
    t_time = int(transient_times[i]*sample_rate)
    r_time = int(ramp_times[i]*sample_rate)
    r_height = ramp_heights[i]
    input_signal_transient[t_time] = 1
    input_signal_ramp[r_time:t_time] = np.linspace(0,r_height,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.125,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.')

Out[7]:

<matplotlib.text.Text at 0x11555a2d0>

Let's epoch the data and compute mean response:

In [8]:

# 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.')

Out[8]:

<matplotlib.text.Text at 0x1156caed0>

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:

In [9]:

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

Out[9]:

<matplotlib.text.Text at 0x1156ff1d0>

In [13]:

# filter the data:
def butter_lowpass(lowcut, fs, order=3):
    nyq = 0.5 * fs
    low = lowcut / nyq
    b, a = signal.butter(order, low, btype='low')
    return b, a

def butter_lowpass_filter(data, lowcut, fs, order=3):
    b, a = butter_lowpass(lowcut, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y

lowcut = 2.0 * sample_rate
fs = convolved_signal_noise.shape[0] / sample_rate
convolved_signal_filtered = butter_lowpass_filter(convolved_signal_noise, lowcut, fs, order=3)

convolved_signal_filtered[convolved_signal==0] = 0


# Fp = 1.0
# convolved_signal_filtered = mne.filter.low_pass_filter(convolved_signal_noise, Fs=sample_rate, Fp=Fp)

# plt.plot(convolved_signal_noise)
plt.figure()
plt.plot(convolved_signal_noise)
plt.plot(convolved_signal_filtered, 'r', lw=2)
plt.xlabel('time (s)')
plt.xlim(xmax=120)
plt.ylabel('a.u.')
# plt.legend(['true convolved signal', 'filtered signal, {} hz'.format(Fp)])

# freqs = np.fft.fftfreq(convolved_signal_filtered.shape[0], 1.0/sample_rate)
# plt.figure()
# plt.plot(freqs, abs(np.fft.fft(convolved_signal_filtered)))

Out[13]:

<matplotlib.legend.Legend at 0x115cf8150>

In [14]:

# 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_filtered[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.')

Out[14]:

<matplotlib.text.Text at 0x112ac8150>

In [82]:

convolved_signal_fft = np.fft.fft(mean_response)
IRF_fft = np.fft.fft(IRF, mean_response.shape[0])
D = np.fft.ifft(convolved_signal_fft / IRF_fft)

plt.figure()
timepoints = np.linspace(0,mean_response.shape[0]/mean_response.shape[0],mean_response.shape[0])
plt.plot(timepoints, mean_input, '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.')

Out[82]:

<matplotlib.text.Text at 0x114680cd0>

In [83]:

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.')

Out[83]:

<matplotlib.text.Text at 0x114907d90>

In [83]:

In [ ]: