The characterization of the spatiotemporal hemodynamic response (stHR) in the brain is important for understanding the interaction between neighboring brain voxels and regions. In this paper, we design an identification algorithm for the characterization of the cerebral stHR which is modeled by a system of coupled hyperbolic partial differential equation (PDE) and infinite-dimensional ordinary differential equation (ODE). The proposed algorithm provides estimates of the hemodynamic variables (cerebral blood flow and mass density contributed by blood) and physiological parameters using non-invasive Blood Oxygenation Level Dependent (BOLD) data measured with functional Magnetic Resonance Imaging (fMRI) modality. The proposed solution concept follows three main steps: (i) discretization of the stHR model using Galerkin-based finite element method; (ii) estimation of the output derivative using high-order sliding mode differentiator; and (iii) estimation of the state, input, and parameters from sampled-in-space measurements using the reduced-order approximation model and a constrained extended Kalman filter with unknown input algorithm. In addition, sufficient conditions that depend on the chosen discretization scheme, and which guarantee the structural identifiability of the input and parameters, and also the observability of the system are provided. The performance of the proposed algorithm is assessed using both synthetic and real data. The set of the used real data represents the 1-D BOLD signal collected from the visual cortex and acquired in 3 Tesla fMRI scanner.