Claude Codeで振動データをFFT解析 — 基準期間との差で兆候を探す

設備の振動波形は、時間軸のまま眺めても変化に気づきにくいデータです。周波数軸に直すと、回転数の整数倍に出る成分や、これまで無かった高周波のピークが見えてきます。Claude Codeには、この「波形をスペクトルにして基準期間と比べる」スクリプトを書かせ、手元で走らせて結果の表を読ませるのが向いています。

この記事では、SciPyのWelch法で振動CSVをパワースペクトル密度(PSD)に変換し、帯域ごとの差を表にするところまでを通しで動かします。判定そのものは保全担当者の仕事で、Claude Codeが受け持つのは計算と表づくりです。

Claude Codeに任せる範囲と、人が決める範囲

Welch法(Welch's method)とは、波形を重なりのある短い区間に分け、区間ごとのパワースペクトルを平均してスペクトル密度を推定する方法です。SciPyのscipy.signal.welchが実装していて、区間の長さ・重なり・窓関数・平均の取り方を引数で選べます。

単発のFFTと違い、区間平均でばらつきが減るため、「先月と今月でどこが変わったか」の比較に向きます。

PLCの時系列ログ(閾値超過や徐々に進む変化)の見方はClaude CodeでPLCのログ解析にあります。ここで扱うのは、加速度センサの生波形のような高いサンプリング周波数のデータを周波数領域で見る場合です。

使うデータと環境

動かした環境はPython 3、SciPy 1.18.1、NumPy 2.5.4です。実機のデータは載せられないので、次の条件の合成データで確かめました。

「基準より2次成分が増え、見慣れない高周波ピークが出た」という状況を人工的に作ったものです。実測データではないので、表の数値そのものに意味はありません。見るべきは、既知の変化をスクリプトが拾えるかどうかです。

import numpy as np
 
FS = 5000
rng = np.random.default_rng(0)
t = np.arange(60 * FS) / FS
 
def wave(a2x, tone):
    x = 1.0 * np.sin(2 * np.pi * 25 * t)          # 回転1次成分 25 Hz
    x += a2x * np.sin(2 * np.pi * 50 * t)         # 2次成分
    x += tone * np.sin(2 * np.pi * 1180 * t)      # 高周波の単独ピーク
    x += rng.normal(0, 0.2, t.size) + 0.5         # ノイズ + DCオフセット
    return x
 
for name, a2x, tone in [("baseline", 0.3, 0.0), ("current", 0.9, 0.25)]:
    np.savetxt(f"{name}.csv", np.c_[t, wave(a2x, tone)],
               delimiter=",", header="time_s,accel_g", comments="", fmt="%.6f")

Welch法で基準との差を出すスクリプト

比較スクリプトは次のとおりです。Claude Codeに書かせた後で、人が引数の意味を読んで確かめる前提で、短く保ちます。

import sys
import numpy as np
from scipy import signal
 
FS = 5000
BANDS = [(0, 20), (20, 40), (40, 60), (60, 500), (500, 2000)]  # Hz
 
def psd(path):
    x = np.loadtxt(path, delimiter=",", skiprows=1)[:, 1]
    return signal.welch(x, fs=FS, window="hann", nperseg=4096,
                        detrend="constant", average="median")
 
def band_power(f, p, lo, hi):
    m = (f >= lo) & (f < hi)
    return np.trapezoid(p[m], f[m])
 
def peaks(f, p, prominence_db=10):
    idx, _ = signal.find_peaks(10 * np.log10(p), prominence=prominence_db)
    return f[idx]
 
fb, pb = psd(sys.argv[1])
fc, pc = psd(sys.argv[2])
print(f"周波数分解能 {fb[1]-fb[0]:.3f} Hz")
print("帯域(Hz)      基準     現在     差(dB)")
for lo, hi in BANDS:
    b, c = band_power(fb, pb, lo, hi), band_power(fc, pc, lo, hi)
    print(f"{lo:>4}-{hi:<5} {b:9.4f} {c:9.4f} {10*np.log10(c/b):+7.1f}")
 
base_pk = peaks(fb, pb)
for fpk in peaks(fc, pc):
    if not np.any(np.abs(base_pk - fpk) <= 2 * (fb[1] - fb[0])):
        print(f"基準に無いピーク: {fpk:.1f} Hz")

引数の意味は、SciPyのリファレンスにある既定値と見比べると読みやすくなります。

• windowは既定が周期的なHann窓で、ここでは明示しています
• npersegは区間の長さです。窓が文字列のときの既定は256で、今回は4096にしています
• noverlapは未指定で、nperseg // 2の50%重なりになります。Hann窓なら50%が妥当な折衷だと説明されています
• detrendの既定は'constant'で、区間ごとに平均を引きます。サンプルに入れた0.5の直流オフセットはこれで消えます
• averageは'mean'と'median'を選べます(SciPy 1.2.0で追加)。medianは突発的なノイズの影響を受けにくい平均の取り方です
• scalingの既定は'density'で、単位は(入力の単位)の2乗をHzで割ったものになります

帯域のパワーは、PSDを帯域内で台形積分したものです。波形全体のRMSと見比べるときは、scaling='spectrum'の出力を使う方法もあります。

実行結果:どの帯域がどれだけ変わったか

2つのCSVを渡して走らせます。

python spectrum_compare.py baseline.csv current.csv

手元での出力は次のとおりでした。

周波数分解能 1.221 Hz
帯域(Hz)      基準     現在     差(dB)
   0-20       0.0004    0.0004    +0.0
  20-40       0.7183    0.7205    +0.0
  40-60       0.0653    0.5817    +9.5
  60-500      0.0070    0.0070    -0.0
 500-2000     0.0242    0.0687    +4.5
基準に無いピーク: 1180.4 Hz

読み取れることは3つです。

1. 回転1次を含む20〜40 Hzは変わらず、2次を含む40〜60 Hzが約9.5 dB(パワー比で約9倍)増えています
2. 60〜500 Hzはほぼ同じで、500〜2000 Hzだけ4.5 dB増えています
3. 基準期間のスペクトルに対応するピークが無い1180.4 Hzが、新しいピークとして挙がりました

帯域パワーの差(dB)は、広い範囲の底上げや総量の変化を見るのに向いています。一方、単独の狭いピークは広い帯域に埋もれやすいため、ピーク検出の出力と併用しています。find_peaksのprominenceは、ピークがまわりの谷からどれだけ突き出しているかを条件にする引数です。ここでは対数にしたPSDで10 dBを指定しました。

基準側のピークとの照合は、周波数分解能の2ビン分(約2.4 Hz)の許容幅で行っています。回転数が毎回ぴったり同じにならない設備では、この幅の置き方が結果を左右します。

npersegを既定のままにすると何が起きるか

区間の長さnpersegは周波数分解能を決めます。分解能はfs / npersegなので、5,000 Hzで4096点なら約1.22 Hzです。

SciPyの既定の256点にすると、同じcurrent.csvで分解能は19.53 Hzになりました。ピーク検出の結果を並べます。

256点では、25 Hzと50 Hzの回転成分がピークとして立たず、1171.9 Hzという、実際の1180 Hzから離れた周波数が返ります。「1次と2次の比が変わった」という、振動解析でいちばん見たい情報が落ちるわけです。

区間を長くすると平均できる区間の数は減ります。今回の60秒・4096点・50%重なりでは約145区間を平均しています。短い録波しか取れない設備では、分解能と平均の安定性の間でnpersegを選ぶことになります。

CLAUDE.mdとpermissionsで条件を固定する

スペクトル解析は、サンプリング周波数や基準ファイルの取り違えで、もっともらしい別の表が出ます。毎回口頭で指示するより、プロジェクトのCLAUDE.mdに書いておくほうが安定します。Claude Codeは、CLAUDE.mdをセッションの開始時に読み込みます。置き場所はプロジェクト直下の./CLAUDE.mdか./.claude/CLAUDE.mdです。

# 振動解析プロジェクト
 
## データの前提
- 加速度センサ、単位は g、サンプリング周波数は 5000 Hz(スクリプトの FS と一致させる)
- CSV の列は time_s, accel_g。基準期間は data/baseline/ に置く
 
## 作業ルール
- 生波形 CSV は cat や head で全文を表示しない。要約はスクリプトの出力だけを読む
- welch の nperseg を変えたら、周波数分解能(fs / nperseg)を出力に含める
- 帯域の境界と閾値は reference/bands.md の値を使い、勝手に変えない
- 結果の表には「異常」「故障」と書かず、基準との差(dB)と新規ピークだけを書く
- スクリプトを変更したら make_sample.py の合成データで既知の変化が検出されることを確かめる

CSVを全文表示させない規則は、トークン消費だけでなく、データの送信範囲を絞る意味もあります。ただしCLAUDE.mdは強制設定ではなく文脈として渡されるので、確実に止めたい操作にはフックなどの別の仕組みが要ります。送られる範囲の考え方はClaude Codeのデータは学習に使われるかで整理しています。

解析スクリプトの実行は、許可ルールで事前に通しておくと確認プロンプトが減ります。.claude/settings.jsonの例です。

{
  "permissions": {
    "allow": [
      "Bash(python spectrum_compare.py *)",
      "Bash(python make_sample.py)"
    ]
  }
}

ルールの書式はBash(npm run *)のようにコマンドにワイルドカードを使う形式で、deny、ask、allowの順に評価されます。解析スクリプト以外のコマンドは従来どおり確認が入ります。

Claude Codeへの依頼文と、検証の回し方

依頼は「作って」で終わらせず、確かめ方まで含めます。

data/baseline/ と data/current/ の振動CSVを Welch 法で比較する
spectrum_compare.py を作って。条件は CLAUDE.md に従うこと。
作ったら make_sample.py で合成データを2つ作り、
2次成分を増やした側で 40-60 Hz の差が出ること、
1180 Hz の新規ピークが出ることを実行して確かめて。

最後の一文が効きます。人工的に変化を入れたデータで、スクリプトがその変化を拾うか、拾った量が想定に近いかを実行で確かめさせるためです。たとえば窓関数やnpersegを直した後に、差が極端に小さくなったり、無関係な帯域が動いたりすれば、そこで検出できます。

反復の流れは次のとおりです。

実データに当てるときのつまずき

合成データではうまくいっても、実機では別の理由で差が出ます。

このうち回転数の違いは、帯域パワーの表だけでは見分けにくい問題です。回転数を別の列で持っているなら、基準と現在の回転数を表の先頭に並べて出させると、読む人が最初に気づけます。

表が示すのは「差」までで、診断ではない

Welch法の比較表は、基準期間と比べてどの周波数の成分がどれだけ変わったかを示します。その変化が軸受の損傷なのか、取付のゆるみなのか、単に負荷の違いなのかは、表からは決まりません。

そのため、表の見出しやClaude Codeの要約に「異常」「劣化」と書かせない運用が安全です。CLAUDE.mdに入れた「異常や故障と書かない」という一行は、この線引きのためのものです。警報の閾値を決めるのは、設備の保全基準と過去の故障履歴を持つ人になります。

寿命データから交換時期の目安を出す方法はClaude Codeでワイブル解析に分けてあります。振動の差は兆候の発見に、寿命分布は交換計画に、と役割が別です。

まとめ

振動CSVの解析は、SciPyのwelchにnpersegとaverageを明示して、基準期間との帯域別の差(dB)と新規ピークを出せば、Claude Codeに任せられる形になります。実行して確かめた合成データでは、増やした2次成分は約9.5 dB、足した1180 Hzの成分は新規ピークとして拾えました。

落とし穴は、既定の256点のままだと分解能が19.53 Hzまで粗くなり、25 Hzと50 Hzの成分が見えなくなる点です。CLAUDE.mdにサンプリング周波数と出力規則を書き、合成データで既知の変化を拾えるか確かめる習慣を付けると、スクリプトを直しても結果の信頼性が保てます。