【Python実践】軸受け（回転系）の異常振動を見つけ出せ！実用エンベローブ解析４つの活用シーン　

import numpy as np

def generate_bearing_fault_data(fs, duration, rpm, fault_freq, resonance_freq=3000):
    """
    ベアリングの異常（衝撃）シミュレーションデータを生成する
    
    Parameters:
    -----------
    fs : int          サンプリング周波数 (Hz)
    duration : float  データ長 (秒)
    rpm : float       軸の回転数 (毎分)
    fault_freq : float 損傷周波数 (Hz) - 1秒間に衝撃が起きる回数
    resonance_freq : float 機械の共鳴周波数 (Hz) - 叩いた時の「キーン」という音
    """
    # 時間軸の作成
    t = np.linspace(0, duration, int(fs * duration))
    
    # 1. 回転に伴う低周波成分（軸のうねりや振動のベース）
    fr = rpm / 60
    vibration_base = 0.5 * np.sin(2 * np.pi * fr * t)
    
    # 2. 周期的な衝撃波（損傷によるバースト信号）の生成
    impacts = np.zeros_like(t)
    # 損傷周波数に基づいた発生タイミング（秒）を算出
    impact_times = np.arange(0, duration, 1/fault_freq)
    
    for start_t in impact_times:
        idx = int(start_t * fs)
        if idx < len(t):
            # 衝撃発生時の減衰振動（共鳴音）をシミュレート
            # 1回あたりの衝撃の長さを約20ms程度とする
            length = min(int(fs * 0.02), len(t) - idx)
            time_segment = np.linspace(0, length/fs, length)
            
            # 指数関数的に減衰する高い周波数のサイン波
            impact_wave = np.sin(2 * np.pi * resonance_freq * time_segment) * np.exp(-500 * time_segment)
            impacts[idx:idx+length] += impact_wave

    # 3. 全体的な背景ノイズの加算
    np.random.seed(42) # 再現性のためにシードを固定
    noise = np.random.normal(0, 0.2, len(t))
    
    # すべてを合成して「生の振動データ」を模倣する
    raw_data = vibration_base + impacts + noise
    
    return t, raw_data

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'MS Gothic'

# --- 動作確認 ---
if __name__ == "__main__":
    FS = 10000    # 10kHz
    DUR = 0.5     # 表示用に短めの0.5秒
    RPM = 1800    # 30Hz
    BPFO = 105.5  # 外輪損傷の理論値
    
    t, data = generate_bearing_fault_data(FS, DUR, RPM, BPFO)
    
    plt.figure(figsize=(10, 4))
    plt.plot(t, data, color='gray', linewidth=0.5)
    plt.title("生成されたベアリング異常のシミュレーション波形")
    plt.xlabel("時間 [s]")
    plt.ylabel("振幅")
    plt.grid(True, alpha=0.3)
    plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
plt.rcParams['font.family'] = 'MS Gothic'

# ===================================================================================== 
# --- 1. シミュレーションデータ生成関数 ---
# ここに前述した generate_bearing_fault_data 関数を貼り付けて実行してください
# =====================================================================================

# --- 2. エンベロープ解析ロジック ---
def run_envelope_analysis(raw_data, fs, lowcut=2000, highcut=4500):
    """
    バンドパスフィルタ + ヒルベルト変換によるエンベロープ抽出
    """
    nyq = 0.5 * fs
    b, a = signal.butter(4, [lowcut/nyq, highcut/nyq], btype='band')
    filtered = signal.filtfilt(b, a, raw_data)
    
    # ヒルベルト変換で包絡線を抽出
    envelope = np.abs(signal.hilbert(filtered))
    
    # エンベロープFFT用のスペクトル計算
    envelope_detrend = envelope - np.mean(envelope)
    n = len(raw_data)
    f_axis = np.fft.fftfreq(n, d=1/fs)
    spec = np.abs(np.fft.fft(envelope_detrend)) / (n / 2)
    
    return envelope, f_axis[:n//2], spec[:n//2]

# --- 3. グラフ可視化関数 ---
def plot_diagnostic_results(f_axis, spec, theoretical_freq, title="エンベロープFFT解析結果"):
    """
    解析結果をプロットし、理論値（BPFO等）と比較する
    """
    plt.figure(figsize=(10, 5))
    plt.plot(f_axis, spec, color='red', label='エンベロープスペクトル')
    
    # 理論値の補助線
    plt.axvline(theoretical_freq, color='blue', linestyle='--', alpha=0.7, 
                label=f'理論値: {theoretical_freq}Hz')
    
    plt.title(title)
    plt.xlabel("周波数 [Hz]（＝衝撃のリズム）")
    plt.ylabel("振幅")
    plt.xlim(0, 500) # 診断に重要な低周波域を拡大
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

# --- メイン処理 ---
if __name__ == "__main__":
    # パラメータ設定
    FS = 10000
    DUR = 1.0
    RPM = 1800
    BPFO_TARGET = 105.5 # 解析の標的となる外輪損傷周波数
    
    # Step 1: データ生成
    t, raw_sig = generate_bearing_fault_data(FS, DUR, RPM, BPFO_TARGET)
    
    # Step 2: 解析実行
    env, f, s = run_envelope_analysis(raw_sig, FS, lowcut=2000, highcut=4500)
    
    # Step 3: 結果表示
    plot_diagnostic_results(f, s, BPFO_TARGET)

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
plt.rcParams['font.family'] = 'MS Gothic'

# ===================================================================================== 
# --- 1. シミュレーションデータ生成関数 ---
# ここに前述した generate_bearing_fault_data 関数を貼り付けて実行してください
# =====================================================================================

# --- 2. エンベロープ解析ロジック ---
def run_envelope_analysis(raw_data, fs, lowcut=2000, highcut=4500):
    """
    バンドパスフィルタ + ヒルベルト変換によるエンベロープ抽出
    """
    nyq = 0.5 * fs
    b, a = signal.butter(4, [lowcut/nyq, highcut/nyq], btype='band')
    filtered = signal.filtfilt(b, a, raw_data)
    
    # ヒルベルト変換で包絡線（エンベロープ）を抽出
    envelope = np.abs(signal.hilbert(filtered))
    
    return envelope

# --- 3. グラフ可視化関数（周波数軸：特定用） ---
def plot_envelope_spectrum(envelope, fs, theoretical_freq):
    """
    エンベロープFFTを実行し、理論値（基本波・高調波）と照合する
    """
    n = len(envelope)
    
    # 直流成分（平均値）を引いてからFFTをかけるのがスペクトルを綺麗にするコツ
    envelope_detrend = envelope - np.mean(envelope)
    
    # FFT実行
    f_axis = np.fft.fftfreq(n, d=1/fs)
    spec = np.abs(np.fft.fft(envelope_detrend)) / (n / 2)
    
    # 可視化（正の周波数領域のみ。診断に必要な低周波域 0~500Hz に限定）
    mask = (f_axis >= 0) & (f_axis <= 500)
    
    plt.figure(figsize=(10, 5))
    plt.plot(f_axis[mask], spec[mask], color='red', label='エンベロープスペクトル')
    
    # 理論値（基本波）に補助線を引く
    plt.axvline(theoretical_freq, color='blue', linestyle='--', alpha=0.7, 
                label=f'理論値: {theoretical_freq}Hz')
    
    # 高調波（2倍, 3倍）にも線を引く（これが並んでいれば異常の確証が持てる）
    plt.axvline(theoretical_freq * 2, color='blue', linestyle=':', alpha=0.5, label='高調波 (2x, 3x)')
    plt.axvline(theoretical_freq * 3, color='blue', linestyle=':', alpha=0.5)
    
    plt.title("エンベロープFFTによる損傷部位の特定")
    plt.xlabel("周波数 [Hz]（＝衝撃のリズム）")
    plt.ylabel("振幅")
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)
    plt.show()

# --- 4. メイン処理 ---
if __name__ == "__main__":
    # パラメータ設定
    FS = 10000          # サンプリング周波数
    DUR = 1.0           # 観測時間
    RPM = 1800          # 回転数（30Hz）
    BPFO_TARGET = 105.5 # 解析の標的となる外輪損傷周波数の理論値
    
    # Step 1: シミュレーションデータの生成
    # ※ generate_bearing_fault_data が定義されている必要があります
    t, raw_sig = generate_bearing_fault_data(FS, DUR, RPM, BPFO_TARGET)
    
    # Step 2: エンベロープ（包絡線）の抽出
    env = run_envelope_analysis(raw_sig, FS)
    
    # Step 3: 周波数スペクトルの表示と理論値照合
    plot_envelope_spectrum(env, FS, BPFO_TARGET)

import numpy as np
from scipy.stats import kurtosis

# ===================================================================================== 
# --- 1. シミュレーションデータ生成関数 ---
# ここに前述した generate_bearing_fault_data 関数を貼り付けて実行してください
# ===================================================================================== 

# ===================================================================================== 
# --- 2. エンベロープ解析ロジック ---
# ここに前述した run_envelope_analysis 関数を貼り付けて実行してください
# ===================================================================================== 

# --- 3. 特徴量抽出関数（数値化） ---
def extract_envelope_features(envelope):
    """
    エンベロープ波形から「異常のひどさ」を示す特徴量を抽出する
    """
    # 1. 尖度（Kurtosis）: 3に近いほど正常、大きいほど鋭い衝撃（異常）がある
    kurt_val = kurtosis(envelope)
    
    # 2. クレストファクタ: 衝撃の鋭さを測る（ピーク値 / 実効値）
    rms_val = np.sqrt(np.mean(envelope**2))
    crest_factor = np.max(envelope) / rms_val if rms_val > 0 else 0
    
    # 3. 衝撃指数（Impulse Factor）: 平均に対するピークの比率
    impulse_factor = np.max(envelope) / np.mean(envelope) if np.mean(envelope) > 0 else 0
    
    return {
        "Kurtosis": kurt_val,
        "CrestFactor": crest_factor,
        "ImpulseFactor": impulse_factor
    }

# --- 4. メイン処理 ---
if __name__ == "__main__":
    # パラメータ設定
    FS = 10000
    DUR = 1.0
    RPM = 1800
    BPFO_TARGET = 105.5
    
    # Step 1: 異常データの生成
    _, raw_sig = generate_bearing_fault_data(FS, DUR, RPM, BPFO_TARGET)
    
    # Step 2: エンベロープ抽出
    env = run_envelope_analysis(raw_sig, FS)
    
    # Step 3: 特徴量の算出
    features = extract_envelope_features(env)
    
    # Step 4: 結果の出力
    print("-" * 30)
    print("エンベロープ特徴量 解析結果")
    print("-" * 30)
    for key, val in features.items():
        print(f"{key:15}: {val:.4f}")
    print("-" * 30)
    print("※ 正常データなら Kurtosis は 3 前後になります。")

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# ===================================================================================== 
# --- 1. シミュレーションデータ生成関数 ---
# ここに前述した generate_bearing_fault_data 関数を貼り付けて実行してください
# ===================================================================================== 

# --- 2. 合わせ技解析・可視化関数 ---
def plot_advanced_combined_analysis(data, fs):
    """
    STFT（スペクトログラム）とエンベロープを並べて表示し、
    異常帯域の特定と包絡線抽出を同時に行う
    """
    # 1. STFTの計算（時間・周波数の同時解析）
    f, t_stft, Zxx = signal.stft(data, fs, nperseg=256, noverlap=128)
    
    # 2. 比較用のエンベロープ計算
    # STFTで見つけた「光る縦線（衝撃）」の帯域（2000-4500Hz）を狙い撃ちします
    nyq = 0.5 * fs
    b, a = signal.butter(4, [2000/nyq, 4500/nyq], btype='band')
    filtered = signal.filtfilt(b, a, data)
    envelope = np.abs(signal.hilbert(filtered))
    
    # 3. 可視化
    plt.figure(figsize=(12, 8))
    
    # 上段：STFTスペクトログラム
    plt.subplot(2, 1, 1)
    # パワーをdB表記にすると、微小な衝撃の光る線が見やすくなります
    plt.pcolormesh(t_stft, f, 20*np.log10(np.abs(Zxx) + 1e-6), shading='gouraud', cmap='magma')
    plt.title("STFTスペクトログラム（縦に走る光る線が衝撃の証拠）")
    plt.ylabel("周波数 [Hz]")
    plt.colorbar(label="強度 [dB]")
    
    # 下段：抽出されたエンベロープ波形
    plt.subplot(2, 1, 2)
    t_env = np.linspace(0, len(data)/fs, len(data))
    plt.plot(t_env, envelope, color='red', linewidth=1, label="抽出された包絡線")
    plt.title("特定帯域（2000-4500Hz）から抽出したエンベロープ波形")
    plt.xlabel("時間 [s]")
    plt.ylabel("振幅")
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

# --- 3. メイン処理 ---
if __name__ == "__main__":
    # パラメータ設定
    FS = 10000
    DUR = 1.0
    RPM = 1800
    BPFO_TARGET = 105.5
    
    # Step 1: データの生成
    t, raw_sig = generate_bearing_fault_data(FS, DUR, RPM, BPFO_TARGET)
    
    # Step 2: 合わせ技解析の実行
    plot_advanced_combined_analysis(raw_sig, FS)