振動波形のヒルベルト変換から対数減衰率を求めるPythonコード

減衰自由振動現象の重要な特性に減衰率があります。減衰特性は物性から求めることが難しいので、実際に系を振動させて測定することが一般的です。ここでは減衰自由振動の波形に対してヒルベルト変換を行い、対数減衰率を計算するPythonコードを紹介します。

こんにちは。wat(@watlablog)です。ここではPythonを使って振動波形から減衰率を測定してみます！

目次（項目クリックでジャンプできます）

減衰自由振動のおさらい

減衰自由振動とは？

減衰自由振動とは、振動を繰り返すにつれて徐々に振幅が小さくなっていく現象を指します。下図に減衰自由振動の例を示します。

図に示すように、 $ζ$ の数値によって振動波形の様相が変わります。 $ζ$ は減衰率（減衰比）と呼ばれるパラメータで、大きいほど振動がはやく止まる系となります。

対数減衰率とは？

減衰自由振動時の各振動毎における振幅は、隣り合う山（周期 $T$ 後の振幅値）が指数関数曲線に従って減少する特徴を持ちます。

そのため、振幅値の対数（デシベル変換も同義）をとることで、その減少曲線を直線にできます。

対数減衰率 $δ$ とは、隣り合う振幅同士の比に対し自然対数をとった式（1）で示されます $^{[1]}$ 。

\begin{matrix} (1) & δ = \ln \frac{a_{n}}{a_{n + 1}} \end{matrix}

また1周期 $T$ の振幅比ではなく、 $m$ 周期 $m T$ との振幅比を使って精度を高める式（2）の方法もあります $^{[2]}$ 。

\begin{matrix} (2) & δ = \frac{1}{m} \ln \frac{a_{n}}{a_{n + m}} \end{matrix}

対数減衰率と減衰比の関係

振動系は1自由度や多自由度の質点、ばね、ダンパーでモデル化することがあります。ダンパーの減衰係数と系の振動現象そのものに密接に関係する減衰比 $ζ$ と対数減衰率 $δ$ には式（3）の関係があります $^{[1]}$ 。

\begin{matrix} (3) & ζ = \frac{δ}{2 π} \end{matrix}

つまり対数減衰率を計測から求めることができれば、減衰比を算出でき、振動モデルのパラメータ同定が可能となります。

本記事では振動波形を計測した結果を使って自分で対数減衰率を求められるようになることを目標とします。

また、上式（1）や（2）はピーク値を手動でピックアップするかピーク検出プログラムを書く必要がありますが、今回はより連続的かつ自動的に振幅を抽出可能なヒルベルト変換を使います。

詳しくはコードを書きながら説明しますので、是非記事後半もご覧ください。

振動波形から対数減衰率と減衰比を求めるPythonコード

動作環境

この記事のコードは以下のPython環境で動作確認を行いました。

Python	Python 3.9.6
PyCharm (IDE)	PyCharm CE 2020.1
Numpy	1.21.1
Scipy	1.4.1
matplotlib	3.4.3

サンプルの減衰自由振動波形を作る

コピペ動作可能な全コードは最後に示しますが、まずはどのような手順でコードが書かれているのかを一つずつ説明します。

まずは自分で書いたコードを検証できるように、サンプルの減衰自由振動波形を作ってみましょう。

理論式： $0 < ζ < 1$

コードを書く前に、理論式を確認します。既に当ブログでも紹介しましたが、減衰比の値により振動の様相が変化します。

$0 < ζ < 1$ の場合、固有角振動数 $ω$ 、減衰固有角振動数 $ω_{d}$ を定義すると式（4）-式（9）を使って応答変位 $x$ を求めることができます $^{[3]}$ 。

\begin{matrix} (4) & ω = \sqrt{\frac{k}{m}} \end{matrix}

\begin{matrix} (5) & ω_{d} = ω \sqrt{1 - ζ^{2}} \end{matrix}

\begin{matrix} (6) & σ = ω ζ \end{matrix}

\begin{matrix} (7) & X = \sqrt{x_{0}^{2} + {(\frac{v_{0} + σ x_{0}}{ω_{d}})}^{2}} \end{matrix}

\begin{matrix} (8) & ϕ = \tan^{- 1} (\frac{v_{0} + σ x_{0}}{x_{0} ω_{d}}) \end{matrix}

\begin{matrix} (9) & x = e^{- σ t} X \cos (ω_{d} t - ϕ) \end{matrix}

理論式： $ζ = 1$

同様に、式（10）の無次元化速度 $σ_{v 0}$ を導入することで $ζ = 1$ の場合の応答変位は式（11）で計算可能です $^{[4]}$ 。

\begin{matrix} (10) & σ_{v 0} = \frac{v_{0}}{x_{0} ω} \end{matrix}

\begin{matrix} (11) & x = x_{0} e^{- ω t} {(σ_{v 0} + 1) ω t + 1} \end{matrix}

理論式： $ζ > 1$

さらに $ζ > 1$ の場合は式（12）が応答変位の式になります $^{[4]}$ （これはかなりややこしや）。

\begin{matrix} (12) & x = x_{0} e^{- ζ ω t} {\cosh (ω t \sqrt{ζ^{2} - 1}) + \frac{σ_{v 0} + ζ}{\sqrt{ζ^{2} - 1}} \sinh (ω t \sqrt{ζ^{2} - 1})} \end{matrix}

理論波形作成を関数にまとめたPythonコード

今回は式（9）までの内容しか検証で使いませんが、一連の式を関数damp_eq()にまとめてみました。質量 $m$ 、剛性 $k$ や初期変位 $x_{0}$ 、初期速度 $v_{0}$ を受け取って波形を返す仕様です。減衰係数 $c$ ではなく減衰比 $ζ$ を与える仕様にしています。

import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt

def damp_eq(m, k, zeta, x0, v0, t):
    ''' 減衰自由振動の理論値を計算 '''

    omega = np.sqrt(k / m)
    sigma = omega * zeta
    sigma_v0 = v0 / (x0 * omega)

    # 減衰比の場合分け
    if zeta < 1:
        # 減衰振動
        omega_d = omega * np.sqrt(1 - zeta ** 2)
        X = np.sqrt(x0 ** 2 + ((v0 + sigma * x0) / omega_d))
        phi = np.arctan2((v0 + sigma * x0), (x0 * omega_d))
        x = np.exp(-sigma * t) * X * np.cos(omega_d * t - phi)
    elif zeta == 1:
        # 臨界減衰
        omega_d = omega
        x = x0 * np.exp(-omega * t) * ((sigma_v0 + 1) * omega * t + 1)
    else:
        # 過減衰
        omega_d = omega
        x = x0 * np.exp(- zeta * omega * t) * (np.cosh(omega * t * np.sqrt(zeta ** 2 - 1))
                                               + (sigma_v0 + zeta) / (np.sqrt(zeta ** 2 + 1))
                                               * np.sinh(omega * t * np.sqrt(zeta ** 2 - 1)))

    return x, omega_d

if __name__ == '__main__':
    ''' メイン処理 '''

    # サンプルの減衰振動を計算(m:質量[kg], k:ばね定数[N/m], zeta:減衰率, x0:初期変位[m], v0:初速度[m/s], dt:時間刻み[s])
    m = 0.5
    k = 1000
    zeta = 0.02
    x0 = 1
    v0 = 0
    dt = 0.001
    t = np.arange(0, 10, dt)
    x, omega_d = damp_eq(m, k, zeta, x0, v0, t)

import numpy as np

from scipy import fftpack

from matplotlib import pyplot as plt

def damp_eq(m, k, zeta, x0, v0, t):

''' 減衰自由振動の理論値を計算 '''

omega = np.sqrt(k / m)

sigma = omega * zeta

sigma_v0 = v0 / (x0 * omega)

# 減衰比の場合分け

if zeta < 1:

# 減衰振動

omega_d = omega * np.sqrt(1 - zeta ** 2)

X = np.sqrt(x0 ** 2 + ((v0 + sigma * x0) / omega_d))

phi = np.arctan2((v0 + sigma * x0), (x0 * omega_d))

x = np.exp(-sigma * t) * X * np.cos(omega_d * t - phi)

elif zeta == 1:

# 臨界減衰

omega_d = omega

x = x0 * np.exp(-omega * t) * ((sigma_v0 + 1) * omega * t + 1)

else:

# 過減衰

omega_d = omega

x = x0 * np.exp(- zeta * omega * t) * (np.cosh(omega * t * np.sqrt(zeta ** 2 - 1))

+ (sigma_v0 + zeta) / (np.sqrt(zeta ** 2 + 1))

* np.sinh(omega * t * np.sqrt(zeta ** 2 - 1)))

return x, omega_d

if __name__ == '__main__':

''' メイン処理 '''

# サンプルの減衰振動を計算(m:質量[kg], k:ばね定数[N/m], zeta:減衰率, x0:初期変位[m], v0:初速度[m/s], dt:時間刻み[s])

m = 0.5

k = 1000

zeta = 0.02

x0 = 1

v0 = 0

dt = 0.001

t = np.arange(0, 10, dt)

x, omega_d = damp_eq(m, k, zeta, x0, v0, t)

プロットを含めたコードは最後に示しますが、実行すると以下の波形を得ます。THE・減衰振動、って感じの波形にしてみました。

ヒルベルト変換で瞬時振幅を求める

次にヒルベルト変換を使って包絡線（Envelope）を求めます。包絡線を求めることで、ピーク検出アルゴリズムを使わないで自動的に減衰自由振動の振幅変化を取得できます。

ヒルベルト変換の考え方

ヒルベルト変換は波形を実部と虚部に分離して考えます。このように考えることで時間が進行していく最中の瞬時振幅、瞬時位相、瞬時周波数をそれぞれ求めることができ、時間領域の信号処理で重宝されています。

難しそうに聞こえますが、やっていることはフーリエ変換ですので、当ブログの読者であればすぐにできます。詳しく知りたい方は「Pythonでヒルベルト変換!時間波形の包絡線を求める方法」をご覧ください。

ヒルベルト変換を行うPythonコード

以下にヒルベルト変換を行うPythonコードを記載します。SciPyで1行で書けたりしますが、せっかくフルスクラッチしたので今回はこちらの方を紹介します。詳しくは参考文献 $^{[5]}$ をどうぞ。

def hilbert(yt, dt):
    ''' ヒルベルト変換 '''

    # フーリエ変換
    spectrum = fftpack.fft(yt)

    # 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）
    zeros = np.zeros(int(len(spectrum) / 2))
    spectrum[int(len(spectrum) / 2):len(spectrum)] = zeros

    # 正周波数域（ナイキスト周波数まで）を2倍する
    spectrum *= 2

    # 実部と虚部から瞬時振幅と瞬時位相を計算する
    zt = fftpack.ifft(spectrum)
    re = zt.real
    im = zt.imag

    # 瞬時振幅
    instantaneous_amp = np.sqrt(re ** 2 + im ** 2)

    # 瞬時位相とアンラップ処理
    instantaneous_phase = np.arctan2(im, re)
    instantaneous_phase = np.unwrap(instantaneous_phase)

    # 瞬時周波数（位相角の微分と時間刻み影響を考慮）
    instantaneous_frequency = ((1 / dt) / (2 * np.pi)) * np.gradient(instantaneous_phase)

    return instantaneous_amp

def hilbert(yt, dt):

''' ヒルベルト変換 '''

# フーリエ変換

spectrum = fftpack.fft(yt)

# 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）

zeros = np.zeros(int(len(spectrum) / 2))

spectrum[int(len(spectrum) / 2):len(spectrum)] = zeros

# 正周波数域（ナイキスト周波数まで）を2倍する

spectrum *= 2

# 実部と虚部から瞬時振幅と瞬時位相を計算する

zt = fftpack.ifft(spectrum)

re = zt.real

im = zt.imag

# 瞬時振幅

instantaneous_amp = np.sqrt(re ** 2 + im ** 2)

# 瞬時位相とアンラップ処理

instantaneous_phase = np.arctan2(im, re)

instantaneous_phase = np.unwrap(instantaneous_phase)

# 瞬時周波数（位相角の微分と時間刻み影響を考慮）

instantaneous_frequency = ((1 / dt) / (2 * np.pi)) * np.gradient(instantaneous_phase)

return instantaneous_amp

このコードを実行することで包絡線（瞬時振幅）を得ます。これを元の波形と重ね描きすることでしっかりと包絡線になっていることがわかります。ただ、波形の両端部が跳ねているという点も確認できます。

解析範囲を自動調整する

そのままカーブフィットすると精度が落ちる

瞬時振幅が求められたので、あとは直線近似（カーブフィット）するだけと思いきや、素直にヒルベルト変換した結果を使って対数減衰率を算出しようとすると精度が悪くなっていることに気づきます。

これは図解すると一目瞭然であり、瞬時振幅の計算時にできる後端部の跳ね上がりや、ほぼ振動しなくなった振幅極小部分で対数振幅が振動することに起因します。

これは意外とやっかいで、系の固有振動数や最大振幅値、データ長にも依存します。

そのため元波形を使って事前に「安定して減衰振動している範囲」を自動抽出する関数を作って実行します。

今回はstable_check()という関数を作り、まずデータを半分捨てます（ここはちょっと改良の余地あり？）。次に最大振幅の10[%]以上の振幅を持つ成分のみを残すことで、対数変換後の精度悪化要因を排除します。

データの全範囲で最大振幅の10[％]に到達しない場合はデータ長の半分までが安定した振動をしていると判断されます。

解析範囲を自動調整するPythonコード

以下の関数がデータ範囲自動調整コードです。np.whereを使うことでindex探査しています。

def stable_check(t, y, i_amp, rate):
    ''' 安定して減衰振動している範囲を抽出する '''

    # データの後半を切り捨て（包絡線の端部は瞬時振幅が高くなってしまうため）
    t_half = t[0:int(len(t) / 2)]
    y_half = y[0:int(len(y) / 2)]

    # 最大値に対するrate割合を超える振幅成分のみ抽出する
    peak = np.max(np.abs(y_half))
    threshold = peak * rate
    index_max = np.max(np.where(np.abs(y_half) > threshold))

    # 安定して減衰振動している範囲を抽出
    t_censored = t_half[0:index_max]
    y_censored = y_half[0:index_max]
    i_amp = i_amp[0:index_max]

    return t_censored, i_amp

def stable_check(t, y, i_amp, rate):

''' 安定して減衰振動している範囲を抽出する '''

# データの後半を切り捨て（包絡線の端部は瞬時振幅が高くなってしまうため）

t_half = t[0:int(len(t) / 2)]

y_half = y[0:int(len(y) / 2)]

# 最大値に対するrate割合を超える振幅成分のみ抽出する

peak = np.max(np.abs(y_half))

threshold = peak * rate

index_max = np.max(np.where(np.abs(y_half) > threshold))

# 安定して減衰振動している範囲を抽出

t_censored = t_half[0:index_max]

y_censored = y_half[0:index_max]

i_amp = i_amp[0:index_max]

return t_censored, i_amp

対数減衰率と減衰比を計算する

振幅をデシベルに変換する

あとは式に沿って計算を進めるだけです。まずは応答変位（ここは加速度や速度でもOK）振幅 $x$ を式（13）でデシベル変換します。デシベル基準値 $dBref$ は任意値でも問題ありません。

\begin{matrix} (13) & x_{dB} = 20 \log_{10} \frac{x}{dBref} \end{matrix}

カーブフィットする

自動範囲調節された瞬時振幅に対してカーブフィットすると、下図になります。

計算にはこの傾き情報を使います。カーブフィットの方法は「Pythonでカーブフィット!最小二乗法で直線近似する方法」 $^{[6]}$ をご覧ください。

前端部も振幅が跳ねていますので、この部分を考慮して解析範囲を調整できればもう少し精度が上がりそうです。しかし実際はそれほど結果に影響していないようです。

対数減衰率と減衰比を計算する

減衰固有振動数 $f_{d}$ を式（14）、カーブフィットした結果の傾きを $a$ とすると、式（15）で対数減衰率 $δ$ を求めることが可能です $^{[2]}$ 。ここまで求まれば式（3）で減衰比 $ζ$ も計算できます。

\begin{matrix} (14) & f_{d} = \frac{1}{2 π} ω_{d} \end{matrix}

\begin{matrix} (15) & δ = \frac{a}{f_{d}} \frac{\ln 10}{20} \end{matrix}

対数減衰率と減衰比を求めるPythonコード

calc_damping_factor()関数に対数減衰率と減衰比を計算する一連のコードをまとめました。

def calc_damping_factor(t, y, omega_d):
    ''' 対数減衰率と減衰比を計算（ヒルベルト変換結果を用いて）'''

    # デシベル変換（dBrefは何でも良い）
    y_dB = 20 * np.log10(y / 1e-6)

    # カーブフィット
    coefficients = np.polyfit(t, y_dB, 1)
    damp_curve = coefficients[0] * t + coefficients[1]
    plot([t, t], [y_dB, damp_curve], ['envelope[dB]', 'curvefit'], 't[s]', 'Displacement[dB]')

    # 対数減衰率deltaと減衰比zetaを計算
    a = -coefficients[0]
    fd = omega_d * (1 / (2 * np.pi))
    delta = (a / fd) * (np.log(10) / 20)
    zeta = delta / (2 * np.pi)

    return zeta

def calc_damping_factor(t, y, omega_d):

''' 対数減衰率と減衰比を計算（ヒルベルト変換結果を用いて）'''

# デシベル変換（dBrefは何でも良い）

y_dB = 20 * np.log10(y / 1e-6)

# カーブフィット

coefficients = np.polyfit(t, y_dB, 1)

damp_curve = coefficients[0] * t + coefficients[1]

plot([t, t], [y_dB, damp_curve], ['envelope[dB]', 'curvefit'], 't[s]', 'Displacement[dB]')

# 対数減衰率deltaと減衰比zetaを計算

a = -coefficients[0]

fd = omega_d * (1 / (2 * np.pi))

delta = (a / fd) * (np.log(10) / 20)

zeta = delta / (2 * np.pi)

return zeta

全コード（コピペ可：プロット含む）

以下に全コードを示します。このコードにはmatplotlibによるプロットを含みます。実行すると順次グラフが出てきます。一つずつ×印でウィンドウを消すと次の処理に進みます。

import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt

def damp_eq(m, k, zeta, x0, v0, t):
    ''' 減衰自由振動の理論値を計算 '''

    omega = np.sqrt(k / m)
    sigma = omega * zeta
    sigma_v0 = v0 / (x0 * omega)

    # 減衰比の場合分け
    if zeta < 1:
        # 減衰振動
        omega_d = omega * np.sqrt(1 - zeta ** 2)
        X = np.sqrt(x0 ** 2 + ((v0 + sigma * x0) / omega_d))
        phi = np.arctan2((v0 + sigma * x0), (x0 * omega_d))
        x = np.exp(-sigma * t) * X * np.cos(omega_d * t - phi)
    elif zeta == 1:
        # 臨界減衰
        omega_d = omega
        x = x0 * np.exp(-omega * t) * ((sigma_v0 + 1) * omega * t + 1)
    else:
        # 過減衰
        omega_d = omega
        x = x0 * np.exp(- zeta * omega * t) * (np.cosh(omega * t * np.sqrt(zeta ** 2 - 1))
                                               + (sigma_v0 + zeta) / (np.sqrt(zeta ** 2 + 1))
                                               * np.sinh(omega * t * np.sqrt(zeta ** 2 - 1)))

    return x, omega_d

def hilbert(yt, dt):
    ''' ヒルベルト変換 '''

    # フーリエ変換
    spectrum = fftpack.fft(yt)

    # 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）
    zeros = np.zeros(int(len(spectrum) / 2))
    spectrum[int(len(spectrum) / 2):len(spectrum)] = zeros

    # 正周波数域（ナイキスト周波数まで）を2倍する
    spectrum *= 2

    # 実部と虚部から瞬時振幅と瞬時位相を計算する
    zt = fftpack.ifft(spectrum)
    re = zt.real
    im = zt.imag

    # 瞬時振幅
    instantaneous_amp = np.sqrt(re ** 2 + im ** 2)

    # 瞬時位相とアンラップ処理
    instantaneous_phase = np.arctan2(im, re)
    instantaneous_phase = np.unwrap(instantaneous_phase)

    # 瞬時周波数（位相角の微分と時間刻み影響を考慮）
    instantaneous_frequency = ((1 / dt) / (2 * np.pi)) * np.gradient(instantaneous_phase)

    return instantaneous_amp

def stable_check(t, y, i_amp, rate):
    ''' 安定して減衰振動している範囲を抽出する '''

    # データの後半を切り捨て（包絡線の端部は瞬時振幅が高くなってしまうため）
    t_half = t[0:int(len(t) / 2)]
    y_half = y[0:int(len(y) / 2)]

    # 最大値に対するrate割合を超える振幅成分のみ抽出する
    peak = np.max(np.abs(y_half))
    threshold = peak * rate
    index_max = np.max(np.where(np.abs(y_half) > threshold))

    # 安定して減衰振動している範囲を抽出
    t_censored = t_half[0:index_max]
    y_censored = y_half[0:index_max]
    i_amp = i_amp[0:index_max]
    plot([t, t_half, t_censored], [y, y_half, y_censored], ['original', 'half data', 'censored data'], 't[s]', 'Displacement[m]')

    return t_censored, i_amp

def calc_damping_factor(t, y, omega_d):
    ''' 対数減衰率と減衰比を計算（ヒルベルト変換結果を用いて）'''

    # デシベル変換（dBrefは何でも良い）
    y_dB = 20 * np.log10(y / 1e-6)

    # カーブフィット
    coefficients = np.polyfit(t, y_dB, 1)
    damp_curve = coefficients[0] * t + coefficients[1]
    plot([t, t], [y_dB, damp_curve], ['envelope[dB]', 'curvefit'], 't[s]', 'Displacement[dB]')

    # 対数減衰率deltaと減衰比zetaを計算
    a = -coefficients[0]
    fd = omega_d * (1 / (2 * np.pi))
    delta = (a / fd) * (np.log(10) / 20)
    zeta = delta / (2 * np.pi)

    return zeta

def plot(t, x, label, xlabel, ylabel):
    ''' シミュレーション結果をプロットする '''

    # フォントの種類とサイズを設定する。
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    # 目盛を内側にする。
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    # グラフの上下左右に目盛線を付ける。
    fig = plt.figure()

    ax1 = fig.add_subplot(111)
    ax1.yaxis.set_ticks_position('both')
    ax1.xaxis.set_ticks_position('both')

    # 軸のラベルを設定する。
    ax1.set_xlabel(xlabel)
    ax1.set_ylabel(ylabel)

    # スケールを設定する。
    #ax1.set_xlim(0, 4)

    # プロットを行う。
    for i in range(len(x)):
        ax1.plot(t[i], x[i], label=label[i], lw=1)
    ax1.legend()

    # レイアウト設定
    fig.tight_layout()

    # グラフを表示する。
    plt.show()
    plt.close()

    return

if __name__ == '__main__':
    ''' メイン処理 '''

    # サンプルの減衰振動を計算(m:質量[kg], k:ばね定数[N/m], zeta:減衰率, x0:初期変位[m], v0:初速度[m/s], dt:時間刻み[s])
    m = 0.5
    k = 1000
    zeta = 0.02
    x0 = 1
    v0 = 0
    dt = 0.001
    t = np.arange(0, 10, dt)
    x, omega_d = damp_eq(m, k, zeta, x0, v0, t)
    plot([t], [x], ['Damping vibration'], 't[s]', 'Displacement[m]')

    # ヒルベルト変換から包絡線（瞬時振幅）を計算
    i_amp = hilbert(x, dt)
    plot([t, t], [x, i_amp], ['disp', 'envelope'], 't[s]', 'Displacement[m]')

    # 減衰係数を計算する範囲を自動抽出(rate:データを抽出する最大振幅に対する割合)
    t, i_amp = stable_check(t, x, i_amp, 0.1)

    # 減衰係数を計算
    zeta_measure = calc_damping_factor(t, i_amp, omega_d)

    # 理論値との精度を確認
    print('Theoretical zeta =', zeta)
    print('Measured zeta =', zeta_measure)
    print('Accuracy =', zeta_measure / zeta)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

import numpy as np

from scipy import fftpack

from matplotlib import pyplot as plt

def damp_eq(m, k, zeta, x0, v0, t):

''' 減衰自由振動の理論値を計算 '''

omega = np.sqrt(k / m)

sigma = omega * zeta

sigma_v0 = v0 / (x0 * omega)

# 減衰比の場合分け

if zeta < 1:

# 減衰振動

omega_d = omega * np.sqrt(1 - zeta ** 2)

X = np.sqrt(x0 ** 2 + ((v0 + sigma * x0) / omega_d))

phi = np.arctan2((v0 + sigma * x0), (x0 * omega_d))

x = np.exp(-sigma * t) * X * np.cos(omega_d * t - phi)

elif zeta == 1:

# 臨界減衰

omega_d = omega

x = x0 * np.exp(-omega * t) * ((sigma_v0 + 1) * omega * t + 1)

else:

# 過減衰

omega_d = omega

x = x0 * np.exp(- zeta * omega * t) * (np.cosh(omega * t * np.sqrt(zeta ** 2 - 1))

+ (sigma_v0 + zeta) / (np.sqrt(zeta ** 2 + 1))

* np.sinh(omega * t * np.sqrt(zeta ** 2 - 1)))

return x, omega_d

def hilbert(yt, dt):

''' ヒルベルト変換 '''

# フーリエ変換

spectrum = fftpack.fft(yt)

# 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）

zeros = np.zeros(int(len(spectrum) / 2))

spectrum[int(len(spectrum) / 2):len(spectrum)] = zeros

# 正周波数域（ナイキスト周波数まで）を2倍する

spectrum *= 2

# 実部と虚部から瞬時振幅と瞬時位相を計算する

zt = fftpack.ifft(spectrum)

re = zt.real

im = zt.imag

# 瞬時振幅

instantaneous_amp = np.sqrt(re ** 2 + im ** 2)

# 瞬時位相とアンラップ処理

instantaneous_phase = np.arctan2(im, re)

instantaneous_phase = np.unwrap(instantaneous_phase)

# 瞬時周波数（位相角の微分と時間刻み影響を考慮）

instantaneous_frequency = ((1 / dt) / (2 * np.pi)) * np.gradient(instantaneous_phase)

return instantaneous_amp

def stable_check(t, y, i_amp, rate):

''' 安定して減衰振動している範囲を抽出する '''

# データの後半を切り捨て（包絡線の端部は瞬時振幅が高くなってしまうため）

t_half = t[0:int(len(t) / 2)]

y_half = y[0:int(len(y) / 2)]

# 最大値に対するrate割合を超える振幅成分のみ抽出する

peak = np.max(np.abs(y_half))

threshold = peak * rate

index_max = np.max(np.where(np.abs(y_half) > threshold))

# 安定して減衰振動している範囲を抽出

t_censored = t_half[0:index_max]

y_censored = y_half[0:index_max]

i_amp = i_amp[0:index_max]

plot([t, t_half, t_censored], [y, y_half, y_censored], ['original', 'half data', 'censored data'], 't[s]', 'Displacement[m]')

return t_censored, i_amp

def calc_damping_factor(t, y, omega_d):

''' 対数減衰率と減衰比を計算（ヒルベルト変換結果を用いて）'''

# デシベル変換（dBrefは何でも良い）

y_dB = 20 * np.log10(y / 1e-6)

# カーブフィット

coefficients = np.polyfit(t, y_dB, 1)

damp_curve = coefficients[0] * t + coefficients[1]

plot([t, t], [y_dB, damp_curve], ['envelope[dB]', 'curvefit'], 't[s]', 'Displacement[dB]')

# 対数減衰率deltaと減衰比zetaを計算

a = -coefficients[0]

fd = omega_d * (1 / (2 * np.pi))

delta = (a / fd) * (np.log(10) / 20)

zeta = delta / (2 * np.pi)

return zeta

def plot(t, x, label, xlabel, ylabel):

''' シミュレーション結果をプロットする '''

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(111)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel(xlabel)

ax1.set_ylabel(ylabel)

# スケールを設定する。

#ax1.set_xlim(0, 4)

# プロットを行う。

for i in range(len(x)):

ax1.plot(t[i], x[i], label=label[i], lw=1)

ax1.legend()

# レイアウト設定

fig.tight_layout()

# グラフを表示する。

plt.show()

plt.close()

return

if __name__ == '__main__':

''' メイン処理 '''

# サンプルの減衰振動を計算(m:質量[kg], k:ばね定数[N/m], zeta:減衰率, x0:初期変位[m], v0:初速度[m/s], dt:時間刻み[s])

m = 0.5

k = 1000

zeta = 0.02

x0 = 1

v0 = 0

dt = 0.001

t = np.arange(0, 10, dt)

x, omega_d = damp_eq(m, k, zeta, x0, v0, t)

plot([t], [x], ['Damping vibration'], 't[s]', 'Displacement[m]')

# ヒルベルト変換から包絡線（瞬時振幅）を計算

i_amp = hilbert(x, dt)

plot([t, t], [x, i_amp], ['disp', 'envelope'], 't[s]', 'Displacement[m]')

# 減衰係数を計算する範囲を自動抽出(rate:データを抽出する最大振幅に対する割合)

t, i_amp = stable_check(t, x, i_amp, 0.1)

# 減衰係数を計算

zeta_measure = calc_damping_factor(t, i_amp, omega_d)

# 理論値との精度を確認

print('Theoretical zeta =', zeta)

print('Measured zeta =', zeta_measure)

print('Accuracy =', zeta_measure / zeta)

最後に理論値として与えた減衰比と測定によって求められた減衰比の比較結果がコンソールに出力されます。Accuracyは1が完全一致なので、今回の結果は0.2[%]程度のずれと読み取れます。

Theoretical zeta = 0.02
Measured zeta = 0.019964782847878126
Accuracy = 0.9982391423939063

Theoretical zeta = 0.02

Measured zeta = 0.019964782847878126

Accuracy = 0.9982391423939063

精度検証：減衰比を振って理論値と比較する

最後に上記全コードを使って、減衰比を色々と振って理論値と比較してみましょう。手動でやるのはすごく面倒なので、上のコードのメイン部をforループで囲い、パラメータをリスト形式で与える等小変更して回してみます。

以下が検証結果です。理論波形で与えた減衰率と測定された減衰率は、 $ζ < 0.2$ あたりまでは2[%]程度のずれであるのに対し、 $ζ > 0.2$ は精度が落ちていることを確認しました。

以下は精度の悪かった $ζ = 0.3$ の場合の計算結果ですが、減衰比が高くなると安定して減衰振動する山の数がそもそも少なくなっています。

簡単な構造物でハンマリング試験を行うと金属材料で $ζ$ ≒0.01程度、ゴム材料で $ζ$ ≒0.1程度の印象です。実用には問題ないと思いますが、さらなる精度UPにはヒルベルト変換を用いた方法以外を考える必要がありそうです（理論式はわかっているので、時間波形に直接フィッティングさせる等）。

参考までに、精度検証に使用した全コードも記載しておきます。

import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt

def damp_eq(m, k, zeta, x0, v0, t):
    ''' 減衰自由振動の理論値を計算 '''

    omega = np.sqrt(k / m)
    sigma = omega * zeta
    sigma_v0 = v0 / (x0 * omega)

    # 減衰比の場合分け
    if zeta < 1:
        # 減衰振動
        omega_d = omega * np.sqrt(1 - zeta ** 2)
        X = np.sqrt(x0 ** 2 + ((v0 + sigma * x0) / omega_d))
        phi = np.arctan2((v0 + sigma * x0), (x0 * omega_d))
        x = np.exp(-sigma * t) * X * np.cos(omega_d * t - phi)
    elif zeta == 1:
        # 臨界減衰
        omega_d = omega
        x = x0 * np.exp(-omega * t) * ((sigma_v0 + 1) * omega * t + 1)
    else:
        # 過減衰
        omega_d = omega
        x = x0 * np.exp(- zeta * omega * t) * (np.cosh(omega * t * np.sqrt(zeta ** 2 - 1))
                                               + (sigma_v0 + zeta) / (np.sqrt(zeta ** 2 + 1))
                                               * np.sinh(omega * t * np.sqrt(zeta ** 2 - 1)))

    return x, omega_d

def hilbert(yt, dt):
    ''' ヒルベルト変換 '''

    # フーリエ変換
    spectrum = fftpack.fft(yt)

    # 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）
    zeros = np.zeros(int(len(spectrum) / 2))
    spectrum[int(len(spectrum) / 2):len(spectrum)] = zeros

    # 正周波数域（ナイキスト周波数まで）を2倍する
    spectrum *= 2

    # 実部と虚部から瞬時振幅と瞬時位相を計算する
    zt = fftpack.ifft(spectrum)
    re = zt.real
    im = zt.imag

    # 瞬時振幅
    instantaneous_amp = np.sqrt(re ** 2 + im ** 2)

    # 瞬時位相とアンラップ処理
    instantaneous_phase = np.arctan2(im, re)
    instantaneous_phase = np.unwrap(instantaneous_phase)

    # 瞬時周波数（位相角の微分と時間刻み影響を考慮）
    instantaneous_frequency = ((1 / dt) / (2 * np.pi)) * np.gradient(instantaneous_phase)

    return instantaneous_amp

def stable_check(t, y, i_amp, rate):
    ''' 安定して減衰振動している範囲を抽出する '''

    # データの後半を切り捨て（包絡線の端部は瞬時振幅が高くなってしまうため）
    t_half = t[0:int(len(t) / 2)]
    y_half = y[0:int(len(y) / 2)]

    # 最大値に対するrate割合を超える振幅成分のみ抽出する
    peak = np.max(np.abs(y_half))
    threshold = peak * rate
    index_max = np.max(np.where(np.abs(y_half) > threshold))

    # 安定して減衰振動している範囲を抽出
    t_censored = t_half[0:index_max]
    y_censored = y_half[0:index_max]
    i_amp = i_amp[0:index_max]
    #plot([t, t_half, t_censored], [y, y_half, y_censored], ['original', 'half data', 'censored data'], 't[s]', 'Displacement[m]')

    return t_censored, i_amp

def calc_damping_factor(t, y, omega_d):
    ''' 対数減衰率を計算（ヒルベルト変換結果を用いて）'''

    # デシベル変換（dBrefは何でも良い）
    y_dB = 20 * np.log10(y / 1e-6)

    # カーブフィット
    coefficients = np.polyfit(t, y_dB, 1)
    damp_curve = coefficients[0] * t + coefficients[1]
    #plot([t, t], [y_dB, damp_curve], ['envelope[dB]', 'curvefit'], 't[s]', 'Displacement[dB]')

    # 減衰係数zetaを計算
    dX = 1
    dY = -coefficients[0]
    fd = omega_d * (1 / (2 * np.pi))
    p = dX * fd
    delta = (1 / p) * np.abs(dY / dX) * (np.log(10) / 20)
    zeta = delta / (2 * np.pi)

    return zeta

def plot(t, x, label, xlabel, ylabel):
    ''' シミュレーション結果をプロットする '''

    # フォントの種類とサイズを設定する。
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    # 目盛を内側にする。
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    # グラフの上下左右に目盛線を付ける。
    fig = plt.figure()

    ax1 = fig.add_subplot(111)
    ax1.yaxis.set_ticks_position('both')
    ax1.xaxis.set_ticks_position('both')

    # 軸のラベルを設定する。
    ax1.set_xlabel(xlabel)
    ax1.set_ylabel(ylabel)

    # スケールを設定する。
    #ax1.set_xlim(0, 4)

    # プロットを行う。
    for i in range(len(x)):
        ax1.scatter(t[i], x[i], label=label[i], lw=1)
    ax1.grid()
    ax1.legend()

    # レイアウト設定
    fig.tight_layout()

    # グラフを表示する。
    plt.show()
    plt.close()

    return

if __name__ == '__main__':
    ''' メイン処理 '''
    zetas = np.arange(0.001, 0.5, 0.01)
    zetas_measured = []
    accuracies = []
    for i in range(len(zetas)):
        # サンプルの減衰振動を計算(m:質量[kg], k:ばね定数[N/m], zeta:減衰率, x0:初期変位[m], v0:初速度[m/s], dt:時間刻み[s])
        m = 0.5
        k = 1000
        zeta = zetas[i]
        x0 = 1
        v0 = 0
        dt = 0.001
        t = np.arange(0, 10, dt)
        x, omega_d = damp_eq(m, k, zeta, x0, v0, t)

        # ヒルベルト変換から包絡線（瞬時振幅）を計算
        i_amp = hilbert(x, dt)
        #plot([t, t], [x, i_amp], ['disp', 'envelope'], 't[s]', 'Displacement[m]')

        # 減衰係数を計算する範囲を自動抽出(rate:データを抽出する最大振幅に対する割合)
        t, i_amp = stable_check(t, x, i_amp, 0.1)

        # 減衰係数を計算
        zeta_measure = calc_damping_factor(t, i_amp, omega_d)
        zetas_measured.append(zeta_measure)

        # 精度を計算
        accuracy = 1 - (zeta_measure / zeta)
        accuracies.append(accuracy)

        # 理論値との精度を確認
        print('Theoretical zeta =', zeta)
        print('Measured zeta =', zeta_measure)
        print('Accuracy =', accuracy)

    plot([zetas], [zetas_measured], ['$\zeta$'], 'Teoretical $\zeta$', 'Measured $\zeta$')
    plot([zetas], [accuracies], ['Accuracy'], '$\zeta$', 'Accuracy')

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

import numpy as np

from scipy import fftpack

from matplotlib import pyplot as plt

def damp_eq(m, k, zeta, x0, v0, t):

''' 減衰自由振動の理論値を計算 '''

omega = np.sqrt(k / m)

sigma = omega * zeta

sigma_v0 = v0 / (x0 * omega)

# 減衰比の場合分け

if zeta < 1:

# 減衰振動

omega_d = omega * np.sqrt(1 - zeta ** 2)

X = np.sqrt(x0 ** 2 + ((v0 + sigma * x0) / omega_d))

phi = np.arctan2((v0 + sigma * x0), (x0 * omega_d))

x = np.exp(-sigma * t) * X * np.cos(omega_d * t - phi)

elif zeta == 1:

# 臨界減衰

omega_d = omega

x = x0 * np.exp(-omega * t) * ((sigma_v0 + 1) * omega * t + 1)

else:

# 過減衰

omega_d = omega

x = x0 * np.exp(- zeta * omega * t) * (np.cosh(omega * t * np.sqrt(zeta ** 2 - 1))

+ (sigma_v0 + zeta) / (np.sqrt(zeta ** 2 + 1))

* np.sinh(omega * t * np.sqrt(zeta ** 2 - 1)))

return x, omega_d

def hilbert(yt, dt):

''' ヒルベルト変換 '''

# フーリエ変換

spectrum = fftpack.fft(yt)

# 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）

zeros = np.zeros(int(len(spectrum) / 2))

spectrum[int(len(spectrum) / 2):len(spectrum)] = zeros

# 正周波数域（ナイキスト周波数まで）を2倍する

spectrum *= 2

# 実部と虚部から瞬時振幅と瞬時位相を計算する

zt = fftpack.ifft(spectrum)

re = zt.real

im = zt.imag

# 瞬時振幅

instantaneous_amp = np.sqrt(re ** 2 + im ** 2)

# 瞬時位相とアンラップ処理

instantaneous_phase = np.arctan2(im, re)

instantaneous_phase = np.unwrap(instantaneous_phase)

# 瞬時周波数（位相角の微分と時間刻み影響を考慮）

instantaneous_frequency = ((1 / dt) / (2 * np.pi)) * np.gradient(instantaneous_phase)

return instantaneous_amp

def stable_check(t, y, i_amp, rate):

''' 安定して減衰振動している範囲を抽出する '''

# データの後半を切り捨て（包絡線の端部は瞬時振幅が高くなってしまうため）

t_half = t[0:int(len(t) / 2)]

y_half = y[0:int(len(y) / 2)]

# 最大値に対するrate割合を超える振幅成分のみ抽出する

peak = np.max(np.abs(y_half))

threshold = peak * rate

index_max = np.max(np.where(np.abs(y_half) > threshold))

# 安定して減衰振動している範囲を抽出

t_censored = t_half[0:index_max]

y_censored = y_half[0:index_max]

i_amp = i_amp[0:index_max]

#plot([t, t_half, t_censored], [y, y_half, y_censored], ['original', 'half data', 'censored data'], 't[s]', 'Displacement[m]')

return t_censored, i_amp

def calc_damping_factor(t, y, omega_d):

''' 対数減衰率を計算（ヒルベルト変換結果を用いて）'''

# デシベル変換（dBrefは何でも良い）

y_dB = 20 * np.log10(y / 1e-6)

# カーブフィット

coefficients = np.polyfit(t, y_dB, 1)

damp_curve = coefficients[0] * t + coefficients[1]

#plot([t, t], [y_dB, damp_curve], ['envelope[dB]', 'curvefit'], 't[s]', 'Displacement[dB]')

# 減衰係数zetaを計算

dX = 1

dY = -coefficients[0]

fd = omega_d * (1 / (2 * np.pi))

p = dX * fd

delta = (1 / p) * np.abs(dY / dX) * (np.log(10) / 20)

zeta = delta / (2 * np.pi)

return zeta

def plot(t, x, label, xlabel, ylabel):

''' シミュレーション結果をプロットする '''

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(111)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel(xlabel)

ax1.set_ylabel(ylabel)

# スケールを設定する。

#ax1.set_xlim(0, 4)

# プロットを行う。

for i in range(len(x)):

ax1.scatter(t[i], x[i], label=label[i], lw=1)

ax1.grid()

ax1.legend()

# レイアウト設定

fig.tight_layout()

# グラフを表示する。

plt.show()

plt.close()

return

if __name__ == '__main__':

''' メイン処理 '''

zetas = np.arange(0.001, 0.5, 0.01)

zetas_measured = []

accuracies = []

for i in range(len(zetas)):

# サンプルの減衰振動を計算(m:質量[kg], k:ばね定数[N/m], zeta:減衰率, x0:初期変位[m], v0:初速度[m/s], dt:時間刻み[s])

m = 0.5

k = 1000

zeta = zetas[i]

x0 = 1

v0 = 0

dt = 0.001

t = np.arange(0, 10, dt)

x, omega_d = damp_eq(m, k, zeta, x0, v0, t)

# ヒルベルト変換から包絡線（瞬時振幅）を計算

i_amp = hilbert(x, dt)

#plot([t, t], [x, i_amp], ['disp', 'envelope'], 't[s]', 'Displacement[m]')

# 減衰係数を計算する範囲を自動抽出(rate:データを抽出する最大振幅に対する割合)

t, i_amp = stable_check(t, x, i_amp, 0.1)

# 減衰係数を計算

zeta_measure = calc_damping_factor(t, i_amp, omega_d)

zetas_measured.append(zeta_measure)

# 精度を計算

accuracy = 1 - (zeta_measure / zeta)

accuracies.append(accuracy)

# 理論値との精度を確認

print('Theoretical zeta =', zeta)

print('Measured zeta =', zeta_measure)

print('Accuracy =', accuracy)

plot([zetas], [zetas_measured], ['

ζ

'], 'Teoretical

ζ

', 'Measured

ζ

plot([zetas], [accuracies], ['Accuracy'], '

ζ

', 'Accuracy')

まとめ

本記事では減衰自由振動の対数減衰率と減衰比をPythonで測定する方法を学びました。

式自体は比較的有名ですが、実際に測定しようと思うとデータ処理を色々やらないと精度が著しく悪化するということもわかりました。

今回は理想波形における測定を行いましたが、次回は実際の音振波形に対して適用する方法を検討します。

おそらくフィルタ処理等を色々しないと実問題には適用できないと考えられるので、この記事は第1歩といったところ。

参考文献

[1]:金子成彦, 大熊政明, 機械力学ハンドブック, 朝倉書店, 2015(初版第1刷), pp87

[2]: 小野測器, 減衰をあらわす係数の意味と求め方

[3]:Pythonで1自由度減衰系の自由振動シミュレーション

[4]:Wikipedia, 減衰振動

[5]:Pythonでヒルベルト変換！時間波形の包絡線を求める方法

[6]:Pythonでカーブフィット!最小二乗法で直線近似する方法

これで減衰特性を測定できるようになりました！また、ヒルベルト変換も十分使いこなせるようになった気がします！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

振動波形のヒルベルト変換から対数減衰率を求めるPythonコード

減衰自由振動のおさらい

減衰自由振動とは？

対数減衰率とは？

対数減衰率と減衰比の関係

振動波形から対数減衰率と減衰比を求めるPythonコード

動作環境

サンプルの減衰自由振動波形を作る

理論式： $0 < ζ < 1$

理論式： $ζ = 1$

理論式： $ζ > 1$

理論波形作成を関数にまとめたPythonコード

ヒルベルト変換で瞬時振幅を求める

ヒルベルト変換の考え方

ヒルベルト変換を行うPythonコード

解析範囲を自動調整する

そのままカーブフィットすると精度が落ちる

解析範囲を自動調整するPythonコード

対数減衰率と減衰比を計算する

振幅をデシベルに変換する

カーブフィットする

対数減衰率と減衰比を計算する

対数減衰率と減衰比を求めるPythonコード

全コード（コピペ可：プロット含む）

精度検証：減衰比を振って理論値と比較する

まとめ

参考文献

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

減衰自由振動のおさらい

減衰自由振動とは？

対数減衰率とは？

対数減衰率と減衰比の関係

振動波形から対数減衰率と減衰比を求めるPythonコード

動作環境

サンプルの減衰自由振動波形を作る

理論式：0<ζ<1

理論式：ζ=1

理論式：ζ>1

理論波形作成を関数にまとめたPythonコード

ヒルベルト変換で瞬時振幅を求める

ヒルベルト変換の考え方

ヒルベルト変換を行うPythonコード

解析範囲を自動調整する

そのままカーブフィットすると精度が落ちる

解析範囲を自動調整するPythonコード

対数減衰率と減衰比を計算する

振幅をデシベルに変換する

カーブフィットする

対数減衰率と減衰比を計算する

対数減衰率と減衰比を求めるPythonコード

全コード（コピペ可：プロット含む）

精度検証：減衰比を振って理論値と比較する

まとめ

参考文献

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

理論式： $0 < ζ < 1$

理論式： $ζ = 1$

理論式： $ζ > 1$

コメントを残すコメントをキャンセル