What you want is exactly the -iresp output of 3dDeconvolve. Use the -stim_times option with tent or csplin(b,c,n) where n = c-b+1 to get the point-by-point impulse response between times s+b and s+c. See the 3dDeconvolve pdf manual on the documentation site... best thing I've read in the last three years. :)

Impulse response, overly general and probably unhelpful:
[en.wikipedia.org]

3dDeconvolve is modeling the BOLD signal as the output of a linear system. The raw regressors (i.e. the actual stimulus timing) are the input to the linear system.

Getting into the math a little bit, if you take any linear system and input a delta function, then the output you get is called the Green's function of that system, and contains enough information to describe the entire response of the system. The response of the system to an arbitrary input can be calculated by taking the convolution of that arbitrary input function with the Green's function.

In practical terms, if you have an event-type design where the input is a series of essentially instantaneous events, then the output is a series of some kind of function starting at the times of those events. We call that function the HRF (using terminology specific to this application) or the impulse response (using terminology general to the mathematical form).

Put even another way, the impulse response is the characteristic average response of the system to a single stimulus event. Which is what you're looking for.

OK, so to tie it all together, you are going to want the following syntax:

-stim_times k stim_file_name csplin(0,16,9) \
-iresp k response_output_file_name

Then use 3dmaskave on the response_output_file_name to get the spatially-averaged response

I recommend reading the math part of the manual, then the stim_times option, and the iresp option, and the examples... well that is basically the whole manual I guess.

-dave