Zenn
🦆

振動解析の備忘録

に公開
numpy
FFT
振動
tech

はじめに

仕事で振動解析に触れる機会があったので、分析方法に関する備忘録を残しました。
一般的な手法を学べたと思うので、ご参考になれば幸いです。

数式とコードで内容を整理していこうと思います。

修正履歴

  • 2025-05-12: 関数の誤りを修正した。
  • 2025-05-14: librosa.load()がv1.0で廃止されるとのメッセージが出てきたので、pydubによる読み取りに変更した。
    • メルフィルタバンクと呼び周波数を作成する関数の引数を変更した。

使用するライブラリのインポート

ライブラリをimportする
# ---- ライブラリをimportする ----

# ---- DataFrameライブラリをimportする ----
import polars as pl
import polars.selectors as cs
#import pandas as pd
import numpy as np

#  ---- DataFrameの設定を初期化する ----
pl.Config(fmt_str_lengths = 100, tbl_cols = 100, tbl_rows = 10)
np.set_printoptions(precision = 3, floatmode = 'fixed')
#pd.set_option('display.max_columns', 100)
#pd.set_option('display.max_rows', 100)

# ---- プロットライブラリをimportする ----
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

# ※日本語でプロットしたい場合、英語フォントはコメントアウトする
import japanize_matplotlib
#plt.rcParams['font.family'] = 'Times New Roman'

# ---- プロットの設定を初期化する ----
plt.rcParams['mathtext.fontset'] = 'stix'
plt.rcParams["font.size"] = 15
plt.rcParams['xtick.labelsize'] = 10
plt.rcParams['ytick.labelsize'] = 10
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['axes.linewidth'] = 1.0
plt.rcParams['axes.grid'] = True

# ---- 波形分析ライブラリをimportする ----
# 音声データを読み込む
from pydub import AudioSegment
from pydub.utils import mediainfo

# 音声データを前処理する
import librosa
import librosa.display

# ---- 科学計算ライブラリをimportする ----
from scipy.signal.windows import hamming# ハミング窓関数

時間波形

音や振動など実世界の時間波形 y(t) は時間的に連続でその値も連続なアナログ値である。
波動をコンピュータで解析するためには時間と値の両方を離散的な値に変換する必要がある。
一定間隔でデータを記録することを標本化と呼び、波動の値を浮動小数点数(float等)で近似することを量子化と呼ぶ。
標本化は時間の離散化で、量子化は波動の値の離散化である。

サンプリング周波数 sr (Hz) で T 秒間記録した波動の値を配列 y[n] に保存する。
n が離散的な時間のインデックスで、 y[n] が離散的な波動の値である。
このとき、標本化された時間波形 y[n] と時刻 t[n] との関係は次のようになる。

\begin{align*} y[n] &\approx y(t[n])\\ t[n] &= \frac{n}{N}T\\ N &\equiv \lfloor \text{sr} \times T \rfloor\\ n &= 0, 1, \dots , N - 1 \end{align*}

保存できる周波数の上限は \text{sr} / 2 (Hz)であり、これは2点のデータの繰り返しまでを保存できることを意味している(標本化定理)。
PCの性能が向上したため、分析において標本化定理が課題になることはほとんど無いと感じる。
組込プログラミングのような性能に制限がある場面では標本化定理を考慮する必要がある。

データ数 N はFFTのために2の累乗でなければならないが、Numpyの numpy.fft.fft で処理できてしまうので偶数であればよい。
Arduino等の組込の場合は2の累乗を守らねばならない。

コード

y_{l}(t) = \cos {2 \pi f_{l} t}, (l = 1, 2) の波形の和を分析することにする。
例えば、オンライントーンジェネレーターで次のコードの2種類の正弦波を流すと、うなりがする音を再生できる。
この音を分析していくことにする。

時間波形をプロットするとうなりを視覚的に確認できる。

時間波形をプロットする
# 1. 記録の設定
# サンプリングレート
sr = 44100# Hz

# 記録時間と開始・終了時刻
start_time = 0# 秒
T = 1# 秒
end_time = start_time + T# 秒

# データ数
N = int(sr * T)# 個

# 時刻
t = np.linspace(start = start_time, stop = end_time, num = N, endpoint = True, dtype = float)


# 2. 波形の合成
# 基本周波数
f1 = 405# Hz
f2 = 400# Hz

# 基本波形
y1 = np.cos(2 * np.pi * f1 * t )
y2 = np.cos(2 * np.pi * f2 * t )

# 合成波形
yt = (y1 + y2) / 2


# 3. プロット
# Figure
fig = plt.figure(figsize = (12, 11))# (横, 縦)

# 1. 元の波形をプロットする
ax_original = fig.add_subplot(2, 2, 1)
sns.lineplot(x = t, y = yt, ax = ax_original, color = "blue", linewidth = 0.1)
plt.setp(ax_original, title = "Time domain Waveform", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_original, xlim = (start_time, end_time), xticks = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])


# 2. うなり1個分の波形をトリミングする
ax_trim = fig.add_subplot(2, 2, 2)
sns.lineplot(x = t, y = yt, ax = ax_trim, color = "blue", linewidth = 0.2)
plt.setp(ax_trim, title = "Trimmed Time domain Waveform", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_trim, xlim = (0.1, 0.3), xticks = [0.1, 0.2, 0.3], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])


# レイアウト調整
plt.tight_layout()
plt.show()

周波数スペクトル(離散フーリエ変換)

計測された時間波形の周波数成分の強度や位相を調べることを、スペクトル解析と呼ぶ。
よく知られている手法がフーリエ変換であり、時間波形 y[n] を様々な周波数スペクトル Y[k] を振幅とする三角関数の和で表す手法である。
標本化された波形のフーリエ変換を離散フーリエ変換と呼び、定義式や周波数との関係は次のようになる。

\begin{align*} Y[k] &\equiv \sum_{n = 0}^{N - 1} x[n] \exp{(-j\frac{2 \pi nk}{N})}\\ f[k] &= \left \{ \begin{matrix} \frac{\text{sr}}{N}k & (0 \leq k \leq N/2)\\ - \frac{\text{sr}}{N}(N - k) & (N/2 + 1 \leq k \leq N - 1)\\ \end{matrix} \right.\\ k &= 0, 1, \dots , N / 2 \end{align*}

1行目の定義式を簡略化して Y[k] = \hat{F}[x[n]] と書くことが多い。
2行目より周波数軸の分解能は \frac{\text{sr}}{N} であり、ピーク位置を精密に分析するには計測時間を長くする必要があることがわかる。

周波数スペクトルから時間波形への逆変換を逆フーリエ変換と呼ぶ。

\begin{align*} y[n] &= \sum_{k = 0}^{N - 1} Y[k] \exp{(j\frac{2 \pi nk}{N})}\\ &= \hat{F}^{-1}[Y[k]] \end{align*}

振幅スペクトル、パワースペクトル、対数パワースペクトル

周波数スペクトルの絶対値 |Y[k]| を振幅スペクトルと呼ぶ。
周波数スペクトルの絶対値の2乗 P = |Y[k]|^{2} をパワースペクトルと呼ぶ。
どちらも周波数ごとの音の強さを表す。

パワースペクトルの値域が広いことことから、対数パワースペクトル(dB)を使用する測定器メーカーが多いように感じる。
小野測器では明示的に区別するために |Y[k]|^{2} をリニアスペクトルと呼んでいる。
目的に応じて測定時間が異なるとスペクトルの値が変化してしまうため、対数パワースペクトルを標準化しておくのが便利である。
真数の分子の \epsilon は対数を発散させないための十分に小さい正の数、分母の (Y^{2})_{sum} は規格化定数である。

\begin{align*} \log_{10}P &\equiv 10 \log_{10} \frac{|Y[k]|^{2} + \epsilon}{ (Y^{2})_{sum} }\\ (Y^{2})_{sum} &\equiv \sum_{k = 0}^{N - 1}|Y[k]|^{2}\\ 0 &< \epsilon << |Y[k]|^{2} \end{align*}

周波数スペクトルの数学的性質

インデックス [1, N / 2] は正の周波数、 インデックス [N / 2 + 1, N - 1] は負の周波数の周波数スペクトルである。
この性質は解析信号を計算するときに利用する。

この性質を確かめる。
k = 1, 2, \dots, N / 2 の範囲で X[N - k] を計算すると複素共役 X[k]^{*} になる。

\begin{align*} Y[N - k] &= \sum_{n = 0}^{N - 1} y[n] \exp{(-j\frac{2 \pi n(N - k)}{N})}\\ &= \sum_{n = 0}^{N - 1} y[n] \exp{(-j(2 \pi n) +j\frac{2 \pi nk}{N})}\\ &= \sum_{n = 0}^{N - 1} y[n] \exp{(+j\frac{2 \pi nk}{N})}\\ &= X[k]^{*} \end{align*}

次に、定義域を無視して負の周波数のスペクトルを計算する。
同様に複素共役を得られるので、この性質を確認できた。

\begin{align*} Y[- k] &= \sum_{n = 0}^{N - 1} y[n] \exp{(-j\frac{2 \pi n(-k)}{N})}\\ &= \sum_{n = 0}^{N - 1} y[n] \exp{(+j\frac{2 \pi nk}{N})}\\ &= Y[k]^{*}\\ \rightarrow Y[N - k] &= Y[- k] \end{align*}

本節の議論は連続関数のフーリエ変換の数学的性質の議論を省略したために奇妙な展開になっている。

コード

テスト用の時間波形とトリミングした時間波形のパワースペクトル、対数パワースペクトルを計算する。
規格化無しのパワースペクトルの値が1桁以上異なるためそのままでは比較できない。
計測時間によらず周波数スペクトルの強度を比較する場合、規格化が必要であることがわかる。
計測時間を短くすると周波数軸の分解能が落ちることがコードでも確認できた。

うなりの周波数スペクトル(5Hz)はFFTでは分析できなかった。

時間波形の周波数スペクトルと周波数配列を計算する
# ---- 時間波形の周波数スペクトルと周波数配列を計算する ----
# yt: 時間波形
# sr: サンプリング周波数. 書籍によってはfsと書かれていることもある
def fft(yt, sr):
    # FFT
    yf = np.fft.fft(yt)
    # データ数
    N = len(yt)
    # 周波数
    f = np.zeros(N)
    f[1:(N // 2 + 1)] = np.linspace(start = sr / N, stop = sr / 2, num = N // 2, endpoint = True)# 正の周波数
    f[(N // 2 + 1):] = np.linspace(start = - sr / N * (N - 1), stop = -sr / N, num = N // 2 - 1, endpoint = True)# 負の周波数
    # 周波数スペクトルと周波数軸を返す
    return yf, f
振幅スペクトル、パワースペクトル、対数パワースペクトルを計算する
# 振幅スペクトル、パワースペクトル、対数パワースペクトルを計算する
# yf: 周波数スペクトル
def spectrum(yf):
    # 対数を発散させない定数
    epsilon = 1e-20
    #振幅スペクトル
    amplitude_spectrum = np.abs(yf)
    # パワースペクトル
    power_spectrum = amplitude_spectrum ** 2
    # 正規化定数
    normalization_factor = np.sum(power_spectrum)
    # 正規化された対数パワースペクトル
    log_power_spectrum = 10 * np.log10((power_spectrum + epsilon) / normalization_factor)
    # 振幅スペクトル、パワースペクトル、対数パワースペクトルを返す
    return amplitude_spectrum, power_spectrum, log_power_spectrum
時間波形のパワースペクトルや対数パワースペクトルをプロットする
# 1. 時間波形をトリミングする
t_short = t[(t >= 0.1) & (t <= 0.3)]
yt_short = yt[(t >= 0.1) & (t <= 0.3)]
N_short = len(yt_short)

# 2. 周波数スペクトルを計算する
# 元の波形
yf, f = fft(yt = yt, sr = sr)
# トリミングした波形
yf_short, f_short = fft(yt = yt_short, sr = sr)

# 3. パワースペクトル、対数パワースペクトルを計算する
# 元の波形
amplitude_spectrum, power_spectrum, log_power_spectrum = spectrum(yf)
# トリミングした波形
amplitude_spectrum_short, power_spectrum_short, log_power_spectrum_short = spectrum(yf_short)


# 4. プロット
# Figure
fig = plt.figure(figsize = (12, 16))# (横, 縦)

# 1. 元の波形をプロットする
ax_original = fig.add_subplot(3, 2, 1)
sns.lineplot(x = t, y = yt, ax = ax_original, color = "blue", linewidth = 0.1)
plt.setp(ax_original, title = "Time domain Waveform", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_original, xlim = (start_time, end_time), xticks = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2. うなり1個分の波形をトリミングする
ax_trim = fig.add_subplot(3, 2, 2)
sns.lineplot(x = t, y = yt, ax = ax_trim, color = "blue", linewidth = 0.2)
plt.setp(ax_trim, title = "Trimmed Time domain Waveform", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_trim, xlim = (0.1, 0.3), xticks = [0.1, 0.2, 0.3], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3. 元の波形のパワースペクトル
ax_power_spectrum = fig.add_subplot(3, 2, 3)
sns.lineplot(x = f[:(N // 2 + 1)], y = power_spectrum[:(N // 2 + 1)], ax = ax_power_spectrum, color = "blue", linewidth = 1)
plt.setp(ax_power_spectrum, title = "Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Power Spectrum")
plt.setp(ax_power_spectrum, xlim = (0, 500))

# 4. トリミングした波形のパワースペクトル
ax_trim_power_spectrum = fig.add_subplot(3, 2, 4)
sns.lineplot(x = f_short[:(N_short // 2 + 1)], y = power_spectrum_short[:(N_short // 2 + 1)], ax = ax_trim_power_spectrum, color = "blue", linewidth = 1)
plt.setp(ax_trim_power_spectrum, title = "Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Power Spectrum")
plt.setp(ax_trim_power_spectrum, xlim = (0, 500))

# 5. 元の波形の標準化した対数パワースペクトル
ax_log_power_spectrum = fig.add_subplot(3, 2, 5)
sns.lineplot(x = f[:(N // 2 + 1)], y = log_power_spectrum[:(N // 2 + 1)], ax = ax_log_power_spectrum, color = "blue", linewidth = 1)
plt.setp(ax_log_power_spectrum, title = "Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_log_power_spectrum, xlim = (0, 500), ylim = (-150, 0))

# 6. トリミングした波形の標準化した対数パワースペクトル
ax_trim_log_power_spectrum = fig.add_subplot(3, 2, 6)
sns.lineplot(x = f_short[:(N_short // 2 + 1)], y = log_power_spectrum_short[:(N_short // 2 + 1)], ax = ax_trim_log_power_spectrum, color = "blue", linewidth = 1)
plt.setp(ax_trim_log_power_spectrum, title = "Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_trim_log_power_spectrum, xlim = (0, 500), ylim = (-150, 0))

# レイアウト調整
plt.tight_layout()
plt.show()

解析信号

この節のモチベーションは、うなりの周波数スペクトルを計算できるようになることである。
三角関数の和の公式を用いることで、周期が近い2つの波動を足し合わせるとうなりが発生することを確認することができる。

\begin{align*} y(t) &= \cos (2 \pi f_{1} t) + \cos (2 \pi f_{2} t)\\ &= 2 \cos(2 \pi f_{mean}t) \cos (2 \pi f_{\Delta}t)\\ &\text{Where} \ f_{mean} \equiv \frac{f_{1} + f_{2}}{2}, \ f_{\Delta} \equiv \frac{f_{1} - f_{2}}{2}, f_{1} > f_{2} \end{align*}

フーリエ変換は時間波形 y(t) を様々な周波数スペクトル X(f) を振幅とする三角関数の和で表す手法のため、そのままでは時間波形の掛け算を表現できない。
先人の素晴らしいアイデアを採用する。
元の時間波形を実部にとり、位相を -\pi/2 ずらした時間波形 \hat{y}(t) を虚部にとる複素数波形 z(t) を計算する。
複素数波形 z(t) の振幅がうなりの包絡線の変化を、z(t) の位相が合成波の位相を表す式になる。\

\begin{align*} z(t) &\equiv y(t) + j \hat{y}(t)\\ &= \left( \cos(2 \pi f_{1} t) + \cos(2 \pi f_{2} t) \right) + j \left( \cos(2 \pi f_{1} t - \frac{\pi}{2}) + \cos(2 \pi f_{2} t - \frac{\pi}{2}) \right)\\ &= 2 \cos(2 \pi f_{\Delta} t) \cos(2 \pi f_{mean} t) + j \left( 2 \cos(2 \pi f_{\Delta} t) \sin(2 \pi f_{mean} t) \right)\\ &= 2 \cos(2 \pi f_{\Delta} t) \left( \cos(2 \pi f_{mean} t) + j \sin(2 \pi f_{mean} t) \right)\\ &= 2 \cos(2 \pi f_{\Delta} t) \exp j(2 \pi f_{mean} t) \end{align*}

うなりを表現できる複素数波形 z(t) を解析信号と呼ぶ。
解析信号の絶対値 |z(t)| を瞬時振幅と呼ぶ。
解析信号の位相 \theta(t) = \tan^{-1} \frac{\hat{y}(t)}{y(t)} を瞬時位相と呼ぶ。
位相を -\pi/2 ずらす変換のことをヒルベルト変換と呼ぶ。
ヒルベルト変換の定義式を含めて、解析信号の定義式を下記に記す。

\begin{align*} z(t) &\equiv y(t) + j \hat{y}(t)\\ \hat{y}(t) &\equiv \hat{H}[y(t)]\\ &= \hat{F}^{-1}[-j \text{sgn}(f) Y(f) ] \end{align*}

解析信号の計算方法

フーリエ変換を用いることで解析信号を簡単に計算できる。
解析信号の周波数スペクトル Z(f) を計算する。

\begin{align*} Z(f) &= Y(f) + i \hat{Y}(f)\\ &= Y(f) + i(-j \text{sgn}(f) Y(f))\\ &= \left \{ \begin{matrix} 2 Y(f) & (f > 0, \text{sgn} = 1)\\ Y(f) & (f = 0, \text{sgn} = 0)\\ 0 & (f > 0, \text{sgn} = -1)\\ \end{matrix} \right. \end{align*}

解析信号を計算する手順は次の通りである。

  1. 信号をフーリエ変換する
  2. 負の周波数成分を0にする
  3. 正の周波数成分の値を2倍にする
  4. 逆フーリエ変換する

コード

解析信号の瞬時振幅の対数パワースペクトルをプロットした。
うなりの周波数ピークを表現できることを確認した。

解析信号の包絡線と包絡線のパワースペクトル、対数パワースペクトルを計算する
# ---- 解析信号の包絡線と包絡線のパワースペクトル、対数パワースペクトルを計算する ----
# yt: 時間波形
def analytic_signal(yt):
    # 対数を発散させない定数
    epsilon = 1e-20
    # データ数
    N = len(yt)
    # 解析信号の周波数スペクトル
    zf = np.fft.fft(yt)
    zf[1:] = 2.0 * zf[1:]# zf = 2 * xf, f > 0
    zf[(N // 2 + 1):] = 0.0# zf = 0, f < 0
    # 解析信号
    zt = np.fft.ifft(zf)
    # 解析信号の包絡線
    envelope = np.abs(zt)
    # 包絡線のパワースペクトル
    fluctuation_spectrum = np.abs( np.fft.fft( envelope ) ) ** 2
    # 正規化定数
    normalization_factor = np.sum(fluctuation_spectrum)
    # 包絡線の対数パワースペクトル
    log_fluctuation_spectrum = 10 * np.log10( (fluctuation_spectrum + epsilon) / normalization_factor) 
    # 解析信号の包絡線と包絡線のパワースペクトル、対数パワースペクトルを返す
    return envelope, fluctuation_spectrum, log_fluctuation_spectrum
解析信号の包絡線などをプロットする
# 1. 解析信号を計算する
envelope, fluctuation_spectrum, log_fluctuation_spectrum = analytic_signal(yt)


# 2. プロット
# Figure
fig = plt.figure(figsize = (12, 16))# (横, 縦)

# 1. 元の波形をプロットする
ax_original = fig.add_subplot(3, 2, 1)
sns.lineplot(x = t, y = yt, ax = ax_original, color = "blue", linewidth = 0.1)
plt.setp(ax_original, title = "Time domain Waveform", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_original, xlim = (start_time, end_time), xticks = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 5. 元の波形の標準化した対数パワースペクトル
ax_power_spectrum_log = fig.add_subplot(3, 2, 2)
sns.lineplot(x = f[:(N // 2 + 1)], y = log_power_spectrum[:(N // 2 + 1)], ax = ax_power_spectrum_log, color = "blue", linewidth = 1)
plt.setp(ax_power_spectrum_log, title = "Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_log, xlim = (0, 500), ylim = (-150, 0))

# 3. 解析信号の元の波形をプロットする
ax_envelope = fig.add_subplot(3, 2, 3)
sns.lineplot(x = t, y = envelope, ax = ax_envelope, color = "blue", linewidth = 1)
plt.setp(ax_envelope, title = "Envelope of Time domain Waveform", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope, xlim = (start_time, end_time), xticks = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 5. 解析信号の標準化した対数パワースペクトル
ax_fluctuation_spectrum_db = fig.add_subplot(3, 2, 4)
sns.lineplot(x = f[:(N // 2 + 1)], y = log_fluctuation_spectrum[:(N // 2 + 1)], ax = ax_fluctuation_spectrum_db, color = "blue", linewidth = 1)
plt.setp(ax_fluctuation_spectrum_db, title = "Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db, xlim = (0, 500), xticks = [0, 5, 100, 200, 300, 400, 500], ylim = (-150, 0))


# レイアウト調整
plt.tight_layout()
plt.show()

解析にフリーBGMを使用する

解析信号などの処理を数式で確認できたため、前節まで正弦波の和を使用してきた。
この後の処理は複雑で解析的に確認できないため、フリーBGMを使って解析する。

AI学習可能なサイトフリー音源・無料BGM「BGMer(ビージーエマー)」から次の2つの楽曲をダウンロードした。
ロジスティック回帰ですらAI(機械学習)の一手法であるため、素材を使わせていただく際にはきちんと規約を確認したい。
それぞれの楽曲情報は参考URLを参照のこと。

明るい雰囲気の「はじめの一歩」と不穏な雰囲気の「ファイルX」の時間波形や対数パワースペクトル、解析信号をプロットする。
不穏な雰囲気の曲には音のうなりを持たせているのかもしれない。
低周波の対数パワースペクトルにピークがあるため、解析信号で不穏な雰囲気の曲のうなりの情報を取得できていることがわかる。

音楽データを読み込んで時間波形に変換する
# ---- 音声データをnumpy配列に変換する ----
# file_path: ファイルパス
# file_extension: 拡張子
def from_audio_to_ndarray(file_path: str, file_extension: str):
    # 音声データを読み込む
    audio = AudioSegment.from_file(file = file_path, format = file_extension)
    # Numpy配列に変換する
    yt = np.array( audio.get_array_of_samples() )
    # 振幅を[-1, 1]で正規化する
    yt = yt / (2 ** (audio.sample_width * 8 - 1) - 1)
    # ステレオの場合、モノラルに変換する
    if audio.channels == 2:
        yt = np.mean(yt.reshape((-1, 2)), axis = 1)
    # サンプリングレートを取得する
    sr = audio.frame_rate
    # 時刻データ
    t = np.linspace(start = 0, stop = len(yt) / sr, num = len(yt), endpoint = True)
    # 時間波形、時刻、サンプリング周波数を返す
    return yt, t, sr

# ---- 音楽データを読み込んで時間波形に変換する ----
# librosa.load()に警告文が出てきたのでpydubに切り替える。
# path: 音楽データのパス
# sr: サンプリング周波数
def load_audio(path, sr):
    # データを読み込む
    yt, sr = librosa.load(path, sr = sr)
    # 時刻データ
    t = np.linspace(start = 0, stop = len(yt) / sr, num = len(yt), endpoint = True)
    # 時間波形、時刻、サンプリング周波数を返す
    return yt, t, sr
読み込んだ音楽データの時間波形やパワースペクトルをプロットする
# 1. はじめの一歩についてデータの読み込みから解析信号の計算までを行う
# 1.1. データの読み込み(サンプリング周波数48kHzのデータを読み込む)
yt_1s, t_1s, sr_1s = from_audio_to_ndarray(file_path = "bgm/BGMer はじめの一歩.mp3", file_extension = "mp3")

# 1.2. 冒頭の10秒間を抽出する
yt_1s = yt_1s[t_1s <= 10]
t_1s = t_1s[t_1s <= 10]

# 1.3 記録の設定
# データ数
N_1s = len(yt_1s)

# 1.4. 周波数スペクトルを計算する
yf_1s, f_1s = fft(yt = yt_1s, sr = sr_1s)

# 1.5. パワースペクトル、対数パワースペクトルを計算する
amplitude_spectrum_1s, power_spectrum_1s, log_power_spectrum_1s = spectrum(yf = yf_1s)

# 1.6. 解析信号を計算する
envelope_1s, fluctuation_spectrum_1s, log_fluctuation_spectrum_1s = analytic_signal(yt = yt_1s)


# 2. ファイルXについてデータの読み込みから解析信号の計算までを行う
# 2.1. データの読み込み(サンプリング周波数48kHzのデータを読み込む)
yt_X, t_X, sr_X = from_audio_to_ndarray(file_path = "bgm/BGMer ファイルX.mp3", file_extension = "mp3")

# 2.2. 冒頭の10秒間を抽出する
yt_X = yt_X[t_X <= 10]
t_X = t_X[t_X <= 10]

# 2.3 記録の設定
# データ数
N_X = len(yt_X)

# 2.4. 周波数スペクトルを計算する
yf_X, f_X = fft(yt = yt_X, sr = sr_X)

# 2.5. パワースペクトル、対数パワースペクトルを計算する
amplitude_spectrum_X, power_spectrum_X, log_power_spectrum_X = spectrum(yf = yf_X)

# 2.6. 解析信号を計算する
envelope_X, fluctuation_spectrum_X, log_fluctuation_spectrum_X = analytic_signal(yt = yt_X)


# 3. プロット
# Figure
fig = plt.figure(figsize = (6 * 4, 1 + 5 * 2))# (横, 縦)

# 3.1. はじめの一歩
# 3.1.1. 時間波形
ax_1s = fig.add_subplot(2, 4, 1)
sns.lineplot(x = t_1s, y = yt_1s, ax = ax_1s, color = "blue", linewidth = 0.1)
plt.setp(ax_1s, title = "はじめの一歩", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_1s, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3.2.2 標準化した対数パワースペクトル
ax_power_spectrum_1s = fig.add_subplot(2, 4, 2)
sns.lineplot(x = f_1s[:(N_1s // 2 + 1)], y = log_power_spectrum_1s[:(N_1s // 2 + 1)], ax = ax_power_spectrum_1s, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_1s, title = "はじめの一歩: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_1s, xlim = (0, 25000))

# 3.1.3. 解析信号のエンベロープ
ax_envelope_1s = fig.add_subplot(2, 4, 5)
sns.lineplot(x = t_1s, y = envelope_1s, ax = ax_envelope_1s, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_1s, title = "はじめの一歩: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_1s, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3.1.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_1s = fig.add_subplot(2, 4, 6)
sns.lineplot(x = f_1s[:N_1s // 2 + 1], y = log_fluctuation_spectrum_1s[:N_1s // 2 + 1], ax = ax_fluctuation_spectrum_db_1s, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_1s, title = "はじめの一歩: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_1s, xlim = (0, 25000))


# 3.2. ファイルX
# 3.2.1. 時間波形
ax_X = fig.add_subplot(2, 4, 3)
sns.lineplot(x = t_X, y = yt_X, ax = ax_X, color = "blue", linewidth = 0.1)
plt.setp(ax_X, title = "ファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3.2.2 標準化した対数パワースペクトル
ax_power_spectrum_X = fig.add_subplot(2, 4, 4)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_power_spectrum_X[:(N_X // 2 + 1)], ax = ax_power_spectrum_X, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_X, title = "ファイルX: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_X, xlim = (0, 25000))

# 3.2.3. 解析信号のエンベロープ
ax_envelope_X = fig.add_subplot(2, 4, 7)
sns.lineplot(x = t_X, y = envelope_X, ax = ax_envelope_X, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_X, title = "ファイルX: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3.2.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_X = fig.add_subplot(2, 4, 8)
sns.lineplot(x = f_X[:N_X // 2 + 1], y = log_fluctuation_spectrum_X[:N_X // 2 + 1], ax = ax_fluctuation_spectrum_db_X, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_X, title = "ファイルX: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_X, xlim = (0, 25000))


# レイアウト調整
plt.tight_layout()
plt.show()

スペクトログラム

曲全体の周波数スペクトルでも分析は可能だが、音色の時間変化をとらえることはできない。
音声の一部を切り取ってフーリエ変換を行うことを繰り返すことで、音色の時間変化を捉えることができる。
短時間のフーリエ変換を短時間フーリエ変換と呼び、パワースペクトルの時間変化をスペクトログラムと呼ぶ。

長さ n_fft の短い配列のフーリエ変換を計算した後、 hop_length だけインデックスをずらしてまた短い配列のフーリエ変換を行う。
これらのパラメータはユーザーがいかようにも調整できてしまうので、いろいろと試す必要がある。

スペクトログラムの対数パワースペクトルとプロット用のメッシュグリッドを計算する
# スペクトログラムの対数パワースペクトルとプロット用のメッシュグリッドを計算する
# yt: 時間波形
# sr: サンプリング周波数
# n_fft: 時間波形のフレーム数(1ステップ毎に計算する時間波形のデータ数)
# hop_length: フレーム毎に遷移するデータ数
# window: 窓関数(scipy.signal.windowsの関数オブジェクト)
def spectrogram(yt, sr, n_fft, hop_length, window):
    # スペクトログラム
    D = librosa.stft(yt, n_fft = n_fft, hop_length = hop_length, window = window)
    # 対数パワースペクトル
    DB = librosa.amplitude_to_db(np.abs(D), ref = np.sum, amin = 1e-20, top_db = 80.0)
    # 時間軸
    times = librosa.times_like(D, sr = sr, hop_length = hop_length, n_fft = n_fft)
    # 周波数軸
    frequencies = librosa.fft_frequencies(sr = sr, n_fft = n_fft)
    # メッシュグリッド
    T, F = np.meshgrid(times, frequencies)
    # スペクトログラムの対数パワースペクトルとプロット用のメッシュグリッドを返す
    return DB, T, F
読み込んだ音楽データのスペクトログラムをプロットする
# 1. はじめの一歩について、スペクトログラムを計算する
DB_1s, T_1s, F_1s = spectrogram(yt = yt_1s, sr = sr_1s, n_fft = sr_1s // 2, hop_length = sr_1s // 8, window = hamming)


# 2. ファイルXについて、スペクトログラムを計算する
DB_X, T_X, F_X = spectrogram(yt = yt_X, sr = sr_X, n_fft = sr_X // 2, hop_length = sr_X // 8, window = hamming)


# 3. プロット
# Figure
fig = plt.figure(figsize = (6 * 4, 1 + 5 * 1))# (横, 縦)

# 3.1. はじめの一歩
# 3.1.1. 時間波形
ax_1s = fig.add_subplot(1, 4, 1)
sns.lineplot(x = t_1s, y = yt_1s, ax = ax_1s, color = "blue", linewidth = 0.1)
plt.setp(ax_1s, title = "はじめの一歩", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_1s, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3.1.2. スペクトログラム
ax_3d_1s = fig.add_subplot(1, 4, 2, projection = "3d")
ax_3d_1s.plot_surface(F_1s, T_1s, DB_1s, cmap = "viridis")
plt.setp(ax_3d_1s, title = "はじめの一歩: Spectrogram", xlabel = "Frequency (Hz)", ylabel = "Time (sec)", zlabel = "Normalized Logarithm Power Spectrum (dB)")


# 3.2. ファイルX
# 3.2.1. 時間波形
ax_X = fig.add_subplot(1, 4, 3)
sns.lineplot(x = t_X, y = yt_X, ax = ax_X, color = "blue", linewidth = 0.1)
plt.setp(ax_X, title = "ファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 3.2.2. スペクトログラム
ax_3d_X = fig.add_subplot(1, 4, 4, projection = "3d")
ax_3d_X.plot_surface(F_X, T_X, DB_X, cmap = "viridis")
plt.setp(ax_3d_X, title = "ファイルX: Spectrogram", xlabel = "Frequency (Hz)", ylabel = "Time (sec)", zlabel = "Normalized Logarithm Power Spectrum (dB)")


# レイアウト調整
plt.tight_layout()
plt.show()

窓関数

前節の短時間フーリエ変換の関数に window = hamming というパラメータが出てきた。
窓関数(window)は、短い時間波形の最初と最後の形状の影響がFFTに現れることを抑制するための関数である。

窓関数 w[n] の使い方は簡単で、フーリエ変換を行う時間波形 y[n] に窓関数をかけておくだけである。

\begin{align*} Y_{w}[k] &= \hat{F}[w[n]y[n]] \end{align*}

よく使われるのがハミング窓である。

w_{\text{Hamming}}[n] = 0.54 - 0.46 \cos (\frac{2 \pi n}{N - 1})

コード

データ数が多い場合窓関数は無くても問題ない。
端の影響を無視できなくなるデータ数は調査が不足している。

窓関数をかけたパワースペクトルや解析信号をプロットする
# 1. 窓関数をかけたパワースペクトルや解析信号を計算する
# 1.1. 窓関数
w_X = hamming(N_X)# SciPy
w_X = np.hamming(N_X)# NumPy
yt_w_X = w_X * yt_X
# 1.2. 周波数スペクトル
yf_w_X, f_w_X = fft(yt = yt_w_X, sr = sr_X)
# 1.3. パワースペクトル、対数パワースペクトル
amplitude_spectrum_w_X, power_spectrum_w_X, log_power_spectrum_w_X = spectrum(yf = yf_w_X)
# 1.4. 解析信号
envelope_w_X, fluctuation_spectrum_w_X, log_fluctuation_spectrum_w_X = analytic_signal(yt = yt_w_X)


# 2. プロット
# Figure
fig = plt.figure(figsize = (6 * 4, 1 + 5 * 2))# (横, 縦)

# 2.1. ファイルX
# 2.1.1. 時間波形
ax_X = fig.add_subplot(2, 4, 1)
sns.lineplot(x = t_X, y = yt_X, ax = ax_X, color = "blue", linewidth = 0.1)
plt.setp(ax_X, title = "ファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.2 標準化した対数パワースペクトル
ax_power_spectrum_X = fig.add_subplot(2, 4, 2)
sns.lineplot(x = f_X[:(N_X // 2)], y = log_power_spectrum_X[:(N_X // 2)], ax = ax_power_spectrum_X, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_X, title = "ファイルX: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_X, xlim = (0, 25000))

# 2.1.3. 解析信号のエンベロープ
ax_envelope_X = fig.add_subplot(2, 4, 5)
sns.lineplot(x = t_X, y = envelope_X, ax = ax_envelope_X, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_X, title = "ファイルX: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.1.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_X = fig.add_subplot(2, 4, 6)
sns.lineplot(x = f_X[:(N_X // 2)], y = log_fluctuation_spectrum_X[:(N_X // 2)], ax = ax_fluctuation_spectrum_db_X, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_X, title = "ファイルX: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_X, xlim = (0, 25000))


# 2.2. 窓関数をかけたファイルX
# 2.2.1. 時間波形
ax_w_X = fig.add_subplot(2, 4, 3)
sns.lineplot(x = t_X, y = yt_w_X, ax = ax_w_X, color = "blue", linewidth = 0.1)
plt.setp(ax_w_X, title = "窓関数をかけたファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_w_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.2 標準化した対数パワースペクトル
ax_power_spectrum_w_X = fig.add_subplot(2, 4, 4)
sns.lineplot(x = f_X[:(N_X // 2)], y = log_power_spectrum_w_X[:(N_X // 2)], ax = ax_power_spectrum_w_X, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_w_X, title = "窓関数をかけたファイルX: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_w_X, xlim = (0, 25000))

# 2.2.3. 解析信号のエンベロープ
ax_envelope_w_X = fig.add_subplot(2, 4, 7)
sns.lineplot(x = t_X, y = envelope_w_X, ax = ax_envelope_w_X, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_w_X, title = "窓関数をかけたファイルX: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_w_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_w_X = fig.add_subplot(2, 4, 8)
sns.lineplot(x = f_X[:(N_X // 2)], y = log_fluctuation_spectrum_w_X[:(N_X // 2)], ax = ax_fluctuation_spectrum_db_w_X, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_w_X, title = "窓関数をかけたファイルX: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_w_X, xlim = (0, 25000))


# レイアウト調整
plt.tight_layout()
plt.show()

高域協調(プリエンファシス)

一般に高周波帯域の波形ほど音圧レベルが減衰しやすいらしい。
実際に前節までのパワースペクトルも高周波帯域の音圧レベルが低かった。
高周波帯域の成分を強調する方法として高域協調(プリエンファシス)という手段がある。
定義式は畳み込みだが、よく使われるのは差分の式である。

\begin{align*} \text{pre\_emphasis} (x[n]) &= x[n] - \alpha x[n -1]\\ &0 < \alpha < 1 \end{align*}

コード

高周波成分がカギを握る場合は、この処理は必須であると感じた。

高域協調処理を施したパワースペクトルや解析信号をプロットする
# 1. 高域協調処理を施したパワースペクトルや解析信号を計算する
# 1.1. 高域協調処理
yt_pe_X = librosa.effects.preemphasis(yt_X, coef = 0.97)
# 1.2. 高域協調処理を施した周波数スペクトル
yf_pe_X, f_pe_X = fft(yt =yt_pe_X, sr = sr_X)
# 1.3. パワースペクトル、対数パワースペクトル
amplitude_spectrum_pe_X, power_spectrum_pe_X, log_power_spectrum_pe_X = spectrum(yf = yf_w_X)
# 1.4. 解析信号
envelope_pe_X, fluctuation_spectrum_pe_X, log_fluctuation_spectrum_pe_X = analytic_signal(yt = yt_pe_X)


# 2. プロット
# Figure
fig = plt.figure(figsize = (6 * 4, 1 + 5 * 2))# (横, 縦)

# 2.1. ファイルX
# 2.1.1. 時間波形
ax_X = fig.add_subplot(2, 4, 1)
sns.lineplot(x = t_X, y = yt_X, ax = ax_X, color = "blue", linewidth = 0.1)
plt.setp(ax_X, title = "ファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.2 標準化した対数パワースペクトル
ax_power_spectrum_X = fig.add_subplot(2, 4, 2)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_power_spectrum_X[:(N_X // 2 + 1)], ax = ax_power_spectrum_X, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_X, title = "ファイルX: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_X, xlim = (0, 25000))

# 2.1.3. 解析信号のエンベロープ
ax_envelope_X = fig.add_subplot(2, 4, 5)
sns.lineplot(x = t_X, y = envelope_X, ax = ax_envelope_X, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_X, title = "ファイルX: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.1.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_X = fig.add_subplot(2, 4, 6)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_fluctuation_spectrum_X[:(N_X // 2 + 1)], ax = ax_fluctuation_spectrum_db_X, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_X, title = "ファイルX: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_X, xlim = (0, 25000))


# 2.2. 高域協調処理を施したファイルX
# 2.2.1. 時間波形
ax_pe_X = fig.add_subplot(2, 4, 3)
sns.lineplot(x = t_X, y = yt_pe_X, ax = ax_pe_X, color = "blue", linewidth = 0.1)
plt.setp(ax_pe_X, title = "高域協調処理を施したファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_pe_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.2 標準化した対数パワースペクトル
ax_power_spectrum_pe_X = fig.add_subplot(2, 4, 4)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_power_spectrum_pe_X[:(N_X // 2 + 1)], ax = ax_power_spectrum_pe_X, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_pe_X, title = "高域協調処理を施したファイルX: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_pe_X, xlim = (0, 25000))

# 2.2.3. 解析信号のエンベロープ
ax_envelope_pe_X = fig.add_subplot(2, 4, 7)
sns.lineplot(x = t_X, y = envelope_pe_X, ax = ax_envelope_pe_X, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_pe_X, title = "高域協調処理を施したファイルX: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_pe_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_pe_X = fig.add_subplot(2, 4, 8)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_fluctuation_spectrum_pe_X[:(N_X // 2 + 1)], ax = ax_fluctuation_spectrum_db_pe_X, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_pe_X, title = "高域協調処理を施したファイルX: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_pe_X, xlim = (0, 25000))


# レイアウト調整
plt.tight_layout()
plt.show()

フィルタバンク分析

この節のモチベーションは、機械学習のオーバーフィッティングを防ぐことである。
固有振動数の確認などグラフを用いた分析が中心の場合、この節の処理は不要である。

データ数 N の時間波形の周波数スペクトルのデータ数は N/2 + 1 である。
例えば、サンプリング周波数48(kHz)で10秒間録音した時間波形の周波数スペクトルのデータ数は N/2 + 1 = 240,001 である。
24万次元のデータを機械学習モデルに学習させるとオーバーフィッティングして実用的な予測ができなくなってしまう。

機械学習モデルで上手く学習できる程度まで周波数スペクトルを集約させる方法を、フィルタバンク分析と呼ぶ。
周波数スペクトルを集約する方法として、人間の聴覚すなわち音の聞こえ方に基づいた尺度で集約する。
この尺度の1つがメル尺度であり、メル尺度の差が同じであれば人間が感じる音高の差が同一であるという尺度である。

コード

メルフィルタバンクと呼び周波数を作成する関数
# メルフィルタバンクと呼び周波数を作成する関数
# sr: サンプリング周波数
# N: 時間波形のデータ数
# f: 周波数軸
# n_mels: メルフィルタバンクの数
def make_mel_filter_bank(sr, N, f, n_mels):
    # メルフィルタバンク
    mel_filter_bank = librosa.filters.mel(sr = sr, n_fft = N, n_mels = n_mels)
    # 最大値が1になるように正規化する
    mel_filter_bank /= np.max(mel_filter_bank, axis = -1)[:, None]

    # 呼び周波数
    mel_filter_bank_label = list()
    for idx in range(mel_filter_bank.shape[0]):
        idx_argmax = np.argmax(mel_filter_bank[idx])
        mel_filter_bank_label.append(np.round(f[idx_argmax], decimals = 0))
    mel_filter_bank_label = np.array(mel_filter_bank_label)

    return mel_filter_bank, mel_filter_bank_label
メルフィルタバンク特徴量と対数メルフィルタバンク特徴量を計算する
# メルフィルタバンクと呼び周波数を作成する関数
# sr: サンプリング周波数
# f: 周波数軸
# n_mels: メルフィルタバンクの数
def make_mel_filter_bank(sr, f, n_mels):
    # データ数
    N = len(f)
    # メルフィルタバンク
    mel_filter_bank = librosa.filters.mel(sr = sr, n_fft = N, n_mels = n_mels)
    # 最大値が1になるように正規化する
    mel_filter_bank /= np.max(mel_filter_bank, axis = -1)[:, None]
    # 呼び周波数の配列を作成する
    mel_filter_bank_label = list()
    for idx in range(mel_filter_bank.shape[0]):
        idx_argmax = np.argmax(mel_filter_bank[idx])
        mel_filter_bank_label.append(np.round(f[idx_argmax], decimals = 0))
    mel_filter_bank_label = np.array(mel_filter_bank_label)
    # メルフィルタバンクと呼び周波数を返す
    return mel_filter_bank, mel_filter_bank_label
メルフィルタバンクを通した特徴量をプロットする
# 1. メルフィルタバンクを通した特徴量を計算する
# 1.1. メルフィルタ
mel_filter_bank_X, mel_filter_bank_label_X = make_mel_filter_bank(sr = sr_X, f = f_X, n_mels = 16)
# 1.2. メルフィルタを通したパワースペクトル、対数パワースペクトル
mel_power_spectrum_X, fbank_X = make_fbank(yf = yf_X, mel_filter_bank = mel_filter_bank_X)
# 1.3. メルフィルタを通した解析信号のパワースペクトル、対数パワースペクトル
mel_power_spectrum_envelope_X, fbank_envelope_X = make_fbank(yf = np.fft.fft(envelope_X), mel_filter_bank = mel_filter_bank_X)


# 2. プロット
# Figure
fig = plt.figure(figsize = (6 * 4, 1 + 5 * 2))# (横, 縦)

# 2.1. ファイルX
# 2.1.1. 時間波形
ax_X = fig.add_subplot(2, 4, 1)
sns.lineplot(x = t_X, y = yt_X, ax = ax_X, color = "blue", linewidth = 0.1)
plt.setp(ax_X, title = "ファイルX", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.2.2 標準化した対数パワースペクトル
ax_power_spectrum_X = fig.add_subplot(2, 4, 2)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_power_spectrum_X[:(N_X // 2 + 1)], ax = ax_power_spectrum_X, color = "blue", linewidth = 0.1)
plt.setp(ax_power_spectrum_X, title = "ファイルX: Logarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_X, xlim = (0, 25000), ylim = (-160, 0))

# 2.1.3. 解析信号のエンベロープ
ax_envelope_X = fig.add_subplot(2, 4, 5)
sns.lineplot(x = t_X, y = envelope_X, ax = ax_envelope_X, color = "blue", linewidth = 0.1)
plt.setp(ax_envelope_X, title = "ファイルX: エンベロープ", xlabel = "Time (sec)", ylabel = "Amplitude")
plt.setp(ax_envelope_X, xlim = (0, 10), xticks = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ylim = (-1, 1), yticks = [-1, -0.5, 0, 0.5, 1])

# 2.1.4. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_X = fig.add_subplot(2, 4, 6)
sns.lineplot(x = f_X[:(N_X // 2 + 1)], y = log_fluctuation_spectrum_X[:(N_X // 2 + 1)], ax = ax_fluctuation_spectrum_db_X, color = "blue", linewidth = 0.1)
plt.setp(ax_fluctuation_spectrum_db_X, title = "ファイルX: エンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_X, xlim = (0, 25000), ylim = (-160, 0))


# メルフィルタを通したファイルX
# 2.2.1 標準化した対数パワースペクトル
ax_power_spectrum_mel_X = fig.add_subplot(2, 4, 3)
sns.scatterplot(x = mel_filter_bank_label_X, y = fbank_X, ax = ax_power_spectrum_mel_X, color = "black", s = 100)
plt.setp(ax_power_spectrum_mel_X, title = "ファイルX: メルフィルタを通したLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_power_spectrum_mel_X, xlim = (0, 25000), ylim = (-160, 0))

# 2.2.2. エンベロープの対数パワースペクトル
ax_fluctuation_spectrum_db_X = fig.add_subplot(2, 4, 7)
sns.scatterplot(x = mel_filter_bank_label_X, y = fbank_envelope_X, ax = ax_fluctuation_spectrum_db_X, color = "black", s = 100)
plt.setp(ax_fluctuation_spectrum_db_X, title = "ファイルX: メルフィルタを通したエンベロープのLogarithm Power Spectrum", xlabel = "Frequency (Hz)", ylabel = "Normalized Logarithm Power Spectrum (dB)")
plt.setp(ax_fluctuation_spectrum_db_X, xlim = (0, 25000), ylim = (-160, 0))


# レイアウト調整
plt.tight_layout()
plt.show()

終わりに

基本的な振動解析手法をコードつきでまとめました。
ご参考になれば幸いです。

参考書籍・URL

Pythonライブラリ

使用したライブラリの関数URL

Discussion