産業用振動センサーデータを用いた高速フーリエ変換（FFT）と異常周波数検知アルゴリズムの実装

はじめに

こんにちは、NKKTech Global 技術チームです。

工場の稼働停止（ダウンタイム）を防ぐ状態基準保全（CBM: Condition-Based Maintenance）において、振動解析は「機械の脈拍を診る聴診器」と言えます。
ベアリングの内輪・外輪の微小な剥離やシャフトの芯狂い（ミスアライメント）は、時間領域の波形データではノイズに埋もれますが、**周波数領域へ変換（スペクトル解析）**することで、固有のピークとして明確に浮き彫りになります。

本記事では、生データの収集からサンプリング定理の考慮、窓関数処理、FFT変換、そして統計的動的閾値（3-Sigma）を用いた異常周波数の自動検知ロジックまでを実装コード付きで解説します。

---

1. 振動解析の基礎理論と特徴周波数

回転機械の異常は、回転数（回転基本周波数: $f_r$ [Hz]）の整数倍やベアリングの幾何学形状に起因する特徴周波数に現れます。

異常モード	主な周波数特徴	物理的要因
アンバランス (Unbalance)	$1 \times f_r$ (回転1次) が支配的	ローターの重心のズレ
ミスアライメント (Misalignment)	$2 \times f_r$ (回転2次) および高調波	軸芯の平行度・角度の狂い
ガタ・緩み (Looseness)	$0.5 \times f_r$ などの分数次調波、多数の高調波	ボルト緩み、軸受隙間過大
ベアリング外輪損傷 (BPFO)	固有周波数（回転数とボール径・ピッチ径から計算）	軌道面への傷・フレーキング

サンプリング定理（ナイキスト周波数）の遵守

検知したい最高周波数を $f_{max}$ とするとき、サンプリング周波数 $f_s$ は最低でも $2 \times f_{max}$ 以上が必要です。実務では折り返し歪み（エイリアシング）を避けるため、$f_s \ge 2.56 \times f_{max}$ を確保します（例: 2kHzまでの異常を見るなら $f_s = 5,120\text{ Hz}$ 以上）。

---

2. 信号処理パイプラインの設計

FFTを正しく適用するには、以下のパイプラインが不可欠です。

---

3. Pythonによる完全実装コード

以下のスクリプトでは、正常回転信号に「ベアリング傷を模した異常高周波パルス」が混入した合成データを生成し、FFT変換と自動ピーク検知を行います。

必要なライブラリ

pip install numpy scipy matplotlib

実装スクリプト (vibration_fft_analyzer.py)

import numpy as np
from scipy.fft import rfft, rfftfreq
from scipy.signal import find_peaks, windows
import matplotlib.pyplot as plt

# ---------------------------------------------------------
# 1. サンプル振動データの合成生成 (1800 RPM = 30 Hz のモーター)
# ---------------------------------------------------------
np.random.seed(42)
fs = 5000       # サンプリングレート: 5,000 Hz (ナイキスト: 2,500 Hz)
duration = 2.0  # 計測時間: 2秒
N = int(fs * duration)
t = np.linspace(0, duration, N, endpoint=False)

fr = 30.0  # 回転基本周波数: 30 Hz (1X)

# 正常振動成分: 1X (アンバランス) + 2X (アライメント成分) + 背景白色ノイズ
signal_normal = (
    1.2 * np.sin(2 * np.pi * fr * t) +
    0.4 * np.sin(2 * np.pi * (2 * fr) * t) +
    np.random.normal(0, 0.3, N)
)

# 異常成分: ベアリング外輪欠陥を想定した 235 Hz の異常ピーク
defect_freq = 235.0
signal_defect = 0.8 * np.sin(2 * np.pi * defect_freq * t)

# 合成波形 (DCバイアス +0.5V 加算)
raw_vibration = signal_normal + signal_defect + 0.5

# ---------------------------------------------------------
# 2. 信号前処理と高速フーリエ変換 (FFT)
# ---------------------------------------------------------
# (1) DCオフセット除去 (平均値を引いて0を中心にする)
detrended_signal = raw_vibration - np.mean(raw_vibration)

# (2) ハン窓 (Hanning Window) の適用 (スペクトル漏れの抑制)
window = windows.hann(N)
windowed_signal = detrended_signal * window

# 窓関数の振幅補正係数 (コヒーレントゲインの補正)
coherent_gain = np.sum(window) / N

# (3) 実数FFTの計算
fft_values = rfft(windowed_signal)
frequencies = rfftfreq(N, d=1/fs)

# (4) 振幅スペクトル（片側）の正規化 [単位: Acceleration (g)]
# rfftの結果はN/2で正規化し、窓関数の減衰を補正する
amplitude_spectrum = (np.abs(fft_values) / (N / 2)) / coherent_gain

# ---------------------------------------------------------
# 3. 異常周波数の自動検知アルゴリズム (動的統計閾値)
# ---------------------------------------------------------
# 移動窓を用いたローカルなノイズフロア推定 (移動平均 + 3σ)
kernel_size = 50
baseline_mean = np.convolve(amplitude_spectrum, np.ones(kernel_size)/kernel_size, mode='same')
baseline_std = np.sqrt(np.convolve((amplitude_spectrum - baseline_mean)**2, np.ones(kernel_size)/kernel_size, mode='same'))

# 動的閾値: 平均 + 5 * 標準偏差 (または固定下限値)
dynamic_threshold = np.maximum(baseline_mean + 5 * baseline_std, 0.15)

# ピーク検出 (閾値を超える極大値を抽出)
peaks, properties = find_peaks(
    amplitude_spectrum,
    height=dynamic_threshold,
    distance=int(5 / (fs / N)) # 周波数分解能に基づき最低5Hz離れたピークのみ
)

peak_freqs = frequencies[peaks]
peak_amps = amplitude_spectrum[peaks]

# ---------------------------------------------------------
# 4. 結果判定とコンソール出力
# ---------------------------------------------------------
print("=== 振動スペクトル解析結果 ===")
for f, a in zip(peak_freqs, peak_amps):
    # 回転周波数の倍数か判定
    order = f / fr
    if abs(order - round(order)) < 0.05 and round(order) in [1, 2, 3]:
        status = f"正常高調波 ({round(order)}X)"
    else:
        status = f"★ 警告: 異常周波数検知!"
    print(f"周波数: {f:6.1f} Hz | 振幅: {a:5.3f} g | 判定: {status}")

# ---------------------------------------------------------
# 5. 可視化
# ---------------------------------------------------------
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t[:500], raw_vibration[:500], color='tab:blue')
plt.title("生振動波形 (最初の0.1秒間)")
plt.xlabel("時間 [s]")
plt.ylabel("加速度 [g]")
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(frequencies, amplitude_spectrum, label="振幅スペクトル", color='tab:blue', alpha=0.8)
plt.plot(frequencies, dynamic_threshold, label="動的検知閾値 (Baseline + 5σ)", color='tab:red', linestyle='--')
plt.scatter(peak_freqs, peak_amps, color='tab:red', marker='x', s=80, label="検知ピーク")

for f, a in zip(peak_freqs, peak_amps):
    plt.annotate(f"{f:.1f}Hz", (f, a), textcoords="offset points", xytext=(0, 8), ha='center', fontsize=9)

plt.xlim(0, 500) # 関心領域 (0〜500 Hz)
plt.title("FFT 振幅スペクトルと異常ピーク検知")
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅 [g]")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

---

4. エッジデバイス・現場導入における技術的勘所

Raspberry Pi、NVIDIA Jetson、産業用エッジゲートウェイで本アルゴリズムを稼働させる場合、以下の最適化が求められます。

1. データ長（N）は2の冪乗（$2^k$）に固定する:
   FFTアルゴリズム（Cooley-Tukey）はデータ点数が $N = 1024, 2048, 4096$ などの $2^n$ である場合に演算効率が最大化されます。
2. 周波数分解能（$\Delta f$）の設計:
   周波数分解能は $\Delta f = \frac{f_s}{N} = \frac{1}{\text{計測時間}}$ で決まります。近接する特徴周波数（例: 29.5 Hz と 30.0 Hz）を分離したい場合は、計測時間を最低2秒（$\Delta f = 0.5\text{ Hz}$）以上確保する必要があります。
3. エッジでの一次判定とクラウド連携のハイブリッド:
   生波形（数kHzの高サンプリングデータ）をクラウドへ常時アップロードすると通信費が肥大化します。エッジ側でFFTとピーク検知を実行し、「異常ピークの周波数・振幅」および「警告発生時の生波形（前後の数秒分）」のみをMQTT/TLS経由でクラウドへ通知する構成が最もコスト対効果に優れます。

---

まとめ

振動データの周波数解析は、単なる「機械学習モデルへの全投入」に先立って実施すべき、物理モデルに基づいた最も確実で解釈性の高い異常検知手法です。

• **直線性・DCオフセットの除去と窓関数（Hanning）**で正確なスペクトルを生成する。
• **動的閾値（ベースライン統計）**により、設備の個体差や負荷変動によるベースノイズの変化に追従する。
• 回転数 $f_r$ との幾何学的関連性をルールベースで照合し、異常部位（ベアリングか、カップリングか）を特定する。

NKKTech Globalでは、産業用センサー選定からエッジゲートウェイでのDSP（デジタル信号処理）、クラウド側でのダッシュボード構築まで、製造業の予兆保全を一気通貫で支援しています。

---

お問い合わせ先

産業用IoT・振動解析アルゴリズムの実装、エッジAIの現場導入、スマートファクトリー化のご相談は、下記よりお気軽にお問い合わせください。

• Webサイト: https://nkktech.com/
• メール: [e-mail]
• LinkedIn: https://www.linkedin.com/company/nkktech

著者：NKKTech Global 技術チーム
私たちは、デジタル信号処理と先端AI技術を融合させ、製造現場のダウンタイムゼロ化に貢献します。