View source for Scipy.signal: Butterworth low, high and band-pass

Sometimes we need to remove of frequency range from a time series. For this we can use a Butterworth filter [https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html scipy.signal.butter] and the [https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html scipy.signal.filtfilt] command.

Questions to [mailto:davrot@uni-bremen.de David Rotermund]
{| class="wikitable"
|[https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html scipy.signal.filtfilt]
| style="text-align: center;" |Apply a digital filter forward and backward to a signal.
|-
|[https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html scipy.signal.butter]
| style="text-align: center;" |Butterworth digital and analog filter design.
|}

== Example data ==
<syntaxhighlight lang="python">import numpy as np
import matplotlib.pyplot as plt

samples_per_second: int = 1000
dt: float = 1.0 / samples_per_second

# 10 secs
t: np.ndarray = np.arange(0, int(10 * samples_per_second)) * dt

f_low: float = 1  # Hz
f_mid: float = 10  # Hz
f_high: float = 100  # Hz

sin_low = np.sin(2 * np.pi * t * f_low)
sin_mid = np.sin(2 * np.pi * t * f_mid)
sin_high = np.sin(2 * np.pi * t * f_high)

plt.figure(1)
plt.plot(t, sin_low)
plt.plot(t, sin_mid + 3)
plt.plot(t, sin_high + 6)
plt.xlabel("Time [s]")
plt.ylabel("Waveform shifted")
plt.title("unfiltered data")
plt.show()</syntaxhighlight>[[File:11 1.png]]<div class="figure">
</div>

== Low pass ==
<syntaxhighlight lang="python">from scipy import signal

lowpass_frequency: float = 2.0  # Hz

# Nint : The order of the filter.
# Wn : The critical frequency or frequencies. For lowpass and highpass filters, Wn is a scalar; for bandpass and bandstop filters, Wn is a length-2 sequence.
# For a Butterworth filter, this is the point at which the gain drops to 1/sqrt(2) that of the passband (the “-3 dB point”).
# For digital filters, if fs is not specified, Wn units are normalized from 0 to 1, where 1 is the Nyquist frequency (Wn is thus in half cycles / sample and defined as 2*critical frequencies / fs). If fs is specified, Wn is in the same units as fs.
# For analog filters, Wn is an angular frequency (e.g. rad/s).
# btype{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’}, optional
# fs float, optional : The sampling frequency of the digital system.
b_low, a_low = signal.butter(
    N=4, Wn=lowpass_frequency, btype="lowpass", fs=samples_per_second
)
sin_low_lp = signal.filtfilt(b_low, a_low, sin_low)
sin_mid_lp = signal.filtfilt(b_low, a_low, sin_mid)
sin_high_lp = signal.filtfilt(b_low, a_low, sin_high)

plt.figure(2)
plt.plot(t, sin_low_lp)
plt.plot(t, sin_mid_lp + 3)
plt.plot(t, sin_high_lp + 6)
plt.xlabel("Time [s]")
plt.ylabel("Waveform shifted")
plt.title(f"{lowpass_frequency} Hz low-pass filtered data")
plt.show()</syntaxhighlight>[[File:11 2.png]]<div class="figure">
</div>

== High pass ==
<syntaxhighlight lang="python">from scipy import signal

highpass_frequency: float = 20.0  # Hz

b_high, a_high = signal.butter(
    N=4, Wn=highpass_frequency, btype="highpass", fs=samples_per_second
)
sin_low_hp = signal.filtfilt(b_high, a_high, sin_low)
sin_mid_hp = signal.filtfilt(b_high, a_high, sin_mid)
sin_high_hp = signal.filtfilt(b_high, a_high, sin_high)

plt.figure(3)
plt.plot(t, sin_low_hp)
plt.plot(t, sin_mid_hp + 3)
plt.plot(t, sin_high_hp + 6)
plt.xlabel("Time [s]")
plt.ylabel("Waveform shifted")
plt.title(f"{highpass_frequency} Hz high-pass filtered data")
plt.show()</syntaxhighlight>[[File:11 3.png]]<div class="figure">
</div>

== Band pass ==
<syntaxhighlight lang="python">from scipy import signal

lowpass_frequency: float = 2.0  # Hz
highpass_frequency: float = 20.0  # Hz

b_band, a_band = signal.butter(
    N=4,
    Wn=(lowpass_frequency, highpass_frequency),
    btype="bandpass",
    fs=samples_per_second,
)
sin_low_bp = signal.filtfilt(b_band, a_band, sin_low)
sin_mid_bp = signal.filtfilt(b_band, a_band, sin_mid)
sin_high_bp = signal.filtfilt(b_band, a_band, sin_high)

plt.figure(4)
plt.plot(t, sin_low_bp)
plt.plot(t, sin_mid_bp + 3)
plt.plot(t, sin_high_bp + 6)
plt.xlabel("Time [s]")
plt.ylabel("Waveform shifted")
plt.title(f"({lowpass_frequency} Hz, {highpass_frequency} Hz) band-pass filtered data")
plt.show()</syntaxhighlight>[[File:11 4.png]]<div class="figure">
</div>