In functional magnetic resonance imaging (fMRI), local brain activation temporally changes blood oxygenation, which induces a blood oxygenation level dependent (BOLD) contrast in MR images [1]. Given a model of the BOLD response to a stimulus pattern, statistics can be used to quantify the match between the predicted and measured signals in a voxel, and significant activation is assessed via hypothesis testing. Statistical parametric mapping (SPM) estimates parameters of the noise distribution in every voxel to determine a threshold for the computed statistic [2]. It uses the general linear model (GLM): the response to a stimulus pattern is modelled as the output of a linear, time invariant (LTI) system [3]. The response to one type of stimulus is modelled by convolving its temporal distribution with the response function of that type of stimulus. The total response is the sum of responses to all individual stimulus types. The stimulus pattern is known from the experimental setup. However, the haemodynamic response function (HRF), i.e., the temporal BOLD impulse response, is unknown. Essentially, it is a smooth curve that starts to rise two seconds after the stimulus, peaks after six seconds, and returns to baseline within 30 seconds.

The resolution of discrete fMRI time samples yields a coarse description of the HRF, which is a problem for designs with short and/or randomised stimuli. In slice-wise MRI acquisition, there is no optimal sampling for all slices, and slight differences in temporal onset between voxels in different slices must be accurately modelled to obtain good sensitivity. This requires a continuous-time HRF model. It is sometimes preferable to use a common HRF for many experiments and many regions. In many cases though, a common HRF cannot be assumed [4, 5]. Two ways to solve this problem are (*i*) including a HRF model with more basis functions [5], and (*ii*) estimating the shape of the HRF in the statistical analysis [4, 6]. A drawback of the first and simple approach is that it covers only limited variation in the HRF shape, while more advanced approaches suffer from high complexity and long computation times. For previous work on the variability of the BOLD response between brain areas and even across repeated measurements of the same subject, see Neumann et al. [7].

Extracting a HRF from fMRI data requires assumptions about its shape, and is computationally expensive. A simple method described in the literature is selective averaging with a long interstimulus interval (ISI), assuming non-overlapping responses [8, 9]. Trials *with* overlapping responses are also averaged, ignoring the fact that the overlaps introduce errors [3, 10, 11]. Other studies use a functional description of the HRF, whose parameters are determined by curve fitting (examples in [12–14]), one of which [12] uses frequency-domain deconvolution. Another technique based on the GLM expands the HRF into a set of basis functions [15, 16].

Ciuciu et al. [4] use a Bayesian method to extract the HRF. The method can simultaneously extract multiple HRFs from multiple experiments. Woolrich et al. [6] use a fully Bayesian approach to fMRI modelling, including the HRF. Both these approaches use certain prior assumptions for the shape of the HRF: in the first case, causality, smoothness, and starting and ending at baseline level, and Gaussian temporal autocorrelations; in the second case, a set of predefined prior HRF shapes, and many priors for modelling the noise distribution and autocorrelation.

The HRF extraction method described here is *data-driven* instead of model-driven. It is based on Fourier-wavelet regularised deconvolution, *ForWaRD* for short, developed originally for denoising and deblurring of images [17]. *ForWaRD* combines frequency-domain deconvolution for identifying overlapping signals, frequency-domain regularisation for suppressing noise, and wavelet-domain regularisation for separating signal and noise. It is related to recent wavelet-based deconvolution techniques [18–20], with the advantage that the roles of signal (for fMRI: sparse, high frequency) and response (for fMRI: smooth, low frequency) can be interchanged: unlike other wavelet- and vaguelette-based deconvolution methods, *ForWaRD* does the first step entirely in the frequency domain.

The method described in this paper uses an fMRI data set and the stimulus time pattern, and produces a time series of image volumes, which contains the HRF in each voxel. Compared to the simple extraction methods above, this method has the advantage of taking overlapping responses into account. On the other hand, it is much simpler than the Bayesian ones described above, in that it does not rely on shape assumptions/priors of the HRF: the extracted time points are determined only by the fMRI signal and the stimuli: the method is data-driven. This important property means that the extracted HRF is not biased by any a *priori* model. The only requirements for *ForWaRD* are (*i*) LTI responses and (*ii*) separability of signal and noise in the frequency/wavelet domain. Studies of the linearity of the BOLD response indicate that (*i*) is a reasonable assumption, and given the smoothness of the HRF and the availability of smooth wavelet basis functions, (*ii*) is satisfied as well. Another important advantage is the fact that wavelet-based signal processing methods are much better in preserving the 1/*f*-type temporal autocorrelations that are found in fMRI time signals than time-domain methods [21]. Models that do not preserve this structure yield biased results, leading to an increase in false positives [21, 22]. Finally, our method is fast: HRF extraction from a single fMRI data set takes about the same time as spatial resampling.

To avoid possible misunderstandings we would like to stress that our method only concerns the relation between stimulus and haemodynamic response, but does not make any statement about the underlying neural response. Although it is known from the work by Logothetis et al. [23] that the BOLD response in the visual cortex of the monkey brain roughly corresponds to a convolution of the neuronal local field potential (LFP) with some neural response function (NRF), the exact form of this coupling is still unknown and remains an area of active research. In our method, we describe the haemodynamic response as a convolution of the *stimulus* (not the LFP) with an impulse response function, i.e., the HRF. Given the stimulus pattern our method can recover the HRF, but we do not claim to be able to extract the NRF via deconvolution.

The output time series can be used to determine subject-specific, group-specific, and region-specific HRFs, *e.g*., by averaging time signals in the corresponding regions. We demonstrate this in an fMRI experiment using two acquisitions, where the HRF extracted from the first fMRI time series is used to predict responses in a second experiment. A functional description for the HRF is found by fitting a combination of two gamma density functions with variable parameters for height, dilation, peak location and lag. We show that contrast and localisation improve by using the extracted HRF instead of the canonical HRF from the SPM2 program (a sum of two gamma density functions (GDFs) with fixed parameters).

The remainder of this paper is organised as follows. Section 2.2 reviews the *ForWaRD* method for regularised deconvolution. Sect. 2.2.4 describes how *ForWaRD* is used for HRF extraction, and it describes how the method was tested on simulated noisy time series and demonstrates a possible application using event-related fMRI data. Sect. 4 contains some general conclusions.