はじめに #
SciPyを使って、FIR (Finite Impulse Response, 有限インパルス応答) フィルタによる離散信号の波形を整形する。ローパス、ハイパス、バンドパス、バンドエリミネイトの各フィルタの設計から、信号への適用まで行う。
環境 #
| ソフトウェア | バージョン |
|---|---|
| Python | 3.6.5 |
| NumPy | 1.14.2 |
| SciPy | 1.0.1 |
元の離散信号 #
次の3つの信号を合成して、フーリエ変換を行う離散信号とする。
- 周波数30[Hz], 振幅3の正弦波
- 周波数50[Hz], 振幅0.3の正弦波
- 周波数120[Hz], 振幅0.2の正弦波
信号のデータ点数は2048, サンプリング周期は0.001[s]とする。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
N = 1024 # サンプル数
dt = 0.001 # サンプリング周期 [s]
f1, f2, f3 = 10, 60, 300 # 周波数 [Hz]
t = np.arange(0, N*dt, dt) # 時間 [s]
x = 3*np.sin(2*np.pi*f1*t) + 0.3*np.sin(2*np.pi*f2*t) + 0.2*np.sin(2*np.pi*f3*t) # 信号
SciPyの関数 #
フィルタの設計 #
まず、関数scipy.signal.firwin()によりフィルタを設計する。
scipy.signal.firwin(
numtaps, cutoff, width=None,
window='hamming', pass_zero=True,
scale=True, nyq=None, fs=None
)
主な引数の説明は以下の通り。
numtaps: int
FIRフィルタの長さ。
cutoff: float or 1D array_like
カットオフ周波数。
ローパス、ハイパスの場合はfloat,
バンドパス、バンドエリミネイトの場合は1D array_likeとする。
cutoffは0からfs/2の間にしなければならない。
window: string
窓関数を指定する。デフォルトは’hamming’ （ハミング窓）。
pass_zero: bool
周波数0（直流成分）が通過するか指定。
デフォルトはTrue.
fs: float
信号のサンプル周波数。デフォルト値は2.
scipy.signal.firwin()の戻り値は、長さnumtapsのFIRフィルタの係数配列となる。
フィルタの適用 #
次に、関数scipy.signal.lfilter()により信号にフィルタを適用する。
scipy.signal.lfilter(b, a, x, axis=-1, zi=None)
主な引数の説明は以下の通り。
b: array
分子の係数。
a: array
分母の係数。
x: array
入力信号。
戻り値は、フィルタを掛けた信号の配列となる。
この関数では、フィルタが以下の式で表されているものとする。
$$ H(z)=\frac{\sum_{k=0}^{M} b_k z^{-k} } {\sum_{k=0}^{N} a_k z^{-k} } $$
引数bが右辺の分子の係数、引数aが右辺の分母の係数に対応する。今回はFIRフィルタを扱うので、分母のaを1とおく。
波形整形 #
基本的なフィルタを信号に適用して、波形や振幅の変化を確認する。
ローパスフィルタ #
カットオフ周波数を40[Hz]としたローパスフィルタ。
filter1 = signal.firwin(numtaps=21, cutoff=40, fs=1/dt)
y1 = signal.lfilter(filter1, 1, x)
F1 = np.fft.fft(y1)
Amp1 = np.abs(F1/(N/2))
60, 300[Hz]の高周波成分が減衰・除去されている。
ハイパスフィルタ #
カットオフ周波数を100[Hz]としたハイパスフィルタ。
filter2 = signal.firwin(numtaps=51, cutoff=100, fs=1/dt, pass_zero=False)
y2 = signal.lfilter(filter2, 1, x)
F2 = np.fft.fft(y2)
Amp2 = np.abs(F2/(N/2))
10, 60[Hz]の低周波成分が除去されている。
バンドパスフィルタ #
通過周波数を30～100[Hz]としたバンドパスフィルタ。
filter3 = signal.firwin(numtaps=51, cutoff=[30, 100], fs=1/dt, pass_zero=False)
y3 = signal.lfilter(filter3, 1, x)
F3 = np.fft.fft(y3)
Amp3 = np.abs(F3/(N/2))
10, 300[Hz]の信号が減衰・除去されいている。
バンドエリミネイトフィルタ #
通過周波数を30[Hz]以下、100[Hz]以上としたバンドエリミネイトフィルタ。
filter4 = signal.firwin(numtaps=31, cutoff=[30, 100], fs=1/dt)
y4 = signal.lfilter(filter4, 1, x)
F4 = np.fft.fft(y4)
Amp4 = np.abs(F4/(N/2))
60[Hz]の信号が除去されている。
