In the current paper, we introduce a parametric data-driven model for functional near-infrared spectroscopy that decomposes a signal into a series of independent, rescaled, time-shifted, hemodynamic basis functions. Each decomposed waveform retains relevant biological information about the expected hemodynamic behavior. The model is also presented along with an efficient iterative estimation method to improve the computational speed. Our hemodynamic decomposition model (HDM) extends the canonical model for instances when a) the external stimuli are unknown, or b) when the assumption of a direct relationship between the experimental stimuli and the hemodynamic responses cannot hold. We also argue that the proposed approach can be potentially adopted as a feature transformation method for machine learning purposes. By virtue of applying our devised HDM to a cognitive load classification task on fNIRS signals, we have achieved an accuracy of 86.20%+-2.56% using six channels in the frontal cortex, and 86.34%+-2.81% utilizing only the AFpz channel also located in the frontal area. In comparison, state-of-the-art time-spectral transformations only yield 64.61%+-3.03% and 37.8%+-2.96% under identical experimental settings.