Functional magnetic resonance imaging (fMRI) is a key neuroimaging technique. The classic fMRI analysis pipeline is based on the assumption that the hemodynamic response (HR) is the same across brain regions, time, and subjects. Although convenient, there is ample evidence that this assumption does not hold, and that these differences result in inaccuracies in brain activity detection. This article presents a new fMRI processing pipeline that captures the intrinsic intra- and intersubject variability of the HR. At the core of this new pipeline is the definition of a fuzzy hemodynamic response function (HRF). The proposed pipeline includes a new fuzzy general linear model (GLM) able to handle the fuzzy HRF, including a practical realization based on the LR representation of fuzzy numbers. This article also describes how to obtain activation maps from the fuzzy GLM, and a methodology to compute the statistical power of the analysis. The method is evaluated in synthetic and real fMRI data and compared with other state-of-the-art techniques. The experiments based on synthetic data show that the fuzzy GLM approach is more robust under uncertainty regarding the true specific shape of the HR. The experiments based on the real data show an increased volume of the activated brain areas, suggesting that the proposed method is able to prevent false negative errors in the boundaries of target brain regions in which HR should be negligible.