Near Infrared Spectroscopy (NIRS)-based Brain Computer Interface (BCI) systems use feature extraction methods relying mainly on the slope characteristics and mean changes of the hemodynamic responses in respect to certain mental tasks. Nevertheless, spatial patterns across the measurement channels have been detected and should be considered during the feature vector extraction stage of the BCI realization. In this work a Graph Signal Processing (GSP) approach for feature extraction is adopted in order to capture the aforementioned spatial information of the NIRS signals. The proposed GSP-based methodology for feature extraction in NIRS-based BCI systems, namely GNIRS, is applied on a publicly available dataset of NIRS recordings during mental arithmetic task. GNIRS exhibits higher classification rates, up to 92.52%, as compared to the classification rates of two state-of-the-art feature extraction methodologies related to slope and mean values of hemodynamic response, i.e., 90.35% and 82.60%, respectively. In addition, GNIRS leads to the formation of feature vectors with reduced dimensionality in comparison with the baseline approaches. Moreover, it is shown to facilitate high classification rates even from the first second after the onset of the mental task, paving the way for faster NIRS-based BCI systems.