どうも、こざくらです。
今日はPythonを使った業務改善メモです。

背景

回転体に取り付けたセンサーから出力電圧をデータ処理をする機会がありました。

出力電圧には、いろいろなノイズ（電源周波数とか高次成分などなど）が乗っちゃってるので、ノイズ処理が必要です。
また、機械は回転体なので、n次成分に対しても適切な処理をしないといけません。

これらの処理を汎用ソフトでこなしていた（というか人に頼んでいた）のですが、時間がかかるし、何より処理がブラックボックス化されて自分で出来ない＆理解できてない状態が不満だったので、勉強も兼ねてPythonで実装してみました。

コード内容

生データの取得

モジュールのインポートからデータの準備は以下の通り。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os

# データのパラメータ
sec = 10             # 測定秒数 
dt = 0.0001          # サンプリング間隔
N = int(sec / dt)         # サンプル数
rps = 120
# 1, 2次周波数の±20Hzさせたら、漏れなくカットオフできそう
fc1 = rps - 20  # カットオフ周波数1[Hz]
fc2 = rps * 2 + 20  # カットオフ周波数2（1・2次BPF）[Hz]
t = np.arange(0, N*dt, dt)  # 時間軸
freq = np.linspace(0, 1/dt, N)  # 周波数軸

# CSVのロード
filepath = r"C:\Users\hogehoge"
filename = os.path.basename(filepath)  # 拡張子を含むファイル名部分の文字列
dataname = os.path.splitext(os.path.basename(filepath))[0]  # 拡張子を除くファイル名部分の文字列

df = pd.read_csv(filepath, encoding="UTF-8", skiprows=9)

# 保存フォルダ作成
save_dir = os.path.splitext(os.path.basename(filepath))[0]
if not os.path.exists(save_dir):
    os.mkdir(save_dir)

ポイント

・np.arange()
引数はnp.arange(start, stop, step)となり、等差数列を配列ndarrayとして生成する。
・np.linspace()
np.arange()と同様に等差数列を生成するけど、間隔(step)ではなく、np.linespace(start, stop, num)のように要素数を指定する。
・os.path.splitext()
ファイル名やフォルダ名から拡張子(.txtとか)を取得するメソッド。
そのため、ピリオドはパス名の1番右で分割する。下記が実行例。

filename = "/python/python.txt"
ex = os.path.splitext(filename)

print(ex)
print(ex[0])
print(ex[1])

# 実行結果
('/python/python', '.txt')
/python/python
.txt

FFT分析（ノイズ処理前）

ノイズ処理前後での差を分かりやすくするため、まずはノイズ処理前のデータについてFFT分析を行います。　

#必要データの抽出
f = df['voltage']
# numpy.ndarrayに変換
f = f.values

# ここからノイズ処理前のデータ分析
# 高速フーリエ変換（周波数信号に変換）
F = np.fft.fft(f)
# FFTの複素数をを絶対値に変換
F_abs = np.abs(F)
# 振幅をもとの信号に揃える
F_abs_amp = F_abs / N * 2  # 交流成分はデータ数で割って2倍
F_abs_amp[0] = F_abs_amp[0] / 2  # 直流成分は2倍不要

ポイント

・DataFrameからndarrayへの変換
まず、csvの元データをpandasのDataFrameとして取得しました。
df = pd.read_csv(hogehoge)
その後、電圧の列を1次元のSeriesとして抜き出します。
f = df['voltage]
1次元のデータになったので、数値計算に適したnumpyのndarrayで取得します。
後でFFT分析をするということもあってndarrayにする。
f = f.values

・高速フーリエ変換
np.fft.fft(f)でFFT出来ちゃう。めっちゃ便利。後述するけど逆FFTももちろん出来てしまう。numpy恐るべし。

FFT分析（ノイズ処理後）

続いて、ノイズ処理+FFT分析のパターン。

# ノイズフィルタリング用にFFT結果コピー
F2 = np.copy(F)
#バンドパス処理（fc1～fc2の帯域以外を0にする）
F2[(freq < fc1)] = 0
F2[(freq > fc2)] = 0
# FFTの複素数をを絶対値に変換
F2_abs = np.abs(F2)
# 振幅をもとの信号に揃える
F2_abs_amp = F2_abs / N * 2  # 交流成分はデータ数で割って2倍
F2_abs_amp[0] = F2_abs_amp[0] / 2  # 直流成分は2倍不要
# 逆FFTで時間信号に戻す
f2_ifft = np.fft.ifft(F2)
# ナイキスト周波数以降も全部ゼロにしちゃったから、振幅を揃えるため
f2_ifft_real = f2_ifft.real * 2 # 実数部の取得かつ2倍して、振幅を元スケールに戻す

ポイント

・バンドパス処理（fc1～fc2の帯域以外を0にする）
1次の回転数に応じてバンドパスフィルタ(BPF)を行う周波数域を変えています。
例えば1次成分が90rpsであれば、70Hz (= 90Hz-20Hz) 以下および200Hz (= 90Hz× 2 + 20Hz)以上の電圧を0にしています。
200Hz以上が電圧0になっているので、これによって高次ノイズもキャンセルされています。

グラフ化してアウトプット

最後にmatplotlibを使ったグラフ化処理です。

# グラフ表示
fig = plt.figure(figsize=(10.0, 8.0))
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 10

# 時間信号（元）
plt.subplot(221)
plt.plot(t, f, label='f(n)')
plt.xlabel("Time", fontsize=12)
plt.ylabel("Signal", fontsize=12)
plt.grid()
leg = plt.legend(loc=1, fontsize=15)
leg.get_frame().set_alpha(1)

#周波数スペクトル（元）
plt.subplot(222)
plt.plot(freq, F_abs_amp, label='|F(k)|')
plt.xlabel('Frequency', fontsize=12)
plt.ylabel('Amplitude', fontsize=12)
plt.grid()
leg = plt.legend(loc=1, fontsize=15)
leg.get_frame().set_alpha(1)

# 時間信号（処理後）
plt.subplot(223)
plt.plot(t, f2_ifft_real, label='f2(n)')
plt.xlabel("Time", fontsize=12)
plt.ylabel("Signal", fontsize=12)
plt.grid()
leg = plt.legend(loc=1, fontsize=15)
leg.get_frame().set_alpha(1)

#周波数スペクトル（処理後）
plt.subplot(224)
plt.plot(freq, F2_abs_amp, label='|F2(k)|')
plt.xlabel('Frequency', fontsize=12)
plt.ylabel('Amplitude', fontsize=12)
plt.grid()
leg = plt.legend(loc=1, fontsize=15)
leg.get_frame().set_alpha(1)
plt.show()

# 画像保存
img_name = os.path.splitext(os.path.basename(filepath))[0]
fig.savefig(dataname + '.png')