1 Basic princeples

To design a Butterworth band-pass filter, we can combine a low-pass filter and a high-pass filter as below, $$ H(z) = H_{lp}(z) \cdot H_{hp}(z) \tag{1}. $$ In previous blogs, we’ve gone over the Butterworth low-pass and high-pass filters. Here, we use the designed filters to construct a band-pass filter. You can also design a Butterworth band-stop filter according to the designed low-pass and high-pass filters if necessary. $$ \begin{cases} H_{lp}(z) &= H_{\omega_2}(z)\\
H_{hp}(z) &= H_{\omega_1}(z) \end{cases} \tag{2}. $$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

import numpy as np
import matplotlib.pyplot as plt
from obspy import read
from scipy import fft
from scipy.signal import butter, filtfilt

def bandpass(d, fs, f1, f2, order=4, zerophase=False, taper=0.01):
    n = len(d)
    f = np.arange(n) * fs / n
    w = 2 * np.pi * f / fs
    z = np.exp(1j*w)
    w1 = 2 * np.pi * f1 / fs
    w2 = 2 * np.pi * f2 / fs
    k = np.arange(2*order)
    if order % 2 == 0:
        q = np.exp(1j*(k+0.5)/order*np.pi)
    else:
        q = np.exp(1j*k/order*np.pi)
    q1 = np.tan(w1/2) * q
    p1 = (1 + q1) / (1 - q1)
    h1 = np.ones_like(w, dtype=complex)
    for pp in p1[abs(p1)<1]:
        h1 *= ((1+pp)  / (1-pp/z))
    h1 *= ( (1-1/z) ** order / 2**order )
    q2 = np.tan(w2/2) * q
    p2 = (1 + q2) / (1 - q2)
    h2 = np.ones_like(w, dtype=complex)
    for pp in p2[abs(p2)<1]:
        h2 *= ((1-pp)  / (1-pp/z))
    h2 *= ( (1+1/z) ** order / 2**order )
    h = h1 * h2
    dd = fft.ifft(fft.fft(d) * h).real
    if zerophase:
        dd = fft.ifft(fft.fft(dd[::-1]) * h)[::-1].real
    ni = int(n*taper)
    k1 = np.arange(ni)
    k2 = np.arange(-ni, 0)
    ta1 = np.cos(k1*np.pi/2/ni) ** order
    ta2 = np.cos(k2*np.pi/2/ni) ** order
    dd[:ni] *= ta2; dd[-ni:] *= ta1
    return f, h, dd

def main():
    tr = read()[0]
    tr.data -= np.mean(tr.data)
    tr.detrend(); tr.taper(0.01)
    d1 = tr.data.copy(); dt = tr.stats.delta; fs = 1 / dt
    n = len(d1); t = np.arange(n) * dt
    f1 = 1.5; f2 = 2.5; N = 4
    d1 = tr.data.copy()
    f, h, d2 = bandpass(d1, fs, f1, f2, order=N, zerophase=1, taper=0.01)
    f, h, d4 = bandpass(d1, fs, f1, f2, order=N**2, zerophase=1, taper=0.01)
    [b, a] = butter(N, [2*f1/fs, 2*f2/fs], 'bandpass')
    d3 = filtfilt(b, a, d1)

    plt.figure(figsize=(8, 6))
    plt.subplot(211)
    plt.plot(t, d2, 'r', lw=2, label='this')
    plt.plot(t, d3, 'b', lw=1, label='scipy')
    plt.legend(fontsize=15)
    plt.title('Same filter orders N (4)', fontsize=15)
    plt.subplot(212)
    plt.plot(t, d4, 'r', lw=2, label='this')
    plt.plot(t, d3, 'k', lw=1, label='scipy')
    plt.legend(fontsize=15)
    plt.title('Filter orders this: $N^2$ (16)  scipy: N (4)', fontsize=15)
    plt.tight_layout()
    plt.show()
    
if __name__ == '__main__':
    main()

Band-Wave