PythonのSciPyでローパスフィルタをかける！

ローパスフィルタについて

ローパスフィルタとは？

フィルタとは、日本語で言うと濾過器のことをさします。信号処理の世界では、ある信号から必要な成分のみを抽出する時に使います。

ローパスフィルタで重要な因子と用語（図解）

ローパスフィルタは低周波数帯域の信号を通し、高周波帯域の信号を阻止するので、周波数軸で考えると下図のように3つの領域（「通過域」「遷移域」「阻止域」）に分けられます。

フィルタは通過域の信号はほとんど減衰させないものを指しますが、阻止域で直角に減衰するのではなく、遷移域と呼ばれる徐々に減衰量が変わっていく範囲を持ちます。

さらに、通過域の端部である通過域端周波数でもある程度の減衰量を持ちます。

各周波数でどのような減衰特性を持たせるか、どこまで許容させるのかといった要素を決めて行くのが、フィルタ設計の主目的※です。

※フィルタ処理も元の波形を変化させてしまう信号処理技術であることには変わりがないので、減衰量の設計以外にも信号の位相変化を許容するかどうか、安定性はあるかどうか…といった要素もフィルタ設計にはありますが、今回はそこまで細かい話には踏み込まず、工学の現場でよく使われているバターワースフィルタを扱います。

PythonのSciPyローパスフィルタ

ソースコード

メインファイル

ソースコードを以下に示します。フィルタ処理を実際に行う部分はdef文を使って関数化しました。関数の使い方については「Pythonの関数def文の使い方！引数や別ファイル式も解説」の記事に詳細を示しています。

まずはメイン(main_filter.py)のコードです。

サンプルとしてガウシアンノイズを生成していますが、こちらについては「Pythonでガウス分布を持つノイズの作り方と調整方法」の記事を参照下さい。

中身ですが、上の図で説明した通過域端周波数fp[Hz]と阻止域端周波数fs[Hz]を設定します。またそれぞれの周波数における減衰量もgpassとgstopで設定します。

import filter_function
import numpy as np
from matplotlib import pyplot as plt

samplerate = 25600                                   #波形のサンプリングレート
x = np.arange(0, 12800) / samplerate                 #波形生成のための時間軸の作成
data = np.random.normal(loc=0, scale=1, size=len(x)) #ガウシアンノイズを生成

fp = 3000       #通過域端周波数[Hz]
fs = 6000       #阻止域端周波数[Hz]
gpass = 3       #通過域端最大損失[dB]
gstop = 40      #阻止域端最小損失[dB]

data_filt = filter_function.lowpass(data, samplerate, fp, fs, gpass, gstop)

import filter_function

import numpy as np

from matplotlib import pyplot as plt

samplerate = 25600 #波形のサンプリングレート

x = np.arange(0, 12800) / samplerate #波形生成のための時間軸の作成

data = np.random.normal(loc=0, scale=1, size=len(x)) #ガウシアンノイズを生成

fp = 3000 #通過域端周波数[Hz]

fs = 6000 #阻止域端周波数[Hz]

gpass = 3 #通過域端最大損失[dB]

gstop = 40 #阻止域端最小損失[dB]

data_filt = filter_function.lowpass(data, samplerate, fp, fs, gpass, gstop)

関数ファイル

以下の関数を使ってローパスフィルタを通した波形を取得します。ここで、設定した周波数をナイキスト周波数で除していますが、これはデジタルフィルタを使う時に周波数を正規化する必要があるからです。

以下に引用する公式ページに使い方が書いてあります（日本語検索では意味まで説明していた記事は見当たらず…）。

For digital filters, these are normalized from 0 to 1, where 1 is the Nyquist frequency.
SciPy.org:scipy.signal.buttord

僕はここの正規化の使い方、考え方のエビデンスをとるのに時間がかかってしまいました。初心者の壁かも知れません！

1がナイキスト周波数である必要があるので、サンプリングレートの半分であるナイキスト周波数で割っています。公式ページでは周波数が[Hz]ではありませんが、比をとっているだけなので、設定した周波数とナイキスト周波数がともに[Hz]で書いてあれば問題ありませんね。

import numpy as np
from scipy import signal

#バターワースフィルタ（ローパス）
def lowpass(x, samplerate, fp, fs, gpass, gstop):
    fn = samplerate / 2                           #ナイキスト周波数
    wp = fp / fn                                  #ナイキスト周波数で通過域端周波数を正規化
    ws = fs / fn                                  #ナイキスト周波数で阻止域端周波数を正規化
    N, Wn = signal.buttord(wp, ws, gpass, gstop)  #オーダーとバターワースの正規化周波数を計算
    b, a = signal.butter(N, Wn, "low")            #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x)                  #信号に対してフィルタをかける
    return y                                      #フィルタ後の信号を返す

import numpy as np

from scipy import signal

#バターワースフィルタ（ローパス）

def lowpass(x, samplerate, fp, fs, gpass, gstop):

fn = samplerate / 2 #ナイキスト周波数

wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化

ws = fs / fn #ナイキスト周波数で阻止域端周波数を正規化

N, Wn = signal.buttord(wp, ws, gpass, gstop) #オーダーとバターワースの正規化周波数を計算

b, a = signal.butter(N, Wn, "low") #フィルタ伝達関数の分子と分母を計算

y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける

return y #フィルタ後の信号を返す

Wnも同様です。正規化する必要があります。そしてこのWnには-3[dB]ポイントという特別な意味合いを持ちます。バターワースフィルタのWn周波数では減衰量が\(1/\sqrt{2}\)（半値、-3[dB]）になります。

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”).
SciPy.org:scipy.signal.butter

signal.butterで計算しているbとaは伝達関数の分母分子です。そもそもローパスフィルタを始めとした各種フィルタはラプラス変換を利用して伝達関数の形式で表すことができます。伝達関数の分母と分子はそれぞれ多項式になっているため、このb, aはベクトル形式になっています。

参考：フーリエ変換まで含めた全コード

先ほどまではフィルタ処理の部分のみを解説していましたが、ファイルを2つ作ったり面倒だと思いますので、コピペしてすぐ動かせる全コードを以下に示します。

このコードはサンプル波形を生成して、フーリエ変換とデシベル変換を行っており、フィルタ前後の波形を比較する事が可能です。是非動かして遊んでみて下さい。

import numpy as np
from scipy import signal
from scipy import fftpack
from matplotlib import pyplot as plt

#バターワースフィルタ（ローパス）
def lowpass(x, samplerate, fp, fs, gpass, gstop):
    fn = samplerate / 2                           #ナイキスト周波数
    wp = fp / fn                                  #ナイキスト周波数で通過域端周波数を正規化
    ws = fs / fn                                  #ナイキスト周波数で阻止域端周波数を正規化
    N, Wn = signal.buttord(wp, ws, gpass, gstop)  #オーダーとバターワースの正規化周波数を計算
    b, a = signal.butter(N, Wn, "low")            #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x)                  #信号に対してフィルタをかける
    return y                                      #フィルタ後の信号を返す

# オーバーラップ処理
def ov(data, samplerate, Fs, overlap):
    Ts = len(data) / samplerate                                                # 全データ長
    Fc = Fs / samplerate                                                       # フレーム周期
    x_ol = Fs * (1 - (overlap / 100))                                          # オーバーラップ時のフレームずらし幅
    N_ave = int((Ts - (Fc * (overlap / 100))) / (Fc * (1 - (overlap / 100))))  # 抽出するフレーム数（平均化に使うデータ個数）

    array = []                                                                 # 抽出したデータを入れる空配列の定義

    # forループでデータを抽出
    for i in range(N_ave):
        ps = int(x_ol * i)                                                     # 切り出し位置をループ毎に更新
        array.append(data[ps:ps + Fs:1])                                       # 切り出し位置psからフレームサイズ分抽出して配列に追加
    return array, N_ave                                                        # オーバーラップ抽出されたデータ配列とデータ個数を戻り値にする

# 窓関数処理（ハニング窓）
def hanning(data_array, Fs, N_ave):
    han = signal.hann(Fs)                    # ハニング窓作成
    acf = 1 / (sum(han) / Fs)                # 振幅補正係数(Amplitude Correction Factor)

    # オーバーラップされた複数時間波形全てに窓関数をかける
    for i in range(N_ave):
        data_array[i] = data_array[i] * han  # 窓関数をかける

    return data_array, acf

# FFT処理
def fft_ave(data_array, samplerate, Fs, N_ave, acf):
    fft_array = []
    for i in range(N_ave):
        fft_array.append(acf * np.abs(fftpack.fft(data_array[i]) / (Fs / 2)))  # FFTをして配列に追加、窓関数補正値をかけ、(Fs/2)の正規化を実施。

    fft_axis = np.linspace(0, samplerate, Fs)                                  # 周波数軸を作成
    fft_array = np.array(fft_array)                                            # 型をndarrayに変換
    fft_mean = np.mean(fft_array, axis=0)                                      # 全てのFFT波形の平均を計算

    return fft_array, fft_mean, fft_axis

#リニア値からdBへ変換
def db(x, dBref):
    y = 20 * np.log10(x / dBref)     #変換式
    return y                         #dB値を返す

# ここからサンプル波形生成とフィルタ処理をする-------------------------------------------
samplerate = 25600
x = np.arange(0, 12800) / samplerate                    # 波形生成のための時間軸の作成
data = np.random.normal(loc=0, scale=1, size=len(x))    # ガウシアンノイズを生成

fp = 3000                                               # 通過域端周波数[Hz]
fs = 6000                                               # 阻止域端周波数[Hz]
gpass = 3                                               # 通過域端最大損失[dB]
gstop = 40                                              # 阻止域端最小損失[dB]

# ローパスをする関数を実行
data_filt = lowpass(data, samplerate, fp, fs, gpass, gstop)
# ------------------------------------------------------------------------------

# ここから平均化フーリエ変換とデシベル変換を行う----------------------------------------------
Fs = 4096                                               # フレームサイズ
overlap = 90                                            # オーバーラップ率

# 作成した関数を実行：オーバーラップ抽出された時間波形配列
t_array_org, N_ave_org = ov(data, samplerate, Fs, overlap)
t_array_filt, N_ave_filt = ov(data_filt, samplerate, Fs, overlap)

# 作成した関数を実行：ハニング窓関数をかける
t_array_org, acf_org = hanning(t_array_org, Fs, N_ave_org)
t_array_filt, acf_filt = hanning(t_array_filt, Fs, N_ave_filt)

# 作成した関数を実行：FFTをかける
fft_array_org, fft_mean_org, fft_axis_org = fft_ave(t_array_org, samplerate, Fs, N_ave_org, acf_org)
fft_array_filt, fft_mean_filt, fft_axis_filt = fft_ave(t_array_filt, samplerate, Fs, N_ave_filt, acf_filt)

# dB変換する
fft_mean_org = db(fft_mean_org, 2e-5)
fft_mean_filt = db(fft_mean_filt, 2e-5)
# -----------------------------------------------------------------------------------

# ここからグラフ描画---------------------------------------------------------------------
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')
ax2 = fig.add_subplot(212)
ax2.yaxis.set_ticks_position('both')
ax2.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('SP[Pa]')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('SPL[dB]')

# データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。
ax1.plot(x, data, label='original', lw=1)
ax1.plot(x, data_filt, label='filtered', lw=1)
ax2.plot(fft_axis_org, fft_mean_org, label='original', lw=1)
ax2.plot(fft_axis_filt, fft_mean_filt, label='filtered', lw=1)
plt.legend()

# 軸のリミットを設定する。
ax2.set_xlim(0, max(fft_axis_org)/2)
ax2.set_ylim(-20, 100)

# レイアウト設定
fig.tight_layout()

# グラフを表示する。
plt.show()
plt.close()
# -----------------------------------------------------------------------------------

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

import numpy as np

from scipy import signal

from scipy import fftpack

from matplotlib import pyplot as plt

#バターワースフィルタ（ローパス）

def lowpass(x, samplerate, fp, fs, gpass, gstop):

fn = samplerate / 2 #ナイキスト周波数

wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化

ws = fs / fn #ナイキスト周波数で阻止域端周波数を正規化

N, Wn = signal.buttord(wp, ws, gpass, gstop) #オーダーとバターワースの正規化周波数を計算

b, a = signal.butter(N, Wn, "low") #フィルタ伝達関数の分子と分母を計算

y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける

return y #フィルタ後の信号を返す

# オーバーラップ処理

def ov(data, samplerate, Fs, overlap):

Ts = len(data) / samplerate # 全データ長

Fc = Fs / samplerate # フレーム周期

x_ol = Fs * (1 - (overlap / 100)) # オーバーラップ時のフレームずらし幅

N_ave = int((Ts - (Fc * (overlap / 100))) / (Fc * (1 - (overlap / 100)))) # 抽出するフレーム数（平均化に使うデータ個数）

array = [] # 抽出したデータを入れる空配列の定義

# forループでデータを抽出

for i in range(N_ave):

ps = int(x_ol * i) # 切り出し位置をループ毎に更新

array.append(data[ps:ps + Fs:1]) # 切り出し位置psからフレームサイズ分抽出して配列に追加

return array, N_ave # オーバーラップ抽出されたデータ配列とデータ個数を戻り値にする

# 窓関数処理（ハニング窓）

def hanning(data_array, Fs, N_ave):

han = signal.hann(Fs) # ハニング窓作成

acf = 1 / (sum(han) / Fs) # 振幅補正係数(Amplitude Correction Factor)

# オーバーラップされた複数時間波形全てに窓関数をかける

for i in range(N_ave):

data_array[i] = data_array[i] * han # 窓関数をかける

return data_array, acf

# FFT処理

def fft_ave(data_array, samplerate, Fs, N_ave, acf):

fft_array = []

for i in range(N_ave):

fft_array.append(acf * np.abs(fftpack.fft(data_array[i]) / (Fs / 2))) # FFTをして配列に追加、窓関数補正値をかけ、(Fs/2)の正規化を実施。

fft_axis = np.linspace(0, samplerate, Fs) # 周波数軸を作成

fft_array = np.array(fft_array) # 型をndarrayに変換

fft_mean = np.mean(fft_array, axis=0) # 全てのFFT波形の平均を計算

return fft_array, fft_mean, fft_axis

#リニア値からdBへ変換

def db(x, dBref):

y = 20 * np.log10(x / dBref) #変換式

return y #dB値を返す

# ここからサンプル波形生成とフィルタ処理をする-------------------------------------------

samplerate = 25600

x = np.arange(0, 12800) / samplerate # 波形生成のための時間軸の作成

data = np.random.normal(loc=0, scale=1, size=len(x)) # ガウシアンノイズを生成

fp = 3000 # 通過域端周波数[Hz]

fs = 6000 # 阻止域端周波数[Hz]

gpass = 3 # 通過域端最大損失[dB]

gstop = 40 # 阻止域端最小損失[dB]

# ローパスをする関数を実行

data_filt = lowpass(data, samplerate, fp, fs, gpass, gstop)

# ------------------------------------------------------------------------------

# ここから平均化フーリエ変換とデシベル変換を行う----------------------------------------------

Fs = 4096 # フレームサイズ

overlap = 90 # オーバーラップ率

# 作成した関数を実行：オーバーラップ抽出された時間波形配列

t_array_org, N_ave_org = ov(data, samplerate, Fs, overlap)

t_array_filt, N_ave_filt = ov(data_filt, samplerate, Fs, overlap)

# 作成した関数を実行：ハニング窓関数をかける

t_array_org, acf_org = hanning(t_array_org, Fs, N_ave_org)

t_array_filt, acf_filt = hanning(t_array_filt, Fs, N_ave_filt)

# 作成した関数を実行：FFTをかける

fft_array_org, fft_mean_org, fft_axis_org = fft_ave(t_array_org, samplerate, Fs, N_ave_org, acf_org)

fft_array_filt, fft_mean_filt, fft_axis_filt = fft_ave(t_array_filt, samplerate, Fs, N_ave_filt, acf_filt)

# dB変換する

fft_mean_org = db(fft_mean_org, 2e-5)

fft_mean_filt = db(fft_mean_filt, 2e-5)

# -----------------------------------------------------------------------------------

# ここからグラフ描画---------------------------------------------------------------------

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(211)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

ax2 = fig.add_subplot(212)

ax2.yaxis.set_ticks_position('both')

ax2.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('Time [s]')

ax1.set_ylabel('SP[Pa]')

ax2.set_xlabel('Frequency [Hz]')

ax2.set_ylabel('SPL[dB]')

# データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。

ax1.plot(x, data, label='original', lw=1)

ax1.plot(x, data_filt, label='filtered', lw=1)

ax2.plot(fft_axis_org, fft_mean_org, label='original', lw=1)

ax2.plot(fft_axis_filt, fft_mean_filt, label='filtered', lw=1)

plt.legend()

# 軸のリミットを設定する。

ax2.set_xlim(0, max(fft_axis_org)/2)

ax2.set_ylim(-20, 100)

# レイアウト設定

fig.tight_layout()

# グラフを表示する。

plt.show()

plt.close()

# -----------------------------------------------------------------------------------

実行結果

下図上段が時間波形、下段が上段の波形をFFTした後の周波数波形です。それぞれフィルタ前（original）とフィルタ後（filtered）を比較しています。

ここではFFTやデシベル変換も使っていますが、それぞれの詳細は以下の記事にまとめていますので、是非ご覧ください。

PythonでFFT！SciPyのFFTまとめ
 Pythonで音圧のデシベル(dB)変換式と逆変換式

元の時間波形はガウシアンノイズなので周波数全体でフラットな特性を持っていますが、フィルタ後は設定した周波数(fp=3000[Hz]）から綺麗に減衰していく結果を得ました。

ノイズ波形だけではつまらないので、サイン波のノイズ除去を行ってみましょう。上のメインコード（main_filter.py）の波形生成部分とフィルタ周波数を以下のように変更してみます。

samplerate = 25600
x = np.arange(0, 12800) / samplerate  # 波形生成のための間軸の作成
data = np.sin(2.0 * np.pi * 50 * x) + 0.3 * np.random.normal(loc=0, scale=1, size=len(x)) #サイン波にガウシアンノイズを重畳

fp = 200       #通過域端周波数[Hz]
fs = 400       #阻止域端周波数[Hz]
gpass = 3      #通過域端最大損失[dB]
gstop = 40     #阻止域端最小損失[dB]

samplerate = 25600

x = np.arange(0, 12800) / samplerate # 波形生成のための間軸の作成

data = np.sin(2.0 * np.pi * 50 * x) + 0.3 * np.random.normal(loc=0, scale=1, size=len(x)) #サイン波にガウシアンノイズを重畳

fp = 200 #通過域端周波数[Hz]

fs = 400 #阻止域端周波数[Hz]

gpass = 3 #通過域端最大損失[dB]

gstop = 40 #阻止域端最小損失[dB]

すると元の時間波形はノイジーなサイン波、フィルタ後はノイズが軽減している結果を得ます。このページでは極端な例を紹介していますが、これがフィルタ処理の効果です。

ソースコードの周波数や減衰量を変えると、減衰が始まる位置、遷移域の傾きを自在に変化させることができますので、是非トライしてみて下さい。

まとめ

フィルタ処理の一種であるローパスフィルタをPythonで使えるようになりました。

ローパスフィルタは周波数軸で通過域、遷移域、阻止域における周波数と減衰量を決める必要があります。

フィルタをかける時に使用する伝達関数の入力値はSciPyの関数を使うことで簡単に求めることができました。

理解するのに時間はかかったけど、ようやくフィルタ処理ができるようになってきたぞ！ここから各種フィルタに発展させていこう！
また、英語の公式ページを見ると理解が深まることがわかったので、これからは積極的に海外Webページも活用していこう！

Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

PythonのSciPyでローパスフィルタをかける！

ローパスフィルタについて

ローパスフィルタとは？

ローパスフィルタで重要な因子と用語（図解）

PythonのSciPyローパスフィルタ

ソースコード

メインファイル

関数ファイル

参考：フーリエ変換まで含めた全コード

実行結果

まとめ

コメント

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

ローパスフィルタについて

ローパスフィルタとは？

ローパスフィルタで重要な因子と用語（図解）

PythonのSciPyローパスフィルタ

ソースコード

メインファイル

関数ファイル

参考：フーリエ変換まで含めた全コード

実行結果

まとめ

SNSでもご購読できます。

コメント

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル