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

  • このエントリーをはてなブックマークに追加
Advertisements

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

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

減衰自由振動のおさらい

減衰自由振動とは?

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

減衰自由振動の例
Advertisements

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

当ブログでは減衰自由振動のシミュレーションを「Pythonで1自由度減衰系の自由振動シミュレーション」という記事で紹介しました。運動方程式や理論値も説明していますので是非ご確認ください。

対数減衰率とは?

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

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

対数減衰率の図

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

\[ \delta = \ln \frac{a_{n}}{a_{n+1}} \tag{1} \]

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

\[ \delta = \frac{1}{m}\ln \frac{a_{n}}{a_{n+m}} \tag{2} \]

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

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

\[ \zeta = \frac{\delta}{2\pi} \tag{3} \]

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

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

また、上式(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 < \zeta < 1\)

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

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

\[ \omega= \sqrt{{\frac{k}{m}}}  \tag{4} \]
\[ \omega _{d}= \omega \sqrt{1-\zeta ^{2}} \tag{5} \]
\[ \sigma = \omega \zeta \tag{6} \]
\[ X= \sqrt{x_{0}^{2}+\left ( \frac{v_{0}+\sigma x_{0}}{\omega _{d}} \right )^{2}} \tag{7} \]
\[ \phi = \tan^{-1}\left ( \frac{v_{0}+\sigma x_{0}}{x_{0}\omega _{d}} \right ) \tag{8} \]
\[ x= e^{-\sigma t}X \cos\left ( \omega _{d}t-\phi \right )  \tag{9} \]

理論式:\(\zeta = 1\)

同様に、式(10)の無次元化速度\(\sigma_{v0}\)を導入することで\(\zeta=1\)の場合の応答変位は式(11)で計算可能です\(^{[4]}\)。

\[ \sigma_{v0} = \frac{v_{0}}{x_{0} \omega} \tag{10} \]
\[ x=x_{0}e^{- \omega t} \left \{ (\sigma_{v0} + 1) \omega t + 1 \right \} \tag{11} \]

理論式:\(\zeta > 1\)

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

\[ x=x_{0}e^{-\zeta \omega t}\left \{\cosh \left ( \omega t \sqrt{\zeta^{2}-1}\right ) + \frac{\sigma_{v0}+\zeta}{\sqrt{\zeta^{2}-1}}\sinh \left ( \omega t \sqrt{\zeta^2-1}\right ) \right \} \tag{12} \]

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

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

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

理論式から求めた減衰自由振動波形

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

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

ヒルベルト変換の考え方

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

実信号と複素信号

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

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

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

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

減衰振動の包絡線を求めた結果

解析範囲を自動調整する

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

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

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

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

カーブフィット時の問題点

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

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

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

データ打ち切り

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

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

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

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

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

\[ x_{\mathrm{dB}} = 20 \log_{10} \frac{x}{\mathrm{dBref}} \tag{13} \]

カーブフィットする

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

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

カーブフィットした結果

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

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

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

\[ f_{d} = \frac{1}{2\pi} \omega_{d} \tag{14} \]
\[ \delta = \frac{a}{f_{d}}\frac{\ln 10}{20} \tag{15} \]

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

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

全コード(コピペ可:プロット含む)

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

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

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

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

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

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

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

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

まとめ

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

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

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

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

参考文献

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

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

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

[4]:Wikipedia, 減衰振動

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

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

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

SNSでもご購読できます。

コメントを残す

*