Helve’s Python memo

Pythonを使った機械学習やデータ分析の備忘録

<SciPy> FIRフィルタによる波形整形

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
カットオフ周波数。

width: float or None
window: string
窓関数を指定する。デフォルトは'hamming' (ハミング窓)。

pass_zero: bool
周波数0(直流成分)が通過するか指定。
デフォルトは"True".

fs: float
信号のサンプル周波数。
cutoff引数は0からfs/2の間にしなければならない。
デフォルト値は2.

戻り値は、長さnumtapsのFIRフィルタの係数配列となる。

フィルタの適用

次に、scipy.signal.lfilter()により、信号にフィルタを適用する。

scipy.signal.lfilter(b, a, x, axis=-1, zi=None)

主な引数の説明は以下の通り。

b: array
分子の係数。

a: array
分母の係数。

x: array
入力信号。

戻り値は、フィルタの出力の配列となる。

この関数では、フィルタが以下の式で表されているものとする。
\displaystyle  H(z)=\frac{\sum_{k=0}^{M} b_k z^{-k} } {\sum_{k=0}^{N}  a_k z^{-k} }
引数bが右辺の分子の係数、引数aが右辺の分母の係数に対応する。
今回、FIRフィルタを扱うので、分母のaを1とおく。


波形整形

基本的なフィルタを信号に適用して、波形や振幅の変化を確認する。
振幅の計算については、以下のページを参照。

helve-python.hatenablog.jp

ローパスフィルタ

カットオフ周波数を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]の高周波成分が減衰・除去されている。
f:id:Helve:20180617193933p:plain

ハイパスフィルタ

カットオフ周波数を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]の低周波成分が除去されている。
f:id:Helve:20180617193937p:plain

バンドパスフィルタ

通過周波数を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]の信号が減衰・除去されいている。
f:id:Helve:20180617193939p:plain

バンドエリミネイトフィルタ

通過周波数を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]の信号が除去されいている。
f:id:Helve:20180617193943p:plain