【Python実践】FFTが使えなくても大丈夫！低周期データから故障の予兆を見つけるRolling Window手法とは

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import entropy
from matplotlib import rcParams
rcParams['font.family'] = 'Meiryo'

# ---------------------------------------------------------
# 1. 共通のダミーデータ生成関数（3つの異常を仕込む）
# ---------------------------------------------------------
def generate_dummy_data(n=800):
    np.random.seed(42)
    t = np.linspace(0, 100, n)
    
    # ベースライン: 正常な稼働波形 (サイン波 + 軽いノイズ)
    base_wave = np.sin(t)
    noise = np.random.normal(0, 0.2, n)
    data = base_wave + noise

    # 【異常①】突発的なショック (t=200, 250付近)
    data[200] += 4.0
    data[250] -= 3.5

    # 【異常②】部品の劣化によるガタつき増大 (t=400〜550)
    data[400:550] += np.random.normal(0, 0.8, 150)

    # 【異常③】パッキン劣化等による微小なドリフト (t=600以降)
    data[600:] += 0.5  # ベースラインが0.5だけ浮き上がる

    return pd.DataFrame({'sensor_val': data})

# ---------------------------------------------------------
# 2. 共通のグラフ描画関数（生データと特徴量を上下に並べる）
# ---------------------------------------------------------
def plot_feature(df, feature_col, title, color='blue'):
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5), sharex=True)

    # 上段：生データ
    ax1.plot(df['sensor_val'], color='gray', alpha=0.7)
    ax1.set_title("Raw Data (生データ：3つの異常が潜んでいます)")
    
    # 異常区間のハイライト（読者に分かりやすくするため）
    ax1.axvspan(190, 260, color='red', alpha=0.1, label='①突発ショック')
    ax1.axvspan(400, 550, color='orange', alpha=0.1, label='②ガタつき増大')
    ax1.axvspan(600, 800, color='green', alpha=0.1, label='③微小ドリフト')
    ax1.legend(loc='upper left')

    # 下段：計算した特徴量
    ax2.plot(df[feature_col], color=color, linewidth=2)
    ax2.set_title(title)

    plt.tight_layout()
    plt.show()

# ▼▼▼▼▼▼▼▼▼===================▼▼▼▼▼▼▼▼▼==================▼▼▼▼▼▼▼▼▼
# ここに、以降に紹介するサンプルソースコードを貼り付けて実行してください。

# 特徴量の計算関数
def calc_rolling_95pct(df, window_size=30):
    df_calc = df.copy()
    # ノイズを無視して実質的な最大値（高負荷側への底上げ）を捉える
    df_calc['rolling_95pct'] = df_calc['sensor_val'].rolling(window=window_size).quantile(0.95)
    return df_calc

# 実行とグラフ化
df = generate_dummy_data()
df_result = calc_rolling_95pct(df)
plot_feature(df_result, 'rolling_95pct', "移動95%値 (③微小ドリフトによるベース上昇に反応)", color='red')

# 特徴量の計算関数
def calc_rolling_std(df, window_size=30):
    df_calc = df.copy()
    # 波形の「暴れ具合」を捉える
    df_calc['rolling_std'] = df_calc['sensor_val'].rolling(window=window_size).std()
    return df_calc

# 実行とグラフ化
df = generate_dummy_data()
df_result = calc_rolling_std(df)
plot_feature(df_result, 'rolling_std', "移動標準偏差 (②ガタつき増大の区間だけが綺麗に浮き上がる)", color='blue')

# 特徴量の計算関数
def calc_crest_factor(df, window_size=30):
    df_calc = df.copy()
    
    # 最大絶対値 ÷ 実効値(RMS)
    def _crest(x):
        rms = np.sqrt(np.mean(x**2))
        return 0 if rms == 0 else np.max(np.abs(x)) / rms

    df_calc['crest_factor'] = df_calc['sensor_val'].rolling(window=window_size).apply(_crest, raw=True)
    return df_calc

# 実行とグラフ化
df = generate_dummy_data()
df_result = calc_crest_factor(df)
plot_feature(df_result, 'crest_factor', "波高率 (①突発的なショックにだけ鋭く反応)", color='orange')

# 特徴量の計算関数
def calc_rolling_entropy(df, window_size=30):
    df_calc = df.copy()
    
    # ヒストグラムからシャノンエントロピーを計算
    def _entropy(x, bins=10):
        counts, _ = np.histogram(x, bins=bins)
        return entropy(counts)

    df_calc['rolling_entropy'] = df_calc['sensor_val'].rolling(window=window_size).apply(_entropy, raw=True)
    return df_calc

# 実行とグラフ化
df = generate_dummy_data()
df_result = calc_rolling_entropy(df)
plot_feature(df_result, 'rolling_entropy', "移動エントロピー (②ガタつきによる波形の乱れを検知)", color='purple')

# 特徴量の計算関数
def calc_cusum(df):
    df_calc = df.copy()
    
    # 正常時(最初の100データ)の平均値を基準に、ズレを累積していく
    baseline_mean = df_calc['sensor_val'][:100].mean()
    df_calc['cusum'] = (df_calc['sensor_val'] - baseline_mean).cumsum()
    return df_calc

# 実行とグラフ化
df = generate_dummy_data()
df_result = calc_cusum(df)
plot_feature(df_result, 'cusum', "CUSUM (③微小ドリフトが蓄積し、右肩上がりのトレンドになる)", color='green')