docs: def _hrf_kernel(hrf_model, tr, oversampling=50, fir_delays=None):
      """ Given the specification of the hemodynamic model and time parameters,
      return the list of matching kernels
  
      Parameters
      ----------
      hrf_model : string or None,
          identifier of the hrf model
  
      tr : float
          the repetition time in seconds
  
      oversampling : int, optional
          temporal oversampling factor to have a smooth hrf
  
      fir_delays : list of floats,
          list of delays for finite impulse response models
  
      Returns
      -------
      hkernel : list of arrays
          samples of the hrf (the number depends on the hrf_model used)
      """
      acceptable_hrfs = [
          'spm', 'spm + derivative', 'spm + derivative + dispersion', 'fir',
          'glover', 'glover + derivative', 'glover + derivative + dispersion',
          None]
      if hrf_model == 'spm':
          hkernel = [spm_hrf(tr, oversampling)]
      elif hrf_model == 'spm + derivative':
          hkernel = [spm_hrf(tr, oversampling),
                     spm_time_derivative(tr, oversampling)]
      elif hrf_model == 'spm + derivative + dispersion':
          hkernel = [spm_hrf(tr, oversampling),
                     spm_time_derivative(tr, oversampling),
                     spm_dispersion_derivative(tr, oversampling)]
      elif hrf_model == 'glover':
          hkernel = [glover_hrf(tr, oversampling)]
      elif hrf_model == 'glover + derivative':
          hkernel = [glover_hrf(tr, oversampling),
                     glover_time_derivative(tr, oversampling)]
      elif hrf_model == 'glover + derivative + dispersion':
          hkernel = [glover_hrf(tr, oversampling),
                     glover_time_derivative(tr, oversampling),
                     glover_dispersion_derivative(tr, oversampling)]
      elif hrf_model == 'fir':
          hkernel = [np.hstack((np.zeros(f * oversampling),
                                np.ones(oversampling)))
                     for f in fir_delays]
      elif hrf_model is None:
          hkernel = [np.hstack((1, np.zeros(oversampling - 1)))]
      else:
          raise ValueError('"{0}" is not a known hrf model. Use one of {1}'.
                           format(hrf_model, acceptable_hrfs))
      return hkernel