« iPad miniの激安ケースを購入 | トップページ | 在宅ワークで大容量バッテリー「JVC-BN-RB10」を使ったら思わぬ落とし穴が »

2023年7月26日 (水)

MT法で音声の異常検知をやってみた

久しぶりに「機械学習」的なネタです。今回、MT法というのを使いました。

「MT法とは何ぞや?」という方は、以下を参照。

MT法で異常検知(python) - Qiita

こっちの動画もおすすめです。

要するにですね、正規分布に従っているデータからの外れっぷりを示す「マハラノビス距離」を使って、そのデータの異常度合を測るという方法で、以下のように使います。

(1) まず、正常データのみをそろえて、そのデータの平均値、分散を計算する。
(2) その平均、分散を使って、未知のデータのマハラノビス距離を計算する。

すると、異常なデータほど、このマハラノビス距離ってやつが大きくなるため、その数値を見て異常発生キターッ!と分かるというわけです。

この手法の特徴としては、正常データの分布を元に予測するため、大量の「正常データ」からごくわずかな「異常データ」の出現を測ることができます。
現実問題として、正常がたくさんで異常が滅多に出ないという製品が多いですから、そういう現場では重宝されます。

なので、まともそうな業者の保険修理件数があれば、ビッ〇モーターの保険修理の異常っぷりも検出できた・・・
いや、それは異常検知以前の問題ですね。

さて、今回はこのMT法を使って「音」の異常を検知してみます。

使った音は、うちの充電式シェーバーで髭を剃る時の音、です。

Img_3560

この何の変哲もない髭剃り機の何の異常を検知するのか?

この充電式のシェーバー、充電したてのころは元気がいいんですが、二週間近くなるとだんだんと電池が切れかかってきたのか、音が明らかに弱ってきます。といっても、毎日聞いていると気づかないのですが、充電した直後に明らかに音が変わるので、すぐに分かります。

つまり、徐々にではあるんですが日々、音が弱くなってきているんだろうな、と。

その衰えっぷりを、このマハラノビス距離ってやつで表現できないかと考えた次第。

てことで、今回用意したのは、充電した直後の1日目のデータ(hige01.wav)と、充電切れかけヤバヤバの12日目のデータ(hige12.wav)を使いました。どちらもだいたい50秒程度のデータ。

で、どうやってこの異常を検知しようか?

方針を決めるためにはまず、データを可視化してみましょう。

Hige01

まずは充電ほやほやの元気のいいシェーバーの音。
周波数は0~5000Hzに絞ってます。
だいたい1500~3000Hzに音の中心が集まっているのが分かります。

Hige12
一方こちらは充電切れ寸前のシェーバーの音データ。

明らかに音圧が減ってますが、ところどころ2000~3000Hzの音がガクッと下がってますね。

ひげを剃ってる間、負荷によってモーターの回転数が下がるみたいで、それがこのまだらな模様の原因となっているものと推測されます。

てことで、なんとなく分かったのは、「1500~3000Hz辺りの音データから特徴を取り出してやれば、シェーバーの充電時期が検知できるんじゃないか?」ということです。

ということで、以下のような方針を立てました。

(1) 1500~3000Hzの音を、250Hz刻みで平均化して出力
(2) その際、1秒ごとに音を分類できるように、IDをつけておく

上のスペクトログラムを各プログラムと合わせて、この(1)、(2)を実行するプログラムを以下に載せておきます。

import os
import sys
import numpy as np
import librosa
import matplotlib.pyplot as plt
import librosa.display
import pandas as pd
 
#音声ファイル読み込み
args = sys.argv
wav_filename = "hige01.wav"

data
, rate = librosa.load(wav_filename, sr=None, mono=False)
# フレーム長
fft_size = 2048
# フレームシフト長
hop_length = int(fft_size * 0.2)
# 時間の刻み幅 (s)
step_time = hop_length / rate
# 周波数の刻み幅 (Hz)
step_freq = rate / fft_size
# 短時間フーリエ変換実行
amplitude = np.abs(librosa.core.stft(data, n_fft=fft_size, hop_length=hop_length))
# 振幅をデシベル単位に変換
log_power = librosa.core.amplitude_to_db(amplitude)
fig_filename = os.path.splitext(os.path.basename(wav_filename))[0]
fft_spect = log_power.T
print(fft_spect.shape)
time_ax = []
freq_ax = []
for i in range(fft_spect.shape[0]):
    time_ax.append(step_time * i)
for j in range(fft_spect.shape[1]):
    freq_ax.append(step_freq * j)
time_ax = np.array(time_ax)
freq_ax = np.array(freq_ax)
df_spect = pd.DataFrame(fft_spect, columns=freq_ax)
df_time = pd.Series(time_ax,name="Time[s]")
df_spect_time = pd.concat([df_time,df_spect],axis=1)
df_spect_time.to_csv(fig_filename + '_pd.csv')
# グラフ表示
librosa.display.specshow(log_power, sr=rate, hop_length=hop_length, x_axis='time', y_axis='hz', cmap='jet')
plt.ylim(0,4000)
plt.xlim(0,10)
plt.colorbar(format='%+2.0f dB')  
plt.title('スペクトログラム' ,fontname="MS Gothic")
plt.savefig("10s_4000Hz_" + fig_filename)
#plt.show()
# データ加工
# 1500~3000Hzを250Hzづつ足して平均化し出力する
start_freq = 1500
interval_freq = 250
wave_feat = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat.append(tmp_wave)

start_freq1
= 1750
wave_feat1 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq1 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat1.append(tmp_wave)

start_freq2 = 2000
wave_feat2 = []

for
k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq2 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat2.append(tmp_wave)

start_freq3
= 2250
wave_feat3 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq2 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat3.append(tmp_wave)

start_freq4 = 2500
wave_feat4 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq4 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat4.append(tmp_wave)

start_freq5 = 2750
wave_feat5 = []

for k in range(fft_spect.shape[0]):
    tmp_wave = 0
    for l in range(int(interval_freq / step_freq)):
        tmp_wave += fft_spect[k, int(start_freq5 / step_freq) + l]
        tmp_wave = tmp_wave / (500 / step_freq)
    wave_feat5.append(tmp_wave)

time_id_step = 100 # IDをつけるためのID刻み幅
time_id = []

for
i in range(time_ax.shape[0]):
    time_id.append(i // time_id_step)
wave_feat_all = np.array([time_ax, wave_feat,wave_feat1,
         
wave_feat2,wave_feat3,wave_feat4,wave_feat5]).T

df_wave_feat = pd.DataFrame(wave_feat_all, columns=["time",start_freq,start_freq1,
                            start_freq2,start_freq3,start_freq4,start_freq5])
df_id = pd.Series(time_id,name="id")
df_wave_time = pd.concat([df_id,df_wave_feat],axis=1)
# IDで出力範囲指定する
df_wave_time = df_wave_time.query('0< id < 54')
df_wave_time.to_csv(fig_filename + '_feat.csv')
#
df_wave_feat.plot(subplots=True, sharex=True, figsize=(6,12))
plt.savefig("1500-3000Hz_wave_feat_" + fig_filename)

そういえばですが、予めnumpy、matplotlib、pandas、scikit-learn、scipy、seabornなどをpipコマンドでインストールしておいてください。
そしてこのプログラムでは「librosa」というライブラリを使ってます。
上のプログラム、ちょっとやっつけな部分がありますがご容赦ください。このプログラム名を「sound_spect.py」としておきます。
で、これをまず実行。

> python sound_spect.py

これを、10行目の「wave_filename = 」の後ろを「hige01.wav」「hige12.wav」と変えてそれぞれ実行すると、「hige01_feat.csv」「hige12_feat.csv」が出力されます。
なお、ファイルごとに録音時間がバラバラなので、このプログラムでは1~53秒までの音を抜き出すようにしております。

20230724-200718

こんな感じのファイルが作成されます。なお、列名が1500となっているのは、1500Hzから250Hz分の平均値、という意味です。また2列目のIDは、1~53まで変わります。

さて、この時系列データの特徴のとり方というのもいろいろな方法があるんですが、ここではオーソドックスに、

・ 1秒ごとに(つまり同じIDごとに)、1500~2750Hz帯それぞれの平均、最大、最小、標準偏差 を計算する

という方法を取ります。

で、それをやってくれるプログラムが、以下。

import pandas as pd
import numpy as np
def process_timeseries_data(input_file, output_file):
    # CSVファイルを読み込む
    df = pd.read_csv(input_file)
    # idごとにグループ化し、time、No.1、No.2、No.3の平均、最大、最小、標準偏差を計算
    grouped = df.groupby('id').agg({
        'time': 'first',   # 各グループ内の最初のtimeを取得
        '1500': ['mean', 'max', 'min', 'std'],  # No.1の平均、最大、最小、標準偏差を計算
        '1750': ['mean', 'max', 'min', 'std'],  # No.2の平均、最大、最小、標準偏差を計算
        '2000': ['mean', 'max', 'min', 'std'],   # No.3の平均、最大、最小、標準偏差を計算
        '2250': ['mean', 'max', 'min', 'std'],   # No.4の平均、最大、最小、標準偏差を計算
        '2500': ['mean', 'max', 'min', 'std'],   # No.5の平均、最大、最小、標準偏差を計算
        '2750': ['mean', 'max', 'min', 'std'],   # No.6の平均、最大、最小、標準偏差を計算
    })
    # 列名を整理
    grouped.columns = ['_'.join(col) for col in grouped.columns]
    # 出力ファイルに書き込む
    grouped.to_csv(output_file, index=True)
if __name__ == "__main__":
    input_file = "hige01_feat.csv"    # 入力ファイル名を適切に指定してください
    output_file = "output_hige01.csv"  # 出力ファイル名を適切に指定してください
    process_timeseries_data(input_file, output_file)

実はこのコード、考えるのがめんどくさくなって全部ChatGPTに作らせました。なんと、一発で動きましたよ。賢いでちゅね~、ChatGPT。
このプログラムは「spect_feat.py」と名付けます。
これを実行します。

> python spect_feat.py

なお、読み込むファイル名は下から3行目の「input_file =」の後ろ、出力ファイル名はその下の「output_file =」を変えます。

で、これを使って2つの音源ファイルから作った「hige01_feat.csv」「hige12_feat.csv」を読み込ませ、「output_hige01.csv」「output_hige12.csv」を得ます。

20230724-200735

なお、出力されたファイルの中身はこんな感じ。ID、時間に続いて、1500Hz帯の平均、最大、最小、標準偏差・・・という感じに並んでます。

それぞれの周波数帯に対して、4つの値。これが5つあって、かつ時間が53秒分あるため、全部で20×53のデータになります。
元が1025×6000とかいう巨大なデータなので、これでもかなり圧縮されております。

なお、「output_hige01.csv」のデータを「正常データ」、「output_hige12.csv」を「異常データ」とみなして、これをMT法によってより分けられるかを検証していきます。

まずは、この2つのデータをExcel等を使い「id」「time」の2列を消し、上から01、12の順に結合しておきます。ファイル名は、ここでは「hige01_12.csv」としておきます。

このファイルは、1行目に列名、次の53行は「正常データ」、その後ろの53行は「異常データ」ということになってます。

その上で、次のプログラムを実行します。

# 直交表から任意の因子数についての要因効果図を作る
# 使う直交表はL12~L64まで
# パッケージインポート
import os
import re
import collections
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
import sklearn
from sklearn import preprocessing
import scipy
from scipy.spatial import distance
import functools
import seaborn as sns
import sys
import tqdm as tq
from tqdm import tqdm
import gzip
import glob
import datetime as dt
import gc
import time
sns.set()

# =========== ユーザー定義 ==================
# データ読み込み
data0 = pd.read_csv('hige01_02.csv', encoding='cp932') # 51番目からは異常判定のあるデータ
# 正常データの行数 53行ある場合は、0始まりの52行目、という意味で「52」と入れる
s_data = 52
# ==========================================

# データを表示(表、グラフ)
print(data0.head())
data0.plot(figsize=(10, 6))
plt.show()

# マハラノビス距離計算前に、まずデータの正規化しておく(平均値0、標準偏差1のデータにする)
ss = preprocessing.StandardScaler()
data1 = ss.fit_transform(data0)
data = pd.DataFrame(data=data1, columns=data0.columns.values)
# 以下、要因効果図作成のため使用する直交表の一覧(2水準系)を定義
# L12 → 因子数が11以下
orthogonal_array_L12 = [
    [1,1,1,1,1,1,1,1,1,1,1],
    [2,1,2,1,1,1,2,2,2,1,2],
    [2,2,1,2,1,1,1,2,2,2,1],
    [1,2,2,1,2,1,1,1,2,2,2],
    [2,1,2,2,1,2,1,1,1,2,2],
    [2,2,1,2,2,1,2,1,1,1,2],
    [2,2,2,1,2,2,1,2,1,1,1],
    [1,2,2,2,1,2,2,1,2,1,1],
    [1,1,2,2,2,1,2,2,1,2,1],
    [1,1,1,2,2,2,1,2,2,1,2],
    [2,1,1,1,2,2,2,1,2,2,1],
    [1,2,1,1,1,2,2,2,1,2,2]
    ]
# L32 → 因子数が21以下
orthogonal_array_L32 = [
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2]
    ]
# L64 → 因子数が63以下
orthogonal_array_L64 = [
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,
1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,
2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1],
    [1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,
2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,
1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,
2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,
1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,
2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,
1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,
2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,1,1,
1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2],
    [1,1,1,2,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,1,1,1,1,2,2,2,2,2,2,
2,2,1,1,1,1,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,2,2,2,2,1,1,1,1],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,
2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,
1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,
2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,
1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,
2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,
1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,1,1,
2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2],
    [1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,
1,1,2,2,1,1,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,
2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,
1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,
2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,
1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,
2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,
1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,1,1,
2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1],
    [1,2,2,2,2,1,1,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,1,1,2,2,2,2,1,1,2,2,
1,1,1,1,2,2,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,2,2,1,1,1,1,2,2],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,
1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,
2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,
1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,
2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,
1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,
2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,
1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2],
    [2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,
2,1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,
1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,
2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,
1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,
2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,
1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,
2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,1,2,
1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1],
    [2,1,2,2,1,2,1,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,2,1,2,2,1,2,1,2,1,
2,1,1,2,1,2,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,2,1,2,1,1,2,1,2],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,
2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,
1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,
2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2],
  [2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,
1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,
2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,
1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,1,2,
2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1],
    [2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,
1,2,2,1,1,2,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,
2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,
1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,
2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,2,1,
1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,
2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,
1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,1,2,
2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2],
    [2,2,1,2,1,1,2,2,1,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,2,1,2,1,1,2,2,1,
1,2,1,2,2,1,1,2,2,1,2,1,1,2,1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1]
]

# MDの関数を定義
def get_mean_cov(feat):
    mean = np.mean(feat.astype(float).values, axis=0)
    cov = np.cov(feat.astype(float).values.T)
    return mean, cov
def get_score(feat, mean, cov):
    result = []
    # 分散共分散行列の逆行列を計算
    cov_i = np.linalg.pinv(cov)
    feat = feat.astype(float).values
    for i in range(len(feat)):
        result.append(distance.mahalanobis(feat[i], mean, cov_i))
    return result
# 正常と異常データを分ける
normal_feat = data.iloc[:s_data,:]
anomaly_feat = data.iloc[s_data+1:,:]
print(normal_feat.shape, anomaly_feat.shape)
# 平均、分散を求める
mean, cov = get_mean_cov(normal_feat)
# マハラノビス距離の計算(正常、異常それぞれ)
normal_score = get_score(normal_feat, mean, cov)
anomaly_score = get_score(anomaly_feat, mean, cov)
# 正常、異常のマハラノビス距離をプロット
plt.plot(normal_score)
plt.plot(anomaly_score)
plt.show()

md_ot = []

num = data.shape[1] # データから因子数を取得

# 因子数ごとに使う直交表を選ぶ (最大因子数は127まで)
if num<12:
    orthogonal_array = orthogonal_array_L12
elif num<31:
    orthogonal_array = orthogonal_array_L32
elif num<63:
    orthogonal_array = orthogonal_array_L64
else:
    print("因子数は最大63まで!")
    sys.exit()
k = 0
for c in orthogonal_array:
    b = []
    for i in range(num):
        if c[i]==2:
            b.append(i) # 直交表で「2」の因子を取り除く為のインデックスの配列を作成
    data_orth = data.drop(columns=data.columns[b]) # 直交表で”2”の項目を取り除く
    normal_feat = data_orth.iloc[:s_data,:]     # 正常データ
    anomaly_feat = data_orth.iloc[s_data+1:,:]  # 異常データ
    # 直交表によって選択された因子のみでマハラノビス計算用の平均、分散を計算
    mean, cov = get_mean_cov(normal_feat)
   
    # 上の平均、分散を用いて、すべての異常値についてのマハラノビス距離を計算
    anomaly_score = get_score(anomaly_feat, mean, cov)
    # 異常値のマハラノビス距離の平均値を結合
    md_ot.append(sum(anomaly_score)/len(anomaly_score))
sum_md = np.zeros((2,num))
count_md_p = np.zeros(num)
count_md_n = np.zeros(num)
sum_msd = np.zeros((2,num))
# マハラノビス距離計算後の処理:因子ありなしの差を可視化する
k=0
for c in orthogonal_array:
    for i in range(num):
        if c[i]==1:
            sum_md[0,i] += md_ot[k]
            sum_msd[0,i] += 1 / (md_ot[k])**2
            count_md_p[i] +=1
           
        else:
            sum_md[1,i] += md_ot[k]
            sum_msd[1,i] += 1 / (md_ot[k])**2
            count_md_n[i] +=1
    k +=1
# それぞれの因子のあり(posi)、なし(nega)のMD平均値を格納する
posi_mean_md = np.zeros(num)
nega_mean_md = np.zeros(num)
for i in range(num):
    posi_mean_md[i] = sum_md[0,i] / count_md_p[i]
    nega_mean_md[i] = sum_md[1,i] / count_md_n[i]
# 各々の因子有無のマハラノビス距離
print(posi_mean_md)
print(nega_mean_md)
# マハラノビス距離についての要因効果図を頑張って書く
# Y軸の範囲を固定
ymin = min([np.min(posi_mean_md), np.min(nega_mean_md)])*0.95
ymax = max([np.max(posi_mean_md), np.max(nega_mean_md)])*1.05
t = np.linspace(0, 1, 2)
y =  np.zeros((num,2))
fig , axes = plt.subplots(1 , num , figsize=(10,5))
for i in range(num):
    y[i] = [posi_mean_md[i], nega_mean_md[i]]
    axes[i].set_title(data.columns.values[i], fontsize=7)
    axes[i].set_xlim(-0.3, 1.3)
    axes[i].set_ylim(ymin, ymax)
    axes[i].plot(t, y[i], color="black", marker='.',markersize=10, label=str(i+1))
    if i>0:
        axes[i].tick_params(labelleft=False, labelright=False,
labeltop=False, length=0, labelsize=7)
    else:
        axes[i].tick_params(length=0, labelsize=7)
plt.show()

これが、今回の肝となるコードになります。プログラム名を「MT_method.py」としておきます。
これを、例によって実行します。

> python MT_method.py

このプログラムがやっていることは、以下の3つ。

(1) まず最初の53行の「正常データ」を使って、平均、分散を計算する
(2) 続いて、この「正常データ」とその下53行の「異常データ」それぞれのマハラノビス距離を算出しグラフ化
(3)たくさんの因子を使って計算したのですが、その因子の寄与度を求めるために、正常データを使ってその因子の有り無しの平均、分散を計算し、異常値のマハラノビス距離を計算、その距離値がどれくらい変動するかを「0:あり」「1:なし」で比較する

ということをやってます。

特にこの(3)が肝で、プログラムの途中に見える1と2だけの強烈な配列(直交表)は、この寄与度を算出するために入れております。
このプログラム自体は汎用性があって、最大63因子までのMT法と、その寄与度分析が可能となっております。
先日の東京出張は、まさにこれを作るために行ったようなものです。

では、結果を見ていきましょう。まずは(2)のマハラノビス距離の比較です。

Hige_graph02

こんなのが出ました。青線が正常、オレンジ色が異常データのものです。
若干ラップしてますが、大体5のところに線を引けば、概ね分けられます。

一般にMT法では、マハラノビス距離が4で分ける、というのが定番ですけど、正常データで4を超えているところが多過ぎて・・・あまりいいデータじゃないってことですかね。

で、このマハラノビス距離の寄与が大きい因子を見るため、(3)を実行した結果が以下。

Hige_graph03

ここでは「要因効果図」と呼びますが、この図では後ろから4つ目、5つ目が大きく変動してますね。これはそれぞれ、「2500_std」「2750_mean」、つまり「2500~2750Hzの標準偏差」と「2750~3000Hzの平均値」がマハラノビス距離に対して効いていることを示しています。

Hige_graph04

試しにこの二つを縦横軸にしてプロットしてみると、御覧の通り、この2軸だけでまあまあ正常/異常が分離できてます。
他の因子を縦横にしても、この青とオレンジが結構入り乱れてました。

最初の可視化の時点で、3000Hz付近の音の揺れの違いがあるとみていたので、これはまあ納得の結果、といったところでしょうか。

という具合に、音の異常を身近なもので実験してみました。

もっとも、たまたま分離できただけかもしれません。本来ならもっとたくさんのデータを使ったほうが精度が高いと思われます。

正直言いますと、1日目、12日目だけでなく、2日目、3日目・・・のデータもあるんですが、実はどれも似たようなマハラノビス距離の分離っぷりになってしまいました。1日目から見れば、どいつもこいつも「異常」ということに。
なので、シェーバーの音自体というよりも、録音環境(マイクの位置、髭の剃り方の違い、など)の方が効いているのかもしれません。

また、時系列データの特徴量抽出にはよく「tsfresh」というライブラリが使われます。今回は時系列データの特徴量としては基本の平均、最大、最小、標準偏差だけを使いましたが、tsfreshでは一つの時系列データから779種類の特徴量を作り出してくれます。そちらを使えば、もっといい特徴量が見つかったかも。

まあ、そういうわけで今回の話は、ご参考ということで。


よくわかるMTシステム―品質工学によるパターン認識の新技術

« iPad miniの激安ケースを購入 | トップページ | 在宅ワークで大容量バッテリー「JVC-BN-RB10」を使ったら思わぬ落とし穴が »

数値解析系」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« iPad miniの激安ケースを購入 | トップページ | 在宅ワークで大容量バッテリー「JVC-BN-RB10」を使ったら思わぬ落とし穴が »

無料ブログはココログ

スポンサード リンク

ブログ村