ただPythonでcsvからデジタルフィルタをかけるだけのコード

  • このエントリーをはてなブックマークに追加
Advertisements

今すぐ、何も考えず、とにかくcsvに記録したデータに対しデジタルフィルタをかけたい人向け。ここではPythonを知らない人のための導入を説明してから、デモcsvファイルとコピペ動作するフィルタ処理コードを紹介して目的を最速で達成します。

こんにちは。wat(@watlablog)です。ただだけシリーズ、ここではPythonを知らなくてもとにかくデジタルフィルタをかける事ができるようになる方法を紹介します

この記事の対象者

とにかくフィルタ処理をしたい人

この記事は「理論は後で良い!今はとにかくローパスフィルタやハイパスフィルタをかけなきゃならんのだ!」という人に向けて書いています。

このノイズまみれの信号を今すぐどうにかキレイにしたいけど、プログラミングの学習時間なんてない!

ちょっとcsvデータにフィルタをかけたいだけなのに、社内の高級ソフトをいちいち使うのがダルい…!

Advertisements

…という人、結構いらっしゃると思います。

Advertisements

日々実験業務を担当されている方でも、じっくり信号処理プログラムを書いている時間はほとんど無いのではと思います。

さらに、ちょっと処理したいだけなのに信号処理機能をフルに積んだ商用ソフトを使っている人もいるのではないでしょうか(計測ソフトに多いかも)。商用ソフトは社内のエンジニア同士でライセンスを予約し合って使っている場合が多いと思いますが、ちょっとした処理でライセンス待ちなんて生産性ガタ落ちです。

この記事はそんな人に向けて、比較的ハードルの低いプログラミング言語であるPythonを使ったフィルタ処理の方法を紹介します。

Pythonワカランな人

本記事ではデジタルフィルタ処理としてローパスフィルタハイパスフィルタバンドパスフィルタバンドストップフィルタPythonを使ってかけます。

しかし、Pythonの事を何も知らない人でも最後まで読み進められるように記事を構成してみました。

さらに、会社等でプロキシ設定に阻まれてライブラリインストール出来ない人も対象にしています。インターネットに接続できて、PyPIにアクセスできれば問題ありません。
もちろん全て無料です。

本記事の目標

csvファイルの複数信号を一度にフィルタ処理する

この記事は以下のフォーマットで時間波形が記録されたデータにフィルタをかけます。おそらく色々なデータロガーでcsv出力するとこのような形式になっている事でしょう。

この形式は「ただPythonでcsvから離散フーリエ変換をするだけのコード」と全く同じフォーマットであるため、フィルタをかけたりフーリエ変換したりと時間波形処理を行き来する事が出来ます。

csvファイル

1行目はヘッダです。A列に時間[s]、B列以降は各信号の名称でも書いておきます(わかりやすくするためであって、名前は何でも良いです)。

A列はフィルタ処理する分だけの時間軸を用意しておいて下さい。時間刻みは一定(等ピッチ)である必要があります。但し、フィルタをかける時の周波数が表現できていないとプログラムエラーとなりますので、ご注意下さい。
(例:0.1[s]刻みの粗いデータに1000[Hz]のフィルタをかける…等)

B列以降はA列の各時刻に対応した振幅成分(例えば電圧、加速度…といった物理的な波形)を用意します。ファイルが許す限り列方向に信号を並べておいて構いません

もっと詳しいフィルタ処理の記事を読みたい人は…

本記事は最速で、この記事だけでフィルタ処理をかける事を目標としていますが、その他過去WATLABブログで書いたフィルタ処理の記事を見たい方は以下のリンクにアクセスしてみて下さい。

ローパスフィルタ(LPF)

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

ハイパスフィルタ(HPF)

PythonのSciPyでハイパスフィルタをかける!

バンドパスフィルタ(BPF)

PythonのSciPyでバンドパスフィルタをかける!

バンドストップフィルタ(BSF)

PythonのSciPyでバンドストップフィルタをかける!

検証済の動作環境

PC

僕は以下のWindows環境、Mac環境で本記事のコードを動作検証しています。Linuxやその他OSは対象としていません。

Windows OS Windows10 64bit
CPU 2.4[GHz]
メモリ 4[GB]
Mac OS macOS Catalina 10.15.7
CPU 1.4[GHz]
メモリ 8[GB]

Python環境

この後説明するPython環境に関するバージョン情報は以下表に示す通りです。おそらく最新バージョンでも動くと思いますが、検証したのは下の環境のみ。とにかくはやくフィルタ処理したい場合は揃えておくのが無難かと思います。

PythonはPython本体、PyCharmはプログラムを記述して実行したりデバッグしたりする統合開発環境IDE)、NumpyScipyPandasmatplotlibはPythonにインポートして使う便利な外部ライブラリです。

Python Python 3.7.7
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.19.0
Scipy 1.4.1
Pandas 1.0.5
matplotlib 3.2.2

以上の前置きを確認したら、早速環境構築をしていきましょう!環境が既に構築されている人はコード部分までスクロールして下さい。

環境構築

Pythonのインストール

はじめにプログラミング言語であるPythonをインストールしましょう。

PythonのインストールにはAnacondaを推奨する書籍やサイトが沢山ありますが、2021年現在Anacondaは商用利用に制限がかかっているようです。それ以外にも色々面倒な管理となりそうであるため、筆者はAnacondaを使っていません(いちいちライブラリをインストールするのは面倒ですが)。

インストールの方法はWindowsとMacで以下の記事をご確認下さい。

Windows版:「Pythonのインストール方法とAnacondaを使わない3つの理由

Mac版:「macOSにPython3をインストールする方法をまとめてみた

そのうちもっと良い環境構築方法も試してみたいと思います(Dockerとか?)

統合開発環境(IDE)のインストール

コードを打ち込んでプログラムを実行するだけならテキストエディタを使ってコマンドプロンプトやターミナルで実行する方法でも十分ですが、デバッグやコード記述補助機能を利用するためには統合開発環境IDE)を使うのが良いです。

こちらも以下のWindowsとMacで記事を用意していますので、参照しながらインストールしてみて下さい。

Windows版:「Pythonの統合開発環境(IDE)はPyCharmで良い?

Mac版:「macOSにPython統合開発環境PyCharmをインストール

外部ライブラリのpipインストール

pip概要と外部ライブラリのインストール方法

先ほど紹介したNumpyやScipyといった外部ライブラリはpipインストールするのが一般的です。

pipについての概要は「Pythonのパッケージ管理ツール pip の使い方とコマンド集」、Windows版とMac版のコマンドの違いは「独断!Pythonでpipインストールしたい外部ライブラリまとめ」をご確認下さい。

プロキシエラーになる場合

もしかするとpipインストール時にプロキシエラーが発生するかも知れません。

PyPIにアクセスできるインターネット環境をお持ちであれば、「pip installでプロキシエラーになる問題を解決する方法」に記載の方法で解決可能です。

※もし社内プロキシ等でひっかかる人は念のためネットワーク管理者にお問い合わせした方が良いかもしれませんが。

csvをフィルタ処理するPythonコード

準備するcsvファイル【ダウンロード可】

ここからはいよいよコードを使ってフィルタ処理をしてみます。
まずはサンプルのcsvファイルとして以下の「signals.csv」をダウンロードしてみて下さい。

サンプルは10[Hz], 20[Hz], 30[Hz]のサイン波が0.001[s]の時間刻みで記録されています。

あとはこのファイルの中身を自分のデータに書き換えて下のコードを実行するだけで目的は達成できるはずです。

ただ、書き換える時はエンコードを「SHIFT-JIS」にする事を忘れずに。もし「UTF-8」で作ってもコードの方を変更すれば大丈夫ですが。

全コード【コピペ・2次利用可】

以下にcsvをフィルタ処理するだけの全コードを示します。このコードを実行するとfilter.csvというファイルが生成されます。

 

実行結果

フィルタ後のcsvファイル

生成されたcsvファイルの例を以下に示します。今回はB列に時間(signal.csvのコピー)、以降は対応する振幅のデータが最初に指定したデータ数分順番に並びます。

生成されるcsvファイル

フィルタ前後の時間波形比較

グラフの例は下図です。パッと確認したい時はPython上で見るのが一番ですね。

サンプルのプログラムはcsv_filter関数実行時にtype='lp'とローパスフィルタを指定しています。

また、関数内で通過域端周波数fp_lp=15[Hz]、阻止域端周波数fs_lp=30[Hz]を設定しているため、10[Hz]のサイン波はあまりフィルタの影響を受けませんが、20[Hz]と30[Hz]のサイン波は振幅が大きく減少している結果を得る事を出来ます。

※上段がフィルタ前、下段がフィルタ後です。

フィルタ処理前後の時間波形比較

以上でcsvファイルに記録した時間波形へフィルタ処理をかける事ができました。

しかし、csvに記録されたフィルタ後の波形を周波数軸で確認するためには、出来上がったフィルタ後のcsvファイルに対し、フーリエ変換のコードを適用させる必要があります。

…ここまでやったんだから、ねぇ?

フーリエ変換プロット確認コードも付けますかね!

csvをフィルタ処理するPythonコード(フーリエ変換機能付き)

フィルタ処理は一度設定が確定するまで、フーリエ変換で所望の結果が得られるかどうかを確認する事をよくやります。

ここではフィルタの設定をその場で確かめるためのフーリエ変換機能を追加したコードを紹介します。

…と言っても「ただPythonでcsvから離散フーリエ変換をするだけのコード」の内容と組み合わせただけで特に新しい事は何もありません!

全コード【コピペ・2次利用可】

先ほどのコードに比べ、importでfftpackをインポートしている点、「# フーリエ変換確認用------」と書いてある部分2箇所と、プロット部分を変更しています。

実行結果

ローパスフィルタ後の周波数波形確認

以下はtype='lp'で関数実行した結果です。

右側のブロックにフーリエ変換した波形をプロットしていますが、10[Hz]のピークはほぼ原型を留めているのに対し、その他の次数は振幅低減している事が周波数波形からも確かめられました。想定通りです。

ローパスフィルタ後の周波数波形

ハイパスフィルタ後の周波数波形確認

以下はtype='hp'で関数実行した結果です。

今度は高周波側である30[Hz]の次数を残し、その他の次数を低減させました。想定通りですね。

ハイパスフィルタ後の周波数波形

バンドパスフィルタ後の周波数波形確認

以下はtype='bp'で関数実行した結果です。

バンドパスの場合はデフォルトで20[Hz]が残るようにしてあります。想定通り。

バンドパスフィルタ後の周波数波形

バンドストップフィルタ後の周波数波形確認

以下はtype='bs'で関数実行した結果です。

バンドストップは逆に20[Hz]のみを低減する設定にしています。これも想定通り。

バンドストップフィルタ後の周波数波形

以上でcsvファイルにフィルタをかけるPythonコードの紹介は終了です。関数内の周波数設定を色々と変更して遊んでみて下さい!

まとめ

本ページでは検索から初めて当ブログに辿り付いた「Pythonはよくワカランけど、とにかく最速でフィルタ処理をしたい人」を対象に目標設定、Python環境の導入から説明しました。

csvファイルもサンプルをダウンロード可能としたため、環境さえ整えばすぐにフィルタ処理を試す事ができると思います。

ただ、現在のコードは周波数設定部分がcsv_filter関数の中にあるので、もしかしたらさらなる改善として関数の外から設定するようにした方が良いかも知れません(やってみて下さい!)。

また、実用性を考えフーリエ変換コードと組み合わせたコードも紹介しました。
是非自身のデータに対して色々なフィルタをかける信号処理ライフをお楽しみ下さい!

ただだけシリーズ第2段としてcsvファイルにフィルタをかけるだけのコードを書いてみました!もしただだけ記事のリクエストがありましたらコメント下さい!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメント

  1. hiroko より:

    読みやすいブログでいつも楽しみにしています。
    掲載されていたサンプルデータで試してみると、上手くいきました。ありがとうございます。
    気になったのですが、”dt”はプログラムで読み取られる設定ですが、以下のデータを使う時には”256″に置き替えた方がよいのでしょうか。

    CSVは以下のように記載しています。
    A列は0, 0.0.00390625, 0.0078125, ・・・, 511.9960938
    B列以降の各列の情報
     サンプリング数 :131,072
     サンプリング周期:1/256秒

    1. hiroko より:

      “dt”に256を入れるとエラーになってしまいました。
      なので、そのままでよかったですね。すみません。。。色々試してみたら混乱してました。

      1. wat より:

        いつもご利用ありがとうございます。
        返信遅れました!
        dtはdelta t(time)の略でよく信号処理で使う「時間刻み」、あるいは「時間分解能」と呼ばれるものです。
        データの時間刻みなので0[s], 0.1[s], 0.2[s]…とデータが記録されている場合はdt=0.1[s]です。
        hirokoさんの場合ですと、dt=0.00390625[s]でしょうか。
        プログラムでは自動で読み取っていますが、print(dt)等でご確認頂けるとわかりやすいかも知れません。
        以上、よろしくお願い致します。

コメントを残す

*