Functional near-infrared spectroscopy (fNIRS) is a non-invasive optical neuroimaging method used to investigate functional activity of the cerebral cortex evoked by cognitive, visual, auditory and motor tasks, detecting regional changes of oxy- and deoxy-hemoglobin concentration. Accurate estimation of the stimulus-evoked hemodynamic response (HR) from fNIRS signals in order to quantitatively investigate cognitive functions requires to cope with several noise components. Some of them appear as random disturbances (typically tackled through averaging techniques), while others are due to physiological sources, such as heart beat, respiration, vasomotor waves, and are particularly challenging to be dealt with because they lie in the same frequency band of HR. In this work we present a new two-steps methodology for the HR estimation from fNIRS data. The first step is a pre-processing stage where physiological trends in fNIRS data are reduced by exploiting a mathematical model identified from the signal of a reference channel. In the second step, the pre-processed data of the other channels are filtered with a recently presented non-parametric Bayesian approach (Scarpa et al., Optics Express, 2010). The presented method for HR estimation is compared with widely used methods: conventional averaging, band-pass filtering and principal component analysis (PCA). Results on simulated data reveal the ability of the proposed method to improve the accuracy of the estimates of the functional hemodynamic response, as well as the estimate of peak amplitude and latency. Encouraging preliminary results in a representative real data set showing an improvement of contrast to noise ratio are also reported.