Design narrow bandpass filter for signal with high sampling rate

Question

I have to bandpass-filter a signal which has been sampled with 4000 Hz.

Only frequencies around 15 Hz shall remain after filtering (let's say in a band between 10 Hz and 20 Hz; the narrower the better).

I have several questions here:

• What is the recommended way to perform tasks like this?

• What bandpass filter is suitable?

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

In my case, it is important that the phase of the original signal will not be distorted or shifted in the whole process. Furthermore​, the filter step response should have close to zero overshoot. The settling time could therefore be a little larger.

I do not want to perform a sophisticated filter design which delivers optimal results for exactly this specific signal. 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) which could be reused afterwards for similar tasks.

---

UPDATE:

After the answer of @MarcusMüller and the provided links, I basically tried every FIR filter design method available in Python's SciPy library and deciced to go with remez. I developed the following code which may be used for further discussion:

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import lfilter, remez

F_test = 20.0
duration = 10.0
fs = 8000
samples = int(fs*duration)
t = np.arange(samples) / fs
signal_test = (5.0 * t * np.sin(2.0*np.pi*F_test*t)) + (0.5 * np.sin(2.0*np.pi*5.0*t)) + (0.5 * np.sin(2.0*np.pi*100.0*t))

#design filter

ntaps = 5000
edges = [0, F_test - 5.0, F_test - 2.5, F_test + 2.5, F_test + 5.0, 0.5 * fs]
taps = remez(ntaps, edges, [0, 1, 0], Hz=fs, maxiter=2500)

#apply filter
signal_test_filtered = lfilter(taps, 1, signal_test)

#create plot
fig = plt.figure()
ax0 = fig.add_subplot(111)
ax0.plot(t, signal_test_filtered, label='signal_test_filtered')
ax0.set_xlabel("time [s]")
ax0.legend()
fig.show()

Answer 1

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:

1. Design a low-pass FIR filter
2. 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.