Design narrow bandpass filter for signal with high sampling rate
So, first, to put things into perspective: 4kHz is not a high sampling rate these days (add 5 orders of magnitudes, and things become hard).

Your 15 kHz passband doesn't say anything about the complexity of the filter; what counts is the transition width, ie. the distance between pass- and stopband, as well as the attenuation of the stopband. For a slightly specialized answer, see the answers to http://dsp.stackexchange.com/questions/31066/how-many-taps-does-an-fir-filter-need/31077.


  What is the recommended way to perform tasks like this?


Do a filter design, apply filter.


  What bandpass filter is suitable?


Any bandpass filter that suits your needs – which you haven't specified.

What's missing is:


transition width
stop-band attenuation
acceptable ripple
length constraints


What we know is:


You want a linear phase filter (because linear phase means constant group delay, and that means un-broken phase relationships)
Thus, you probably want a FIR that is symmetric in time



  Furthermore​, the filter step response should have close to zero overshoot. 


Gibb's phenomenon is non-negotiable :) so yeah, use a longish filter with a nice rolloff.


  The settling time could therefore be a little larger.


That will be the effect, yes. You could use a filter design tool that allows you to design with a windowing method and use a window that suits your application well – which I don't know.


  I am more interested in a general solution that delivers solid results using standard IIR or FIR filters (e. g. as available in Python's SciPy library) 


Well, yeah, that's what I'd recommend, but you say just shortly above that you don't want to do a filter design? I'm a bit conflicted. Anyway, there's a lot of functions that will give you a proper design.

I'd recommend the following: 


Design a low-pass FIR filter
convert it to a bandpass filter, by multiplying with a cosine (real-valued data, symmetric filter) or $e^{j2\pi\frac{f_{center}}{f_s}\cdot n}$ (complex-valued data, one-sided filter) as needed.


Designing is easy, something along (untested, straight from the back of my head)

from scipy.signal import fir_filter_design
from math import cos, pi

f_center=ToBeDefined!!!
f_cut = 15
f_s = 4e3
f_rel = f_s/2/f_cut
transition_width = ToBeDefined!!! # e.g. 10 
trans_rel = f_s/2/transition_width
attenuation = ToBeDefined!!! # e.g. 10.0**-5

num_of_taps = estimate_by_reading_my_answer_above(f_rel, trans_rel, attenuation)

lowpass_taps = fir_filter_design.firwin2(numtaps=num_of_taps, cutoff=f_rel, window = "hamming"|"blackmanharris"|"hann"|"chebwin"|…)

bandpass_taps = [lp_tap * cos(2*pi*f_center/f_s*n) for n, lptap in enumerate(lowpass_taps)]



  Do I need to downsample the signal before filtering. If yes, what type of additional low pass filter should I use to avoid aliasing?


Well, considering you won't need most of your signal, yes, that seems advantageous. Typically, you'd just filter away (and decimate in the same step) an 1/N of your initial sampling rate. For example, if your bandpass lies about 100 Hz - 15 / 2 Hz to 100 Hz + 15 / 2 Hz, you won't need anything above 200 Hz – so just decimate to 1/20 of your input sampling rate, with a low pass filter that has 1/40 input sampling rate transition width.

But again, 4kHz is laughably little, and a couple thousand taps won't stress your PC at all. If in doubt,

from scipy import signal
signal.convolve(your_input_signal, filter_taps, method="fft")


does a fast convolution. (you can also just omit the method, or set it to auto, because scipy will just pick fast convolution for long filters, automatically). For a bit of an impression on what's possible with halfway-optimized code on a PC, see my answer to how to implement signal generation on a GPU (and why not); there, I showed that you can do a 107-tap filter at 160 MS/s on my PC back then, so you should be able to do a let's say 700 tap filter at easily 20 MS/s (just a wild guess). That's only 5 thousand times as fast as your sampling.