ついにWATLABブログから書籍「いきなりプログラミングPython」が発売しました！
特徴量抽出に役立つバンドパスフィルタと逆の作用をするバンドストップフィルタも習得していきます。フィルタには実際難しい理論がありますが、PythonのSciPyなら簡単に使うことが可能です。
こんにちは。wat(@watlablog)です。
基本的なフィルタを学んでいますが、今回はバンドストップフィルタを扱います！
バンドストップフィルタについて
バンドストップフィルタとは？
前回の記事「PythonのSciPyでバンドパスフィルタをかける！」では、バンドパスフィルタを使えるようになりました。
今回扱うのはバンドストップフィルタ（band-stop filter）といって、バンドパスの逆の作用をするフィルタを説明します。
バンドストップフィルタとは、ある周波数幅（バンド）のみの信号を阻止（ストップ）するフィルタのことです。
このバンドストップフィルタは音楽の分野で日常的に活用されています。電子楽器は当然電気を用いて駆動していますが、演奏で出したい音の他にも電源ノイズ（50Hz、60Hzの電源周波数に準ずる）が演奏に交じって発生するハム音という（ブーン音）ものがあります。
この音はサイン波に近い形状をしているので、倍音（高調波成分）があまり無く、このバンドストップフィルタで阻止して演奏の質を上げる用途で使われています。
その他、音声分析の分野では、ある音を消した時にどう聞こえるかといった用途で使われることもあり、使い道は様々あるフィルタです。
バンドストップフィルタで重要な因子と用語（図解）
バンドストップフィルタはある周波数幅の帯域だけを阻止させるため、以下の図のように1つの阻止域がありますが、通過域と遷移域が2つあるのが特徴です。
そのため、通過域端と阻止域端が2つずつあることになり、これまでのローパスフィルタとハイパスフィルタと比べるとプログラムする上で設定項目が増えることになります。これはバンドパスフィルタと同じ理屈です。
これまで同様、バターワースフィルタを例にコードの解説をしていきます。
PythonのSciPyバンドストップフィルタ
ソースコード
メインファイル
ソースコードを以下に示します。フィルタ処理を実際に行う部分はdef文を使って関数化しました。関数の使い方については「Pythonの関数def文の使い方！引数や別ファイル式も解説」の記事に詳細を示しています。
まずはメイン(main_filter.py)のコードです。
サンプルとしてガウシアンノイズを生成していますが、こちらについては「Pythonでガウス分布を持つノイズの作り方と調整方法」の記事を参照下さい。
中身はバンドパスフィルタの場合とほぼ同じです。詳しくは「PythonのSciPyでバンドパスフィルタをかける！」を参照下さい。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
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 = np.array([500, 6000]) #通過域端周波数[Hz]※ベクトル fs = np.array([1000, 3000]) #阻止域端周波数[Hz]※ベクトル gpass = 3 #通過域端最大損失[dB] gstop = 40 #阻止域端最小損失[dB] data_filt = filter_function.bandstop(data, samplerate, fp, fs, gpass, gstop) |
関数ファイル
関数ファイル（filter_function.py)もローパスやハイパスとほぼ同じです。signal.butterでbandstopを指定する所が変化点です。
|
1 2 3 4 5 6 7 8 9 |
#バターワースフィルタ（バンドストップ） def bandstop(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, "bandstop") #フィルタ伝達関数の分子と分母を計算 y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける return y #フィルタ後の信号を返す |
参考：フーリエ変換まで含めた全コード
先ほどまではフィルタ処理の部分のみを解説していましたが、ファイルを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 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 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 |
import numpy as np from scipy import signal from scipy import fftpack from matplotlib import pyplot as plt #バターワースフィルタ（バンドストップ） def bandstop(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, "bandstop") #フィルタ伝達関数の分子と分母を計算 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 = np.array([500, 6000]) # 通過域端周波数[Hz]※ベクトル fs = np.array([1000, 3000]) # 阻止域端周波数[Hz]※ベクトル gpass = 3 # 通過域端最大損失[dB] gstop = 40 # 阻止域端最小損失[dB] data_filt = bandstop(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)変換式と逆変換式
これまでと同じように、ノイズ波形だけではつまらないので、サイン波を二つ合成させた波形を作って効果を確認していきます。
上のメインコード（main_filter.py）の波形生成部分とフィルタ周波数を以下のように変更してみます。
|
1 2 3 4 5 6 |
data = np.sin(2.0 * np.pi * 200 * x) + np.sin(2.0 * np.pi * 1000 * x) + 0.5 * np.random.normal(loc=0, scale=1, size=len(x)) #2つのサイン波とノイズ波を合成 fp = np.array([500, 1500]) #通過域端周波数[Hz]※ベクトル fs = np.array([250, 3000]) #阻止域端周波数[Hz]※ベクトル gpass = 3 #通過域端最大損失[dB] gstop = 40 #阻止域端最小損失[dB] |
元波形は1000[Hz]の信号が200[Hz]でうねり成分を持ち、かつ全体的に高周波ノイズを含んでいます。フィルタ後の波形は1000[Hz]の信号のみを減衰させ、うねり成分とノイズ成分のみを評価できる状態にしたことになります。
まとめ
フィルタ処理の一種であるバンドストップフィルタをPythonで使えるようになりました。
バンドストップフィルタはバンドパスフィルタと同じく通過域端と阻止域端が2つずつあるフィルタであるため、周波数の設定も同じくベクトル形式になります。
SciPyを使えばバンドストップフィルタを容易に設計できることを学びました。
ローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタを使えるようになったぞ！信号処理初心者でもこんなに簡単に扱えるとは、Pythonすごい！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！