【Python実践】ウェーブレット変換で実用「異常検知」～閾値設定から誤検知対策まで～」

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pywt
import os
plt.rcParams['font.family'] = 'MS Gothic'

# ---------------------------------------------------------
# 1. テストデータ作成 & CSV保存関数
# ---------------------------------------------------------
def create_test_data_and_save_csv(filename='sensor_data.csv'):
    """
    異常検知テスト用の疑似データを作成し、CSVとして保存する。
    """
    np.random.seed(42)
    t = np.linspace(0, 1, 400)
    
    # 正常な振動（低周波） + ノイズ
    normal_signal = np.sin(2 * np.pi * 5 * t)
    noise = np.random.normal(0, 0.3, len(t))
    data = normal_signal + noise

    # ★異常を混入（0.5秒地点 index=200 に鋭いスパイク）
    data[200] += 4.0 

    # DataFrame化
    df = pd.DataFrame({
        'time': t,
        'value': data
    })
    
    # CSV保存 (indexは保存しない)
    df.to_csv(filename, index=False)
    print(f"[Info] テストデータを保存しました: {filename}")


# ---------------------------------------------------------
# 2. 異常検知関数 (DataFrame対応版)
# ---------------------------------------------------------
def detect_anomalies(df, target_column, wavelet='mexh', scale=1, sigma_n=3):
    """
    DataFrameの指定カラムに対してウェーブレット変換を行い、異常を検知する。
    
    Args:
        df (pd.DataFrame): 解析対象のデータフレーム
        target_column (str): 解析するカラム名
        wavelet (str): マザーウェーブレットの種類 (default: 'mexh')
        scale (int): 抽出するスケール (default: 1, 高周波・スパイク狙い)
        sigma_n (float): 閾値を決める際の標準偏差の倍率 (default: 3)
        
    Returns:
        dict: 解析結果が入った辞書
            - 'energy': ウェーブレット係数のエネルギー(絶対値)
            - 'threshold': 計算された閾値
            - 'anomaly_indices': 閾値を超えたインデックスのリスト
    """
    # データ取り出し
    data = df[target_column].values
    
    # (A) ウェーブレット変換
    # scalesはリストで渡す必要がある
    coefs, freqs = pywt.cwt(data, [scale], wavelet)
    
    # (B) ターゲット抽出 & エネルギー変換
    # coefsは [スケール数, データ長] の形状なので、index 0 を取得
    target_coef = coefs[0]
    energy = np.abs(target_coef)
    
    # (C) 統計的閾値の決定（3シグマ法）
    mean = np.mean(energy)
    std = np.std(energy)
    threshold = mean + sigma_n * std
    
    # (D) 判定 (閾値を超えたインデックスを取得)
    anomaly_indices = np.where(energy > threshold)[0]
    
    return {
        'energy': energy,
        'threshold': threshold,
        'anomaly_indices': anomaly_indices
    }


# ---------------------------------------------------------
# 3. グラフ描画関数
# ---------------------------------------------------------
def plot_results(df, target_column, result_dict):
    """
    元データと異常検知結果をグラフに描画する。
    """
    data = df[target_column].values
    energy = result_dict['energy']
    threshold = result_dict['threshold']
    anomalies = result_dict['anomaly_indices']
    
    fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

    # 上段：生データ
    axes[0].plot(data, color='gray', label='観測データ (Raw)')
    axes[0].set_title(f'元データ: {target_column} (ノイズにより異常が見えにくい)')
    axes[0].axhline(3.0, color='green', linestyle='--', label='固定閾値の例 (3.0)')
    axes[0].legend()
    axes[0].grid(True)

    # 下段：解析結果
    axes[1].plot(energy, color='blue', label='ウェーブレットエネルギー (解析値)')
    axes[1].set_title('解析結果: ノイズが除去され、異常だけが反応')
    axes[1].axhline(threshold, color='red', linestyle='--', label=f'動的閾値 (Mean + 3σ): {threshold:.2f}')
    
    # 異常箇所のプロット
    if len(anomalies) > 0:
        axes[1].scatter(anomalies, energy[anomalies], color='red', s=100, zorder=5, label='検知箇所')
        print(f"[Result] 異常検知インデックス: {anomalies}")
    else:
        print("[Result] 異常は検知されませんでした。")

    axes[1].legend()
    axes[1].grid(True)
    plt.tight_layout()
    plt.show()


# =========================================================
# メイン実行ブロック
# =========================================================
if __name__ == "__main__":
    CSV_FILE = 'factory_sensor_log.csv'
    TARGET_COL = 'value'

    # 1. データの準備（CSV作成）
    create_test_data_and_save_csv(CSV_FILE)

    # 2. データの読み込み（Pandas）
    if os.path.exists(CSV_FILE):
        df_sensor = pd.read_csv(CSV_FILE)
        print(f"[Info] データを読み込みました。行数: {len(df_sensor)}")
        
        # 3. 異常検知の実行
        # 引数にDataFrameとカラム名を渡すだけでOK
        results = detect_anomalies(df_sensor, TARGET_COL)
        
        # 4. 結果の可視化
        plot_results(df_sensor, TARGET_COL, results)
        
    else:
        print("CSVファイルが見つかりません。")