I have a signal that is sampled at 2KHz. The frequency band of interest is from 1/60 Hz to 1 Hz, and I want the final signal to be downsampled to 10 Hz. What is the best way to do this? I read that narrow band filtering can be unstable? The way that I am doing now is to bandpass my signal (1/60 Hz to 1 Hz) using the 2nd order Butterworth filter from the 'scipy' package and then downsample the filtered signal to 10 Hz. Is this the correct way to do this?

  • $\begingroup$ Why do you want to use an analog-originating filter prototype here? It might make a lot of sense to use Butterworth after decimating, but why for the decimating step? $\endgroup$
    – mmmm
    Jul 19 at 1:21
  • $\begingroup$ No. You should first downsample the signal to 10 Hz and then pass it into a bandpass filter. $\endgroup$
    – ZR Han
    Jul 19 at 1:22
  • $\begingroup$ @ZR HAN Can you please tell me why? $\endgroup$ Jul 19 at 1:54
  • $\begingroup$ Because the filtering effort grows with the ratio of total bandwidth to width of your transition from pass to stop band. By reducing the sample rate first, you make it much easier to run a sharp filter. $\endgroup$
    – mmmm
    Jul 19 at 2:39
  • $\begingroup$ @ZRHan: if you down sample without low-pass filtering first you risk significant aliasing. $\endgroup$
    – Hilmar
    Jul 19 at 12:36

Because your passband is too narrow, leading to zeros and poles very close to the unit circle. Quantization error finally results in poles outside the unit circle and consequently an unstable filter.

First design a bandpass filter with $f_0 = 1/60$ Hz, $f_1 = 1$ Hz and $f_s = 10$ Hz.

from scipy.signal import butter
from scipy.signal import freqz
from scipy.signal import tf2zpk
import numpy as np
import matplotlib.pyplot as plt

N = 10
fs = 10
f1 = 1/60
f2 = 1
Wn = np.array([f1, f2]) / (fs / 2)
b, a = butter(N, Wn, 'bandpass')
w, h = freqz(b, a, 512, fs=fs)
plt.semilogx(w, 20 * np.log10(abs(h)))

enter image description here The result shows that the designed filter gives a frequency response that meets your demand.

Let's have a look at the poles and we can notice that the norm of some poles are close to the unit circle.

z, p, _ = tf2zpk(b, a)
print('abs(p) =', abs(p))
abs(p) = [0.91458903 0.91458903 0.76746989 0.76746989 0.65197236 0.65197236
 1.13021488 1.13021488 1.08513426 1.08513426 1.00506678 1.00506678
 0.91213724 0.91213724 0.85982491 0.85982491 0.57031723 0.57031723
 0.52727447 0.52727447]

Now increase the sampling rate to 2 kHz and we can see that the frequency response is not what we desire any more.

enter image description here

Some of the poles lie outside the unit circle due to the quantization error.

abs(p) = [1.39139451 1.36932567 1.36932567 1.3073682  1.3073682  1.21665313
 1.21665313 1.11080068 1.11080068 1.00253899 1.00253899 0.9029426
 0.9029426  0.82070421 0.82070421 0.76114105 0.76114105 0.72592978
 0.72592978 0.71440178]

Btw, when you downsample the signal 10 Hz, you may apply it by several stages, such as 2000 -> 1000 -> 500 -> 250 -> 50 -> 10. This is because downsampling also requires a low-pass (anti-aliasing) filter and if you decimate the signal by 200 at one time, it also suffers the quantization problem.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.