Python 振動センサーを用いた異常検知アルゴリズムの実装【RMS】

import numpy as np
from scipy.signal import butter, sosfiltfilt


# =============================================================================
# 定数
# =============================================================================

# 1 秒間のサンプリング回数 (Hz)
SAMPLE_RATE_HZ = 1000.0

# 継続時間 (s)
SAMPLE_DURATION_S = 2.0

# ハイパスフィルタのカットオフ周波数 (Hz)
HIGHPASS_HZ = 5.0

# 閾値 (mm/s)
THRESHOLD_MM_S = 4.5

# サンプル用の時刻データ (s)
SAMPLE_TIME_S = np.arange(
    0.0,
    SAMPLE_DURATION_S,
    1.0 / SAMPLE_RATE_HZ,
)

# サンプル用の加速度データ (m/s²)
SAMPLE_ACCELERATION_M_S2 = (
    0.5 * np.sin(2.0 * np.pi * 50.0 * SAMPLE_TIME_S)
    + 0.05 * np.sin(2.0 * np.pi * 10.0 * SAMPLE_TIME_S)
)


# =============================================================================
# RMS 速度の算出
# =============================================================================
def calculate_rms_velocity(
    time_s: np.ndarray,
    acceleration: np.ndarray,
    highpass_hz: float = HIGHPASS_HZ,
) -> float:

    time_s = np.asarray(time_s, dtype=float)
    acceleration = np.asarray(acceleration, dtype=float)

    if time_s.ndim != 1 or acceleration.ndim != 1:
        raise ValueError("入力データは 1 次元配列にしてください。")

    if len(time_s) != len(acceleration) or len(time_s) < 16:
        raise ValueError("データ数が不正です。")

    if not (
        np.all(np.isfinite(time_s))
        and np.all(np.isfinite(acceleration))
    ):
        raise ValueError("NaN または無限大が含まれています。")

    dt = np.diff(time_s)

    if np.any(dt <= 0):
        raise ValueError("時刻は単調増加である必要があります。")

    dt_mean = np.mean(dt)

    if not np.allclose(dt, dt_mean, rtol=0.02, atol=1e-12):
        raise ValueError("サンプリング間隔が不均一です。")

    fs = 1.0 / dt_mean
    nyquist = fs / 2.0

    if not 0 < highpass_hz < nyquist:
        raise ValueError("カットオフ周波数が不正です。")

    # 加速度の DC オフセットを除去する
    acceleration = acceleration - np.mean(acceleration)

    # ハイパスフィルタ
    sos = butter(4, highpass_hz, btype="highpass", fs=fs, output="sos")
    filtered_acceleration = sosfiltfilt(sos, acceleration)

    # 台形則で加速度を積分し、速度を算出する
    velocity = np.zeros_like(filtered_acceleration)
    velocity[1:] = np.cumsum(
        (filtered_acceleration[1:] + filtered_acceleration[:-1])
        * 0.5
        * dt
    )

    # 速度の平均値を除去する
    velocity -= np.mean(velocity)

    # m/s から mm/s に変換して RMS を算出する
    rms_velocity = np.sqrt(np.mean(velocity ** 2))

    return float(rms_velocity * 1000.0)


# =============================================================================
# 異常判定
# =============================================================================
def detect_anomaly(
    rms_velocity: float,
    threshold: float = THRESHOLD_MM_S,
) -> bool:

    if threshold <= 0:
        raise ValueError("閾値は正の値にしてください。")

    return rms_velocity >= threshold


# =============================================================================
# メイン処理
# =============================================================================
def main():

    # 定数からデータを取得する
    time_s = SAMPLE_TIME_S
    acceleration = SAMPLE_ACCELERATION_M_S2

    # RMS 速度の算出
    rms_velocity = calculate_rms_velocity(
        time_s,
        acceleration,
    )

    # 異常判定
    anomaly = detect_anomaly(
        rms_velocity,
        THRESHOLD_MM_S,
    )

    # 表示
    print(f"RMS 速度 : {rms_velocity:.3f} mm/s")
    print(f"閾値 : {THRESHOLD_MM_S:.3f} mm/s")

    if anomaly:
        print("判定結果 : 異常")
    else:
        print("判定結果 : 正常")


# =============================================================================
# プログラムの実行開始位置
# =============================================================================
if __name__ == "__main__":
    main()