移動平均や標準偏差によるトレンド監視は、多くの現場で導入されている王道の予兆検知手法です。
しかし、統計量が教えてくれるのは、主に値の大きさやばらつきの変化です。
一方で、制御系データ（車速やエンジン回転数など）に潜む異常は、数値の増減ではなく、制御のリズムや微細な揺らぎとして現れることがあります。
一般に、100ms周期のような低サンプリング環境ではFFT（高速フーリエ変換）は有効ではないと考えられがちです。
確かに、物理的に正確な周波数解析は困難ですが、FFTを「正しい周波数を求めるための手法」ではなく、波形の型崩れを検知するセンサーとして捉え直したらどうでしょうか。
たとえ物理的に厳密ではない成分であっても、そこに一貫したパターンが現れるなら、それは異常の“指紋”になり得ます。
本記事では、従来の統計的手法を補完し、低頻度データからリズムの乱れを抽出するハイブリッドな予兆検知手法を、Pythonによる実装例とともに解説します。
目次
なぜ統計量だけでなく「FFT」を使うのか？
予兆検知の現場で最も普及しているのは、移動平均や移動標準偏差（移動指標）による監視です。これらは「異常の有無」を察知する体温計のような役割を果たし、値の変動（バラつき具合）は把握できますが、「バラつきの質」や「リズムの乱れ」までは捉えられません。
FFTを使うことで、その2つを明確に識別できるようになります。
FFTを使うメリット
「バラつきの質」が区別できる
統計量、特に標準偏差（ばらつき）が教えてくれるのは、データの「激しさ」の総量です。 例えば、以下の2つのケースを想定してみます。
- ケースA（外乱）： 路面の凹凸で、車速がランダムにガサガサと揺れている。
- ケースB（予兆）： 制御系が不安定になり、一定の周期で「波打つ」ような挙動が出始めた。
統計量で見ると、どちらも「標準偏差が増えた」という一つの数値にまとめられてしまいます。しかし、FFTを使えば、ケースAは全帯域のノイズ上昇、ケースBは「特定の周波数へのパワー集中」として明確に区別できます。
制御のリズム（型崩れ）を数値化できる
「数値は正常範囲内だけど、波形がなんとなくカクカクしている」……ベテランが経験則で感じるこうした違和感の正体は、波形に含まれる「高調波」です。
きれいなサイン波に近い挙動が、メカのガタや制御の遅れによって「カクついた波形（矩形波に近い形）」に変化すると、FFTスペクトル上には基本波の整数倍の場所に「山」が立ち並びます。これは統計量では「わずかなノイズ増」として埋もれてしまう程度の変化ですが、FFTなら「リズムの型崩れ」として鮮明に捉えることができます。
FFTを使った具体例
百聞は一見にしかず。実際にデータで確認してみましょう。
以下のグラフは、ある制御系（架空）の出力値を60秒間記録したものです。30秒を境に、正常状態から異常状態へと切り替わっています。
異常区間では「1.2Hzのハンチング（制御発振の予兆）」が密かに混入していますが、波形を目視しても正常と異常の区別はほぼつきません。
では、移動標準偏差を見るとどうでしょう。赤い破線（異常発生ポイント）の前後で、値は多少上下していますが、「これが異常だ」と断言できるほどの変化は読み取れません。
次のグラフは、同じデータをFFTで分析した結果です。正常時は0.5Hzの成分のみが見られるのに対し、異常時には1.2Hz付近と2Hz付近に新たな周波数成分が出現していることが分かります。
次ののグラフは、1.2Hz成分の強度をリアルタイムで追跡した結果です。
薄いピンクの線が生データ（フレームごとの1.2Hz成分の強度）で、そのままではノイズが多くトレンドが読みにくい状態です。そこに移動平均をかけて平滑化したものが、太い赤い線です。
異常発生ポイント（黒い破線・30秒）を境に強度が上昇し始め、40秒付近で正常区間の統計から算出した判定閾値（μ+3σ＝0.087）を超えたことで、異常として自動検知されています（赤い塗りつぶし領域）。
移動標準偏差では「なんとなく上がっているかも」という程度にしか見えなかった変化が、FFTで1.2Hz成分に着目することで、閾値を使った客観的な自動判定へと昇華できることが分かります。
こちらをクリックすると、サンプルプログラムが表示されます。
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.fft import fft, fftfreq
# ===================================================
# 共通設定
# ===================================================
plt.rcParams['font.family'] = 'Meiryo'
plt.rcParams['axes.unicode_minus'] = False
ANOMALY_TIME = 30 # 異常発生ポイント（秒）
# ===================================================
# 共通描画関数
# ===================================================
def plot_timeseries(t, signal, title, ylabel='値', anomaly_time=ANOMALY_TIME):
"""時系列データを描画する"""
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(t, signal, color='gray', alpha=0.6)
ax.axvline(anomaly_time, color='red', ls='--')
ax.set_title(title, fontsize=14)
ax.set_xlabel('時間 [秒]')
ax.set_ylabel(ylabel)
plt.tight_layout()
plt.show()
def plot_moving_stat(t, stat_series, title, label, color='green', anomaly_time=ANOMALY_TIME):
"""移動統計量（移動標準偏差など）を描画する"""
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(t, stat_series, color=color, label=label)
ax.axvline(anomaly_time, color='red', ls='--')
ax.set_title(title, fontsize=14)
ax.set_xlabel('時間 [秒]')
ax.legend()
plt.tight_layout()
plt.show()
def plot_fft_comparison(xf, psd_normal, psd_anomaly,
title_normal='正常時の成分', title_anomaly='異常時の成分'):
"""正常／異常のFFTスペクトルを左右に並べて描画する"""
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5), sharey=True)
ax1.stem(xf, psd_normal, basefmt=" ", markerfmt="b,")
ax1.set_title(title_normal, fontsize=12)
ax1.set_xlabel('周波数 [Hz]')
ax2.stem(xf, psd_anomaly, basefmt=" ", markerfmt="r,")
ax2.set_title(title_anomaly, fontsize=12)
ax2.set_xlabel('周波数 [Hz]')
plt.tight_layout()
plt.show()
def plot_fft_feature(t_feat, features, smoothed, thresh, title,
feat_label='特徴量（生データ）', anomaly_time=ANOMALY_TIME,
ylim=None):
"""FFT特徴量の推移・移動平均・閾値・異常検知領域を描画する"""
fig, ax = plt.subplots(figsize=(12, 7))
ax.plot(t_feat, features, color='salmon', alpha=0.3, label=feat_label)
ax.plot(t_feat, smoothed, color='red', linewidth=3, label='移動平均（平滑化トレンド）')
ax.axhline(thresh, color='blue', linestyle=':', linewidth=2,
label=f'判定閾値 (トレンドのμ+3σ: {thresh:.3f})')
ax.axvline(anomaly_time, color='black', linestyle='--', linewidth=2, label='異常発生ポイント')
ax.fill_between(t_feat, smoothed, thresh,
where=(smoothed > thresh), color='red', alpha=0.2, label='異常検知領域')
ax.set_title(title, fontsize=14)
ax.set_ylabel('異常指紋の強度')
ax.set_xlabel('時間 [秒]')
ax.grid(True, alpha=0.2)
ax.legend(loc='upper left', frameon=True, shadow=True)
if ylim:
ax.set_ylim(ylim)
plt.tight_layout()
plt.show()
# ===================================================
# FFT計算関数
# ===================================================
def get_fft(s, fs):
"""FFTスペクトルを返す（片側）"""
n = len(s)
yf = np.abs(fft(s))[:n//2]
xf = fftfreq(n, 1/fs)[:n//2]
return xf, 2.0/n * yf
def calc_band_energy(signal, fs, freq_low, freq_high, window=50):
"""スライディングウィンドウで特定帯域のエネルギーを計算する"""
features = []
for i in range(len(signal) - window):
segment = signal[i : i + window]
xf, psd = get_fft(segment, fs)
idx = np.where((xf >= freq_low) & (xf <= freq_high))
features.append(np.mean(psd[idx]))
return features
# ===================================================
# データ生成
# ===================================================
fs = 10.0 # サンプリング周波数 100ms
t = np.linspace(0, 60, 600, endpoint=False)
split = 300
np.random.seed(45)
# 正常：0.5Hz + ノイズ
sig_n = np.sin(2 * np.pi * 0.5 * t[:split]) + np.random.normal(0, 0.2, split)
# 異常：0.5Hz + 1.2Hz(ハンチング) + 12Hz(エイリアシング→2Hz) + ノイズ
sig_a = (np.sin(2 * np.pi * 0.5 * t[split:]) +
0.2 * np.sin(2 * np.pi * 1.2 * t[split:]) +
0.15 * np.sin(2 * np.pi * 12.0 * t[split:]) +
np.random.normal(0, 0.15, split))
sig_full = np.concatenate([sig_n, sig_a])
# 移動標準偏差
mv_std = pd.Series(sig_full).rolling(50).std()
# ===================================================
# グラフ描画
# ===================================================
# ① 時系列データ
plot_timeseries(t, sig_full, '① 時系列データ（正常→異常）')
# ② 統計量による監視
plot_moving_stat(t, mv_std, '② 統計量による監視（変化が緩慢で判別困難）', label='移動標準偏差')
# ③ FFTスペクトル比較（正常 vs 異常）
xf, psd_n = get_fft(sig_n, fs)
_, psd_a = get_fft(sig_a, fs)
plot_fft_comparison(xf, psd_n, psd_a,
title_normal='正常時の成分（0.5Hzのみ）',
title_anomaly='異常時の成分（指紋が出現！）')
# ④ FFT特徴量の推移と自動判定
WINDOW = 50
FREQ_LOW = 1.5
FREQ_HIGH = 2.5
SMOOTH_WIN = 30
features = calc_band_energy(sig_full, fs, FREQ_LOW, FREQ_HIGH, window=WINDOW)
smoothed = pd.Series(features).rolling(window=SMOOTH_WIN, center=True).mean()
# 正常区間（先頭250点）から閾値を算出
normal_trend = smoothed.dropna()[:250]
thresh = np.mean(normal_trend) + 3 * np.std(normal_trend)
plot_fft_feature(
t_feat = t[WINDOW:],
features = features,
smoothed = smoothed,
thresh = thresh,
title = '③ FFT特徴量の推移と自動判定（移動平均ベースの3σ閾値）',
feat_label = '1.2Hz成分の強度（生データ）',
ylim = (0.04, 0.12)
)
低サンプリング環境でFFTを使うための「逆転の発想」
低サンプリング環境（100ms周期など）でFFTを活用するための「逆転の発想」とは、一言で言えば「エイリアシング（折り返し雑音）を敵ではなく、異常検知のセンサー（味方）として利用する」ことです。
通常の信号解析では、サンプリング定理に基づき、エイリアシングは「正しく解析できないゴミ」として忌避されます。しかし、現場の制約を突破するためには、以下の3つの視点を持つことが重要です。
1. サンプリング定理の「壁」を「指紋」に変える
教科書的には、10Hzサンプリングなら5Hzまでの振動しか見えません。しかし、異常によって発生する高周波（例：12Hz）は、消えてなくなるのではなく、低い周波数（例：2Hz）に「化けて」出現します。
- 一般的な常識
エイリアシングは情報の欠落であり、解析不能とする。 - 逆転の発想
特定の低周波領域に現れた「謎の山」を、高周波異常の固有の指紋として定義し、監視対象とする。
2. 「正体」不明でも「変化」は真実
何の予備知識もなければ、その山が本物の1.2Hz（制御の乱れ）か、化けて出た2.0Hz（高周波振動）かは区別できません。しかし、予兆検知においては、物理的な正解を当てることよりも「いつもと違う状態」を見つけることが最優先です。
- 多角的な監視
統計量（標準偏差）では埋もれてしまう微細な変化も、複数の「山」の出現として捉えることで、検知の解像度が格段に上がります。
3. 「一目瞭然」にするための後処理の徹底
FFTの計算結果（スカラー値）をそのままプロットしても、ノイズで激しく上下して判定には使えません。そこで、「FFTで特定のリズムを抜き出す」工程の後に、「移動平均でトレンドを浮き彫りにする」という二段構えの処理を行います。
- 生データは嘘をつく
一瞬のノイズで山が立つのに惑わされないようにする。 - トレンドで確信を得る
移動平均後のデータに対して、正常時のバラツキから算出した 3σ 閾値を適用することで、統計的にも明らかな「異常検知領域」を可視化します。
FFTを使った予兆検知の方法
FFTで得られたスペクトルから「予兆」を抜き出す手法は、単にピークを見るだけではありません。監視対象（ギヤ、ベアリング、制御系など）に合わせて、以下のような手法を使い分けます。
| 手法 | 概要 | 主な用途 |
|---|---|---|
| 特定帯域エネルギー | 特定の周波数範囲のパワーを合計・平均して監視する | 異常の指紋（エイリアシング）や特定回転数に起因する振動の監視 |
| 重心周波数（スペクトルセントロイド） | スペクトル全体の重みがどの周波数帯にあるかを計算する | 摩耗進行による振動の「高域寄り」といったマクロな劣化トレンドの把握 |
| ギヤ側帯波（サイドバンド）解析 | ギヤ噛み合い周期の周囲に現れる側帯波の成長を監視する | ギヤの歯こぼれ・摩耗の初期検知 |
| エンベロープ解析 | 包絡線処理により高周波共振に隠れた微細な打撃周期を抽出する | ベアリングの内輪・外輪の傷の検知 |
それぞれの手法ごとのサンプルを実行するために必要な共通関数を紹介しておきます。
実行する場合は、必ず下記のソースコードを先頭に張り付けてください。
ここをクリックすると、共通関数のソースコードが表示されます。
# ==================================================================
# FFT予兆検知 共通関数ライブラリ（トレンド監視機能付き）
#
# 【概要】
# 低サンプリング環境（100ms周期など）において、FFTを用いた予兆検知を
# 行うための共通関数群です。以下の3レイヤーで構成されています。
#
# [1] データ生成 : ダミーデータの生成（検証・デモ用）
# [2] 解析ロジック : FFTスペクトルの算出
# [3] グラフ描画 : 時系列・スペクトル比較・トレンド監視の可視化
#
# 【設計方針】
# monitor_trend() はアルゴリズムをコールバック関数として受け取る構造になっており、
# 特定帯域エネルギー・重心周波数・サイドバンド解析など、手法ごとの特徴量算出ロジックを
# 差し替えるだけで同じ監視パイプラインを再利用できます。
#
# 【依存ライブラリ】
# numpy, matplotlib, pandas, scipy
# ==================================================================
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.fft import fft, fftfreq
# 日本語フォント設定（Windows環境向け。Mac/Linuxは 'IPAexGothic' などに変更）
plt.rcParams['font.family'] = 'Meiryo'
plt.rcParams['axes.unicode_minus'] = False
# ==================================================================
# 1. データ生成
# ==================================================================
def generate_mechanical_fault_data(fs=10.0, duration=60, seed=123):
"""
回転体の機械的ガタによる高周波衝撃が、低サンプリングによって
低周波成分（エイリアシング）として現れるダミーデータを生成する。
【エイリアシングの仕組み】
サンプリング周波数 fs=10Hz の場合、ナイキスト周波数は 5Hz。
18Hz の信号は折り返しにより |18 - 10×2| = 2Hz として観測される。
この「化けて出た 2Hz」を異常の指紋（フィンガープリント）として
監視するのが本手法の核心。
Parameters
----------
fs : float
サンプリング周波数 [Hz]（デフォルト: 10.0 = 100ms周期）
duration : int
計測時間 [秒]（デフォルト: 60秒）
seed : int
乱数シード（再現性確保用）
Returns
-------
t : ndarray
時間軸 [秒]
sig_full : ndarray
正常区間 + 異常区間を結合した時系列信号
split : int
正常→異常の切り替わりインデックス（= duration/2 のサンプル数）
使用例
------
>>> t, sig, split = generate_mechanical_fault_data(fs=10.0, duration=60)
>>> sig_normal = sig[:split] # 正常区間
>>> sig_anomaly = sig[split:] # 異常区間
"""
np.random.seed(seed)
t = np.linspace(0, duration, int(fs * duration), endpoint=False)
split = len(t) // 2
# 正常区間：主回転成分（0.5Hz）+ 穏やかなガウスノイズ（σ=0.1）
sig_n = np.sin(2 * np.pi * 0.5 * t[:split]) + np.random.normal(0, 0.1, split)
# 異常区間：主回転成分（0.5Hz）+ 機械的ガタによる高周波衝撃（18Hz）+ ノイズ
# 18Hz はナイキスト(5Hz)を超えるため、10Hzサンプリングでは 2Hz に折り返される
# 振幅 0.25 は統計量（標準偏差）では埋もれやすいが、FFTスペクトルでは明確に現れる
impact_component = 0.25 * np.sin(2 * np.pi * 18.0 * t[split:])
sig_a = (np.sin(2 * np.pi * 0.5 * t[split:]) +
impact_component +
np.random.normal(0, 0.15, len(t) - split))
return t, np.concatenate([sig_n, sig_a]), split
# ==================================================================
# 2. 解析ロジック
# ==================================================================
def get_fft_spectrum(sig_data, fs):
"""
FFTを実行し、片側の周波数軸と正規化済み振幅スペクトルを返す。
【正規化について】
生のFFT出力は信号長 n に比例して大きくなるため、
2.0/n を乗じて実際の振幅値（Peak振幅）に正規化している。
※ DC成分（0Hz）は 1.0/n が正確だが、予兆検知では交流成分を
主に扱うため、簡略化のため統一して 2.0/n を適用している。
【使用上の注意】
- 入力信号は窓関数なし（矩形窓）で処理される。
スペクトルリーク（隣接周波数への漏れ）が気になる場合は、
事前に np.hanning(n) などの窓関数を乗じること。
- 周波数分解能は fs/n [Hz]。精度を上げるには信号長 n を大きくする。
Parameters
----------
sig_data : array-like
解析対象の時系列信号（1次元）
fs : float
サンプリング周波数 [Hz]
Returns
-------
xf : ndarray
周波数軸 [Hz]（0 〜 ナイキスト周波数）
yf : ndarray
正規化済み振幅スペクトル（xf と同じ長さ）
使用例
------
>>> xf, yf = get_fft_spectrum(sig_data=sig_normal, fs=10.0)
>>> peak_freq = xf[np.argmax(yf)] # 最も強い周波数成分を取得
"""
n = len(sig_data)
yf = np.abs(fft(sig_data))[:n//2]
xf = fftfreq(n, 1/fs)[:n//2]
return xf, 2.0 / n * yf
# ==================================================================
# 3. グラフ描画
# ==================================================================
def plot_raw_timeseries(t, sig_data, split_idx, title="① 時系列データの確認"):
"""
時系列信号を描画し、正常区間と異常区間を視覚的に確認するためのグラフを表示する。
【用途】
データ生成直後の目視確認、またはブログ・報告書用の説明図として使用。
「統計量だけでは異常が分かりにくい」ことを示すための問いかけグラフ。
Parameters
----------
t : ndarray
時間軸 [秒]
sig_data : ndarray
時系列信号（t と同じ長さ）
split_idx : int
正常→異常の切り替わりインデックス
（generate_mechanical_fault_data の戻り値 split をそのまま渡す）
title : str
グラフタイトル
使用例
------
>>> t, sig, split = generate_mechanical_fault_data()
>>> plot_raw_timeseries(t, sig, split)
"""
plt.figure(figsize=(12, 4))
plt.plot(t, sig_data, color='gray', alpha=0.6, label='観測データ')
plt.axvline(t[split_idx], color='red', linestyle='--', label='異常発生ポイント')
plt.title(title, fontsize=14)
plt.ylabel("振幅")
plt.xlabel("時間 [秒]")
plt.legend(loc='upper right')
plt.grid(True, alpha=0.2)
plt.tight_layout()
plt.show()
def plot_spectrum_comparison(xf, yf_normal, yf_anomaly, title="周波数成分の比較"):
"""
正常時と異常時のFFTスペクトルを左右に並べて比較表示する。
【読み方のポイント】
- 正常時：主回転成分（0.5Hz）のみに「山」が立つ
- 異常時：エイリアシングにより 2Hz 付近に新たな「山」が出現する
→ この「余計な山」が異常の指紋（フィンガープリント）
【表示範囲について】
ナイキスト周波数（fs/2 = 5Hz）までを表示対象とし、
それ以上の周波数は意味を持たないため xlim で除外している。
Parameters
----------
xf : ndarray
周波数軸 [Hz]（get_fft_spectrum の戻り値をそのまま渡す）
yf_normal : ndarray
正常時の振幅スペクトル
yf_anomaly : ndarray
異常時の振幅スペクトル
title : str
グラフ全体のタイトル
使用例
------
>>> xf, yf_n = get_fft_spectrum(sig[:split], fs=10.0)
>>> _, yf_a = get_fft_spectrum(sig[split:], fs=10.0)
>>> plot_spectrum_comparison(xf, yf_n, yf_a)
"""
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5), sharey=True)
# 正常時スペクトル（青）
ax1.stem(xf, yf_normal, basefmt=" ", markerfmt="b,")
ax1.set_title("正常時スペクトル（0.5Hz主成分）", fontsize=12)
# 異常時スペクトル（赤）：2Hz付近に指紋（エイリアシング）が現れているはず
ax2.stem(xf, yf_anomaly, basefmt=" ", markerfmt="r,")
ax2.set_title("異常時スペクトル", fontsize=12)
# 両軸共通設定：ナイキスト周波数（5Hz）までを表示
for ax in [ax1, ax2]:
ax.set_xlabel("周波数 [Hz]")
ax.set_xlim(0, 5)
ax.grid(True, alpha=0.2)
plt.suptitle(title, fontsize=14)
plt.tight_layout()
plt.show()
def monitor_trend(sig_data, fs, algo_callback, window_size=50, title="予兆検知モニター"):
"""
スライディングウィンドウでFFT特徴量を時系列算出し、トレンド監視グラフを表示する。
【処理フロー】
1. スライディングウィンドウで信号を切り出す（window_size サンプルずつ）
2. 各ウィンドウに algo_callback を適用して特徴量（スカラー値）を算出
3. 移動平均で平滑化してトレンドを浮き彫りにする
4. 正常区間の統計（μ+3σ）から判定閾値を自動算出
5. plot_trend_monitor() で結果を可視化
【algo_callback の仕様】
algo_callback(segment, fs) -> float の形式で特徴量を1つ返す関数を渡す。
例：特定帯域エネルギーを使う場合
def band_energy(seg, fs):
xf, yf = get_fft_spectrum(seg, fs)
idx = np.where((xf >= 1.5) & (xf <= 2.5))
return np.mean(yf[idx])
monitor_trend(sig, fs=10.0, algo_callback=band_energy)
この構造により、手法（特定帯域・重心周波数・サイドバンドなど）を
差し替えるだけで同じ監視パイプラインを再利用できる。
Parameters
----------
sig_data : ndarray
監視対象の時系列信号（正常＋異常の全区間）
fs : float
サンプリング周波数 [Hz]
algo_callback : callable
特徴量算出関数。シグネチャ: (segment: ndarray, fs: float) -> float
window_size : int
スライディングウィンドウのサンプル数（デフォルト: 50 = 5秒分）
title : str
グラフタイトル
使用例
------
>>> def band_energy(seg, fs):
... xf, yf = get_fft_spectrum(seg, fs)
... return np.mean(yf[np.where((xf >= 1.5) & (xf <= 2.5))])
>>> t, sig, split = generate_mechanical_fault_data()
>>> monitor_trend(sig, fs=10.0, algo_callback=band_energy)
"""
# スライディングウィンドウで特徴量を算出
features = [algo_callback(sig_data[i:i + window_size], fs)
for i in range(len(sig_data) - window_size)]
# 移動平均でトレンドを平滑化（window=30、前後15点を使った中央移動平均）
smoothed = pd.Series(features).rolling(window=30, center=True).mean()
# 正常区間（先頭250点）の統計からμ+3σ閾値を算出
# ※ 250点 = window_size=50 を除いた後の正常区間の概算サンプル数
normal_part = smoothed.dropna()[:250]
thresh = normal_part.mean() + 3 * normal_part.std()
# 時間軸：ウィンドウの開始位置を秒に変換
t_feat = np.arange(len(features)) / fs + (window_size / fs)
plot_trend_monitor(t_feat, features, smoothed, thresh, title)
def plot_trend_monitor(t, raw, smooth, thresh, title):
"""
FFT特徴量の推移・平滑化トレンド・判定閾値・異常検知領域を1枚のグラフに描画する。
【各要素の役割】
- 薄いサーモン色（生特徴量）: フレームごとのFFT特徴量。ノイズが多く判定には使わない。
- 太い赤線（トレンド） : 移動平均後のトレンド。これを閾値と比較して判定する。
- 青い点線（判定閾値） : 正常区間のμ+3σ。統計的に「ここを超えたら異常」の基準線。
- 黒い破線（異常発生ポイント）: データ上の正常→異常の切り替わり時刻（答え合わせ用）。
- 赤い塗りつぶし（異常検知領域）: トレンドが閾値を超えた区間。自動判定の結果。
Parameters
----------
t : ndarray
時間軸 [秒]（monitor_trend 内で算出した t_feat を渡す）
raw : list or ndarray
生特徴量（平滑化前）
smooth : pd.Series
移動平均後の平滑化トレンド
thresh : float
判定閾値（μ+3σ）
title : str
グラフタイトル
Notes
-----
本関数は monitor_trend() から内部的に呼び出されることを想定しているが、
外部から特徴量・閾値を直接渡して単独で使用することも可能。
"""
plt.figure(figsize=(12, 6))
# 生特徴量（薄い色で背景的に表示）
plt.plot(t, raw, color='salmon', alpha=0.3, label='生特徴量')
# 平滑化トレンド（濃い色・太線で強調）
plt.plot(t, smooth, color='red', linewidth=3, label='トレンド（移動平均）')
# 判定閾値（μ+3σ）
plt.axhline(thresh, color='blue', linestyle=':', linewidth=2,
label=f'判定閾値(3σ): {thresh:.3f}')
# 異常発生ポイント（答え合わせ用の垂直線）
plt.axvline(30, color='black', linestyle='--', label='異常発生ポイント')
# 閾値超過区間を塗りつぶして異常検知領域を強調
plt.fill_between(t, smooth, thresh,
where=(smooth > thresh), color='red', alpha=0.2, label='異常検知領域')
plt.title(title, fontsize=14)
plt.xlabel("時間 [秒]")
plt.legend(loc='upper left', frameon=True, shadow=True)
plt.grid(True, alpha=0.2)
plt.tight_layout()
plt.show()
# --生データをグラフ化する場合は下記コメントを外して下さい。
# fs = 10.0
# t_raw, sig_full, split_idx = generate_mechanical_fault_data(fs)
# 時系列表示
# plot_raw_timeseries(t_raw, sig_full, split_idx)
今回のサンプルデータは、振動の「音色の変化（質の変化）」を再現しています。
正常区間 (0〜30秒): 主回転成分（低周波の0.5Hz）に、穏やかなホワイトノイズが乗っている状態です。このとき、スペクトルの重みは低周波側に集中しています。
異常区間 (30〜60秒): 18Hzという高い周波数の衝撃成分（ガタつき）を混入させています。10Hzサンプリング環境では、これが2Hzのエイリアシングとして現れますが、重要なのは「スペクトル全体の分布が右側（高周波側）へ引っ張られる」という点です。
1. 特定帯域エネルギー（今回の手法）
今回実装した手法です。特定の周波数範囲（例：1.5Hz〜2.5Hz）のパワーを合計、または平均して監視します。
- 用途
今回のような「異常の指紋（エイリアシング）」や、特定の回転数に起因する振動の監視に最適です。 - 逆転の発想
正確な周波数が特定できなくても、「この辺りがざわついてきた」という変化を捉えることができます。
# ===================================================================
# ここに共通関数のソースコードを貼り付け下さい
# ===================================================================
# ---------------------------------------------------------
# 4. 個別のアルゴリズム（コールバック関数）
# ---------------------------------------------------------
def algo_band_energy(segment, fs):
xf, psd = get_fft_spectrum(segment, fs)
idx = np.where((xf >= 1.5) & (xf <= 2.5))
return np.mean(psd[idx]) if len(idx[0]) > 0 else 0
# ---------------------------------------------------------
# 5. メイン実行
# ---------------------------------------------------------
fs = 10.0
t_raw, sig_full, split_idx = generate_mechanical_fault_data(fs)
# スペクトル比較
xf, psd_n = get_fft_spectrum(sig_full[:split_idx], fs)
_, psd_a = get_fft_spectrum(sig_full[split_idx:], fs)
plot_spectrum_comparison(xf, psd_n, psd_a)
# トレンド監視
monitor_trend(sig_full, fs, algo_band_energy, title="FFT：特定帯域（2Hz指紋）による自動検知モニター")
2. 重心周波数（スペクトルセントロイド）
スペクトル全体の「重み」がどこにあるかを計算します。
- 用途
摩耗の進行などで、全体的に振動が「高域寄り」になったというマクロな変化を捉えます。 - メリット
個別の山を見る必要がないため、ノイズに強く、全体的な劣化トレンドを追うのに適しています。
# ===================================================================
# ここに共通関数のソースコードを貼り付け下さい
# ===================================================================
# ---------------------------------------------------------
# 4. 個別のアルゴリズム（コールバック関数：重心周波数）
# ---------------------------------------------------------
def algo_spectral_centroid(segment, fs):
"""
【重心周波数】
スペクトル全体の「中心」がどこにあるかを算出する。
摩耗が進んで「ガサガサ」と高域成分が増えた時などに数値が上昇する。
"""
# 1. FFTスペクトルを取得（共通関数を利用）
xf, psd = get_fft_spectrum(segment, fs)
# 2. 重心周波数の計算 formula: Σ(周波数 * 強度) / Σ(強度)
# 分母が0になるのを防ぐため微小値を加算
sum_psd = np.sum(psd)
if sum_psd == 0:
return 0
centroid = np.sum(xf * psd) / sum_psd
return centroid
# ---------------------------------------------------------
# 5. メイン実行
# ---------------------------------------------------------
# データ生成（18Hzの衝撃が混入し、全体的に高域成分が増えるパターン）
fs = 10.0
t_raw, sig_full, _ = generate_mechanical_fault_data(fs)
# 共通基盤を呼び出し
monitor_trend(
sig_data=sig_full,
fs=fs,
algo_callback=algo_spectral_centroid,
window_size=50,
title="重心周波数：スペクトル全体の高域シフトを監視"
)
3. ギヤ側帯波（サイドバンド）解析
ギヤの噛み合い周期の周りに現れる「小さな山」を監視します。
- 用途
ギヤの歯こぼれや摩耗の検知。 - 特徴
ギヤが傷むと、中心の山よりも先にその「横の山（側帯波）」が成長します。この比率を追うことで、初期段階での予兆検知が可能になります。
# ---------------------------------------------------------
# 4. 個別のアルゴリズム（コールバック関数：ギヤ側帯波）
# ---------------------------------------------------------
def algo_gear_sideband_energy(segment, fs):
"""
【ギヤ側帯波解析】
噛み合い周波数（Mesh）の左右に現れるサイドバンド成分を抽出する。
ギヤの歯面のキズや摩耗が進むと、この「横の山」が大きく成長する。
"""
xf, psd = get_fft_spectrum(segment, fs)
# 例：噛み合い周波数が 1.2Hz、回転周期（サイドバンドの間隔）が 0.5Hz の場合
mesh_freq = 1.2
rotation_freq = 0.5
# ターゲット：メッシュ周波数の±回転周期（0.7Hz付近 と 1.7Hz付近）
# 実務では高調波（2次、3次）や複数のサイドバンドを合算することもある
sidebands = [mesh_freq - rotation_freq, mesh_freq + rotation_freq]
total_energy = 0
margin = 0.2 # 抽出する帯域の幅
for fb in sidebands:
idx = np.where((xf >= fb - margin) & (xf <= fb + margin))
if len(idx[0]) > 0:
total_energy += np.mean(psd[idx])
return total_energy
# ---------------------------------------------------------
# 5. メイン実行
# ---------------------------------------------------------
# データ生成（1.2Hzの周辺に0.7Hzや1.7Hzの山が立つパターンを想定）
fs = 10.0
t_raw, sig_full, _ = generate_mechanical_fault_data(fs)
# 共通基盤を呼び出し
monitor_trend(
sig_data=sig_full,
fs=fs,
algo_callback=algo_gear_sideband_energy,
window_size=100, # 解像度が必要なため窓は長めを推奨
title="ギヤ解析：側帯波（サイドバンド）エネルギーの推移"
)
4. エンベロープ解析（ベアリング診断）
高周波の共振に隠れた、ベアリング特有の「コンコン」という微細な打撃周期を抽出します。
- 用途
ベアリングの内輪・外輪の傷。 - 強み
そのままのFFTではノイズに埋もれてしまう小さなキズのサインを、包絡線処理によって「指紋」として浮き彫りにします。
from scipy.signal import hilbert
# ---------------------------------------------------------
# 4. 個別のアルゴリズム（コールバック関数：エンベロープ）
# ---------------------------------------------------------
def algo_envelope_fault_amp(segment, fs):
"""
【エンベロープ解析】
高周波の共振に隠れたベアリングの打撃周期を抽出する
"""
# 1. ヒルベルト変換で包絡線（エンベロープ）を取得
analytic_signal = hilbert(segment)
envelope = np.abs(analytic_signal)
# 2. 包絡線から直流成分（平均値）を除去
envelope_detrended = envelope - np.mean(envelope)
# 3. 包絡線をFFTして、損傷周期（指紋）の強度を確認
xf, psd = get_fft_spectrum(envelope_detrended, fs)
# ベアリングの損傷周期が出現しやすい帯域（例: 2.5Hz〜3.5Hz）を監視
idx = np.where((xf >= 2.5) & (xf <= 3.5))
return np.mean(psd[idx]) if len(idx[0]) > 0 else 0
# ---------------------------------------------------------
# 5. メイン実行
# ---------------------------------------------------------
# データ生成（ベアリング故障を模した衝撃波形データが必要な場合は生成関数を差し替え）
fs = 10.0
t_raw, sig_full, split_idx = generate_mechanical_fault_data(fs)
# 共通基盤を呼び出し
monitor_trend(
sig_data=sig_full,
fs=fs,
algo_callback=algo_envelope_fault_amp,
window_size=100, # 解析精度のため少し長めの窓を推奨
title="エンベロープ解析：ベアリング損傷周期の自動監視"
)
まとめ
本記事では、100ms周期（10Hz）という低サンプリング環境の制約を逆手に取り、FFTを予兆検知に活用する「ハイブリッド手法」を解説しました。
現場で真に役立つシステムを構築するためのポイントは、以下の3点に集約されます。
- 統計量で見えない「異常の質」を可視化する
移動標準偏差などの統計量は、異常の「大きさ」を測る体温計のような役割を果たしますが、その「質（リズムの乱れ）」までは教えてくれません。 単なる外乱によるノイズ増と、故障の前兆である周期的な挙動を明確に区別するには、周波数領域でのアプローチが不可欠です。 - エイリアシングを「異常の指紋」として味方につける
低サンプリング環境では高周波成分を直接観測することはできませんが、それは「消失」したのではなく、エイリアシングによって低周波側に「化けて」現れています。 この現象を敵ではなく、異常を知らせる固有の指紋（SOS）として捉え直すのが本手法の核心です。 - 「移動平均 」＋閾値で判定
FFTで抽出した特徴量は一瞬のノイズに敏感に反応するため、そのままでは誤報の原因となります。移動平均によりトレンドを平滑化し、閾値による判定（正常時の統計から算出した 3σなど）を行うことで検知を行います。
今回の記事が、皆さんの予兆検知のお役に立てれば幸いです。
コメント