code: def spectral_decomposition(
      hrf2decompose, time=None, delta=None, ncomp=2
  ):
      """PCA decomposition of symbolic HRF shifted over time
  
      Perform a PCA expansion of a symbolic HRF, time shifted over the
      values in delta, returning the first ncomp components.
  
      This smooths out the HRF as compared to using a Taylor series
      approximation.
  
      Parameters
      ----------
      hrf2decompose : sympy expression
          An expression that can be lambdified as a function of 't'. This
          is the HRF to be expanded in PCA
      time : None or np.ndarray, optional
          None gives default value of np.linspace(-15,50,3251) chosen to
          match fMRIstat implementation.  This corresponds to a time
          interval of 0.02.  Presumed to be equally spaced.
      delta : None or np.ndarray, optional
          None results in default value of np.arange(-4.5, 4.6, 0.1)
          chosen to match fMRIstat implementation.
      ncomp : int, optional
          Number of principal components to retain.
  
      Returns
      -------
      hrf : [sympy expressions]
          A sequence length `ncomp` of symbolic HRFs that are the
          principal components.
      approx :
          TODO
      """
      if time is None:
          time = np.linspace(-15, 50, 3251)
      dt = time[1] - time[0]
      if delta is None:
          delta = np.arange(-4.5, 4.6, 0.1)
      # make numerical implementation from hrf function and symbol t.
      # hrft returns function values when called with values for time as
      # input.
      hrft = lambdify_t(hrf2decompose(T))
      # Create stack of time-shifted HRFs.  Time varies over row, delta
      # over column.
      ts_hrf_vals = np.array([hrft(time - d) for d in delta]).T
      ts_hrf_vals = np.nan_to_num(ts_hrf_vals)
      # PCA
      U, S, V = npl.svd(ts_hrf_vals, full_matrices=0)
      # make interpolators from the generated bases
      basis = []
      for i in range(ncomp):
          b = Interp1dNumeric(
              time, U[:, i], bounds_error=False, fill_value=0.0
          )
          # normalize components with integral of abs of first component
          if i == 0:
              d = np.fabs((b(time) * dt).sum())
          b.y /= d
          basis.append(b)
      # reconstruct time courses for all bases
      W = np.array([b(time) for b in basis]).T
      # regress basis time courses against original time shifted time
      # courses, ncomps by len(delta) parameter matrix
      WH = np.dot(npl.pinv(W), ts_hrf_vals)
      # put these into interpolators to get estimated coefficients for any
      # value of delta
      coef = [
          Interp1dNumeric(delta, w, bounds_error=False, fill_value=0.0)
          for w in WH
      ]
      # swap sign of first component to match that of input HRF.  Swap
      # other components if we swap the first, to standardize signs of
      # components across SVD implementations.
      if coef[0](0) < 0:  # coefficient at time shift of 0
          for i in range(ncomp):
              coef[i].y *= -1.0
              basis[i].y *= -1.0
  
      def approx(time, delta):
          value = 0
          for i in range(ncomp):
              value += coef[i](delta) * basis[i](time)
          return value
  
      approx.coef = coef
      approx.components = basis
      (
          approx.theta,
          approx.inverse,
          approx.dinverse,
          approx.forward,
          approx.dforward,
      ) = invertR(delta, approx.coef)
      # construct aliased functions from bases
      symbasis = []
      for i, b in enumerate(basis):
          symbasis.append(
              implemented_function("%s%d" % (str(hrf2decompose), i), b)
          )
      return symbasis, approx
  
file_path: nipy/nipy/modalities/fmri/fmristat/hrf.py
imports: import numpy
  import numpy.linalg
  from sympy.utilities.lambdify import implemented_function
  from ..None import hrf
  from ..utils import Interp1dNumeric, T, lambdify_t
  from .invert import invertR
label_desc: Use image node_id 1 for calling a global function with example usage: spectral_decomposition(hrf2decompose, time, delta, ncomp) and returns: symbasis,  approx, value