docs: def S(self):
          """Cross-spectral density.
  
          .. math:: \mathbf{S}(f) = \mathbf{H}(f) \mathbf{C} \mathbf{H}'(f)
          """
          if self.c is None:
              raise RuntimeError('Cross-spectral density requires noise '
                                 'covariance matrix c.')
          H = self.H()
          # TODO: can we do that more efficiently?
          S = np.empty(H.shape, dtype=H.dtype)
          for f in range(H.shape[2]):
              S[:, :, f] = H[:, :, f].dot(self.c).dot(H[:, :, f].conj().T)
          return S