ディジタル信号処理｜周波数解析→フィルタ設計の具体例をPythonで実装

いよいよ！ディジタル信号処理の総集編です！
周波数解析→フィルタ設計→フィルタリングの一連の流れを、サンプルプログラムで確認します。

なお、サンプルプログラムで使用しているaripac.signal_processingモジュールは、Githubに公開してありますので、ぜひチェックしてください🐜

https://github.com/ari23ant/python

フォルダ構成
開発環境
サンプルプログラム
解説｜時間遅れ
おわりに

フォルダ構成

以下のような構成を前提とします。

dev
├─sample_asp.py
└─aripac
　　├─signal_processing.py
　　├─__init__.py

開発環境

後述するサンプルプログラムを実行する開発環境は、次の通りです。

項目	内容
OS	Windows 10
Python	3.7.6(Anaconda 4.8.1)
Numpy	1.17.4
Scipy	1.3.2
matplotlib	3.1.1

サンプルプログラム

Pythonを使ったFFTのサンプルプログラム(sample_asp.py)は以下の通りです。自由にコピペして、実際に動かしてみてください。

#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import aripac.signal_processing as asp

class TestFFT():
    def __init__(self):
        self.msg = 'Test script for aripac.signal_processing'
        # 設定値
        self.fs = 250   # サンプリング周波数
        self.sec = 10   # 時間
        self.start_time = 2.25  # 抽出開始時間
        self.period = 3.5       # 抽出期間

        # 目的の信号
        self.A = 1  # 振幅
        self.f = 5  # 周波数

        # ノイズ｜sin波
        self.An = 0.1
        self.fn = 50

        # （単純）移動平均
        self.ma_num = 5
        #self.ma_num = 55

    def process(self):
        print(self.msg)
        #------- 1次元の信号波形を準備し周波数解析 -------#
        #--- 波形用意 ---#
        # 目的の信号
        t = np.arange(0, self.sec, 1/self.fs)  # 時間軸
        target_signal = self.A * np.sin(2.0 * np.pi * self.f * t)  # 抽出したい信号波形

        # ノイズ｜sin波
        noise_sin = self.An * np.sin(2.0 * np.pi * self.fn * t)

        # 合成波
        conbined_signal = target_signal + noise_sin

        # 波形の一部を抜き出す 実際と同じように両端を半端に切り取る
        start = int(self.start_time * self.fs)  # インデックスは整数
        end = int(self.start_time + self.period) * self.fs
        y = conbined_signal[start:end]
        x = t[start:end]

        #--- 周波数解析 ---#
        # FFT
        asp.fft(y, self.fs)

        #------- ノイズ除去するローパスフィルタを設計 -------#
        #--- フィルタ準備 ---#
        # 移動平均フィルタ
        lowpass = asp.make_MA_filter(self.ma_num)

        #--- 周波数応答 ---#
        asp.fvtool(lowpass, self.fs)

        #------- 信号にフィルタをかけてノイズが除去できているか確認 -------#
        #--- 畳み込み ---#
        # FIRフィルタならnumpy.convolveとscipy.signal.lfilterはおなじ畳み込み式だが、
        # 時間遅れの扱いがだいたい倍異なるので十分注意すること
        lp_y = np.convolve(y, lowpass, 'valid')
        delay = int( (len(lowpass)-1) / 2 )  # 時間遅れ（移動平均はFIRフィルタ）
        lp_x = x[delay:-delay]  # lp_yに合わせて時間遅れを調整

        #--- 周波数解析 ---#
        # FFT
        asp.fft(lp_y, self.fs)

        #------- フィルタリング前後で波形を比較 -------#
        plt.figure()
        plt.plot(x, y, label='raw')
        plt.plot(lp_x, lp_y, label='filt')
        plt.grid()
        plt.legend()

        plt.show()


if __name__ == '__main__':
    proc = TestFFT()
    proc.process()

実行結果は以下の通りです。今回は4枚のグラフがあります。

まず、生波形のFFTの結果です。5Hzと50Hzの周波数が混在していることがわかります。

設計したフィルタの応答の結果です。50Hzの周波数を除去する単純移動平均のローパスフィルタを設計できています。

単純移動平均フィルタをaripac.signal_processingでfvtooした例

ローパスフィルタをかけた波形を改めてFFTした結果です。現実はここまできれいではありませんが、狙い通りフィルタリングできていることがわかります。

フィルタリングした波形をaripac.signal_processingでFFTした例

生波形とフィルタリングした波形を重ねてプロットした結果です。フィルタリングすると、生波形にあったギザギザが取れています。

解説｜時間遅れ

すでに以前の記事で解説したので、そちらを参照してください。

ディジタル信号処理｜周波数解析とFFT（実践） - ari23の研究ノート

ディジタル信号処理｜フィルタの周波数応答をPythonで実装 - ari23の研究ノート

今回はフィルタリングのところだけ取り上げます。

        #--- 畳み込み ---#
        # FIRフィルタならnumpy.convolveとscipy.signal.lfilterはおなじ畳み込み式だが、
        # 時間遅れの扱いがだいたい倍異なるので十分注意すること
        lp_y = np.convolve(y, lowpass, 'valid')
        delay = int( (len(lowpass)-1) / 2 )  # 時間遅れ（移動平均はFIRフィルタ）
        lp_x = x[delay:-delay]  # lp_yに合わせて時間遅れを調整

ソースコード上のコメントにもありますが、時間遅れの扱いに注意します。

FIRフィルタ（移動平均フィルタ）では、以下の図のように順次畳み込みをします。

上記の図より、（フィルタのタップ数-1)/2 分だけ両端のデータがなくなります。つまり、強いFIRフィルタフィルタをかけてしまうと、その分データ数が短くなり、トレードオフの関係にあることがわかります。

もちろん、scipy.signal.lfilterを使っても構いませんが、どのように時間遅れが発生するのかは十分注意して使用することが大切です。
なぜなら、複数の信号をそれぞれフィルタリングしたとき、時間遅れをきちんと扱わないと、本来のタイムスタンプからずれてしまい、時刻が一致しないまま解析にかけてしまうおそれがあります。