Pythonで伝達関数!2つの信号の周波数応答関数を計算する方法

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

信号処理の分野では2つの信号の関係を周波数領域で分析するために周波数応答関数という計算を多用します。Pythonでフーリエ変換が出来るようになったらあと一歩なので、ここではPythonによる周波数応答関数の計算方法を紹介します。

こんにちは。wat(@watlablog)です。ここではPythonを使って信号処理の重要計算である周波数応答関数を求めてみます

周波数応答関数の概要

伝達関数?周波数応答関数?

機械系でも電気系でも、実物でもシミュレーションモデルでも、開発系の仕事をしている人であればあるシステムに何らかの信号を入力したときの出力を分析することが多いと思います。

Advertisements

僕は機械系なので構造物で例えますが、以下のようにある物体を任意の力[\(\rm N\)]で加振させ、その物体の加速度応答[\(\rm m/s^{2}\)]との比を評価することがよくあります。

入力と出力の例
Advertisements

ここで、力や加速度、他にも速度や変位、ひずみや応力といった物体の状態を決定することができる量状態量(State value)と呼び、これらの状態量を使って入力と出力の比をとったものを伝達関数(Transfer function)と呼びます。

この伝達関数の独立変数に周波数をとったものを特に周波数応答関数(FRF:Frequency Response Function)と呼びます。

実際に使う時は、単にFRFと呼ぶことが多いですね。

この周波数応答関数からは様々な情報を得ることが出来ます。ではその目的を探ってみましょう。

FRFは入力と出力の関係を分析するために計算する

先ほど図示した、入力である加振力と応答である加速度の波形は、波の形も異なるので、一見するだけでは2信号間にどのような関係があるのかを調べることは困難です。

振動している信号の分析としてはフーリエ変換を使って周波数振幅位相を調べることがよくありますが、それぞれの信号を別々にフーリエ変換した所で、2信号間の関係性を分析することは依然として困難なままです。

そこで、FRFの出番となります。

FRFを計算することで、
入力の振幅に対して応答の振幅がどれだけ出たか?
入力に対して応答はどれだけ遅れて発生したか?
…といった情報が一目瞭然となります。

さらに、それらの情報からシステムがどのような振動モードを持っているかを同定する実験モーダル解析が出来たり、プロットを工夫して共振点の評価にどれだけ誤差を含んでいるかといった分析も可能となります。

いろんな物理分野の分析に大変重宝するため、是非覚えてみましょう!

周波数応答関数の式

先ほど文章で説明した通り、入力\(I\)と出力\(O\)の比である伝達関数\(H\)で周波数\(f\)(または角振動数\(\omega \))を考慮すれば良いので、FRFは以下の式(1)で定義されます。

\[
H(\omega )=\frac{O(\omega )}{I(\omega )} \tag{1}
\]

ここで、\(I(\omega ), O(\omega )\)は共に複素数です。
一般に入力と出力の各信号は共に位相情報も周波数の関数であるので、式(2)と分母を実数にするために共役複素数\(I^{\star}(\omega )\)(複素数の虚部に-1を乗じたもの、1+2jであれば共役複素数は1-2j)を分母分子に乗じます。

\[
H(\omega )=\frac{O(\omega )I^{\star}(\omega )}{I(\omega )I^{\star}(\omega )} \tag{2}
\]

ちなみに、分母の\(I(\omega )I^{\star}(\omega )\)はパワースペクトル(Power spectrum)、分子の\(O(\omega )I^{\star}(\omega )\)はクロススペクトル(Cross spectrum)と呼ばれます。

クロススペクトル単体は2信号間の相関や位相差を分析する時に使われますが、FRFの振幅成分は入出力の比、位相成分は入出力の差を意味する関数であるとわかります。

周波数応答関数の代表例

機械系でよく使われるFRFの呼称と単位を下表に示します。参考までに。

呼称単位
コンプライアンス(Compliance)[\(\rm m/N\)]
モビリティ(Mobility)[\(\rm m/(Ns)\)]
アクセレランス(Accelerance)[\(\rm m/(Ns^{2})\)]
動剛性(Dynamic stiffness)[\(\rm N/m\)]
機械インピーダンス(Mechanical impedance)[\(\rm Ns/m\)]
動質量(Apparent mass)[\(\rm Ns^{2}/m\)]

参考文献

周波数応答関数については、以下の「モード解析入門」に式(P198くらい)や理論的背景が詳しく書いてあります。本書は振動学を学ぶ人必携の書と思いますので、是非書店等で見てみて下さい。

モード解析手法の決定版。1自由度の振動から多自由度、モードの理論と実験的手法まで網羅された良書。

さて、式が出て概要がわかってきたので、早速Pythonによるコーディングをしていきましょう!

PythonでFRFを計算するコード

高速フーリエ変換(FFT)についてのおさらい

今回はフーリエ変換を既にPythonで使えるものとして話を進めますが、フーリエ変換の概要、計算上の注意点、Pythonコードについては過去に書いた「PythonでFFT!SciPyのFFTまとめ」をご覧下さい。

しかし上記記事はやや複雑な処理を複数使っています。本記事では窓関数やオーバーラップ処理は使わず、最も簡易なフーリエ変換を使ってコーディングします。

上記記事の説明は前提として、コード自体は「PythonでFFTとIFFT!逆フーリエ変換で時間波形を作る」で紹介したものを流用したいと思います。

準備:サンプルの信号を用意する

FRFの計算確認で適当に用意した信号を使うとあまり面白くありません。

ここでは出来るだけ物理的に値を確認できるよう、「工学カテゴリ」で学んだ振動モデルを使います。

シミュレーションモデル

モデルは以下に示す2自由度のばねマスダンパー系を用意しました。強制加振を入力、\(x_{2}\)の変位応答をとってFRFを計算します。

2自由度ばねマスダンパー系

このモデルの運動方程式は式(3)です。全て行列、ベクトル形式で書いています。

\[
\mathbf{M}\mathbf{\ddot{x}}+\mathbf{C}\mathbf{\dot{x}}+\mathbf{K}\mathbf{x}=\mathbf{F} \tag{3}
\]

数値計算コードを別途.pyファイルで作成しておく

まずは上記運動方程式を数値解析するコードを「simulation.py」というファイルで作成します。

このコードは「Pythonの4次ルンゲ・クッタ法で多自由度連成振動を解く方法」でやった時とほとんど同じですが、外力をランダム加振とするために予め外力の離散1D配列を作っておく所、メインの計算実行文もdef関数にしている所が異なります。

メイン:PythonでFRFを計算するコード

ボード線図を描いてみる

以下がメインの.pyファイルの内容です(グラフ化まで含んでいるので少々長いですが)。こちらのメインファイルの方には、FRF計算部分のみを関数にし、先ほどのsimulation.pyをimportしてサンプル信号を生成、その信号を使ってFRFを算出するという流れです。

以下が結果です。左の2プロットが入力としてのランダム加振力、変位の時間波形です。そして右の2プロットが周波数応答関数の振幅(下)と位相(上)です。

振幅と位相を縦にプロットしたグラフをボード線図(Bode plot)と呼び、一般的に広く使われている描画方法です。

ボード線図の計算結果

2自由度なので、振幅にはピークが2つ確認出来ます。そして振幅でピークが見える周波数には位相が90[deg.](力と変位の関係なので)になっています。

このくらいの検証しかしていませんが、モード理論通りなのでうまくいっていると考えます。

コクアド線図を描いてみる

以上でFRFの目的はおおよそ達成したのですが、他にも有用なプロット方法があるため、Pythonで再現してみたいと思います。最初はコクアド線図です。

コクアド線図は共振のピークがボード線図よりも急峻であることから、共振判別の1つの方法として使われます(あまり一般的には見ませんが)。

この線図は横軸に周波数をとり、縦軸に実部と虚部を並べます。英語でCoincident quadrature plotとなるため、Co. quad. plotとなります。

コードはグラフプロット部分のみを以下のように変更します。FRF計算関数には予め実部と虚部が入ったh_ioも出力しているのでこの部分だけでOKです。

以下がコクアド線図の結果です。実部(Re)と虚部(Im)で急峻なピークが確認出来ました。

コクアド線図の結果

ナイキスト線図を描いてみる

次はナイキスト線図(Nyquist plot)です。このナイキスト線図は横軸に実部、縦軸に虚部をとってプロットする図で、共振部分だけにフォーカスした図になり、1つの図だけ見れば良いというメリットがあるのでコクアド線図よりは使われているイメージです。

以下がナイキスト線図のプロット結果です。実部と虚部の値が共振部分で円を描くというプロットです…が、シミュレーションと言ってもランダムノイズを加振に使用しているからか、ちょっといびつな形状をしていますね。ナイキスト線図は円の綺麗さから非線形性や誤差の割合を俯瞰できるため、実験でよく使われていますかね。

ナイキスト線図の結果

まとめ

本記事では周波数応答関数(FRF)の概要を説明し、定義式や機械系で代表的なFRFの種類を紹介しました。

また、PythonによるFRF算出コードを紹介し、実際にシミュレーションモデルを使った事例としてボード線図、コクアド線図、ナイキスト線図をそれぞれ紹介しました。

FRFは実験やシミュレーションで分野を問わず頻繁に使われる強力な信号分析手法です。是非本ページのコードを参考にしてみて下さい。

FRFを覚えたので周波数分析の幅がひろがりました!Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

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

SNSでもご購読できます。

コメント

コメントを残す

*