Hemodynamic Response Function

hrf.py

# Thanks to http://www.jarrodmillman.com/rcsds/lectures/convolution_background.html
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

def hrf(x):
    """ Return values for HRF at given times """
    # Gamma pdf for the peak
    peak_values = gamma.pdf(x, 6)
    # Gamma pdf for the undershoot
    undershoot_values = gamma.pdf(x, 12)
    # Combine them
    values = peak_values - 0.35 * undershoot_values
    # Scale max to 0.6
    return values / np.max(values) * 0.6

# TR = Repetition Time
TR = 1.
# Time vector
x = np.arange(0, 30, TR)
# Start plot
fig, ax = plt.subplots()
ax.plot(x, hrf(x), 'darkslategrey', marker='o')
ax.set_xlabel('Time (s)')
ax.set_ylabel('HRF model of signal')
plt.show()