Pythonでコヒーレンス関数!FRFのノイズ成分を評価する方法

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

周波数応答関数(FRF)はしばしばノイズに悩まされます。正確な評価のためにはノイズを正しく見積もる必要があり、その時はコヒーレンス関数が有用です。ここではPythonでコヒーレンス関数を計算する方法と簡単な分析例を紹介します。

こんにちは。wat(@watlablog)です。ここではFRF分析とセットで使われるコヒーレンス関数の計算方法と計算例を紹介します

コヒーレンス関数の概要

周波数応答関数のおさらい

本記事ではコヒーレンス関数の使い方を習得し、周波数応答関数のノイズを評価することを目標としますが、まずは周波数応答関数(FRF)について知っておく必要があります。

FRFは以下の式(1)で計算することができます。ここで分母がパワースペクトル、分子がクロススペクトルであることは「Pythonで伝達関数!2つの信号の周波数応答関数を計算する方法」で説明した通りです。

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

パワースペクトルを\(P\)、クロススペクトルを\(C\)とし、入力\(i\)と出力\(o\)の関係を添字で表現すると式(2)と書くことも可能です(複素共役で演算することを知っているのが前提ですが、こちらの記法の方がシンプル)。

\[
H(\omega )=\frac{C_{io}(\omega )}{P_{ii}(\omega )} \tag{2}
\]

理想的な条件下であればここまでの知識で2信号の関係を分析することに対して問題はありません。

しかし、リアルワールドデータとしての出力信号には往々にしてノイズが乗ってしまいます。より正確な分析を行うためには、このノイズ成分を評価することが重要です。

ノイズが乗るとどうなる?

式(1)と式(2)の状況は以下の信号系モデルで表現することが可能です。このモデルは任意の伝達系\(h(t)\)に入力\(i(t)\)が入ってきた時に、出力\(o(t)\)を出すという線形モデルです(\(t\)は時間を表します)。

理想的な系

線形で、何のロス(漏れ)もなく、出力に遅延もなく、かつ外乱も入ってこない理想的な状況(理論式ベースのシミュレーション等)では問題無いものの、この状況が1つでも崩れることでそれは出力信号にノイズとなって現れてしまいます。

百聞は一見にしかず…ということで、実際にノイズの有無でFRFがどうなるかを見てみましょう!

理想的なFRFの例

まずは理想的なデータで計算したFRF(ボード線図)の例を以下の図に示します。

入力データはインパルス波形を模擬したForce[N], 出力波形は2自由度ばねマスダンパー系を4次のRunge-Kutta法で解いたシミュレーション結果(Displacement[m])です。
※理想的とは言っても数値解析なので、シミュレーション上の誤差は乗っていますが…。

理想的なデータで計算したFRF

振幅波形(右下)にはピークが2つ、2自由度分綺麗に見えています。これは想定通り。

ノイズが乗ってしまったFRFの例

続いてノイズを重畳させてみましょう。ノイズ成分は出力波形にガウシアンノイズを加算しています。下図がノイズの乗った出力信号で計算したFRFです。

出力信号にノイズが乗ったデータで計算したFRF

振幅波形には…もはや2つ目のピークは目視で確認できず位相波形(右上)もひどい有様です。

このように、ノイズが乗ってしまうと本来あるはずの情報を正確に分析できないといった大きな問題が生じてしまいます。

ノイズを考慮したモデルとコヒーレンス関数

モデルとノイズ評価の意義

ノイズは様々な状況で信号に加算されます。

実験では測定器の測定範囲に対する計測信号の比が悪いことや、電源ノイズといった外乱、シミュレーションでは時間刻みやアルゴリズム、収束性といった所に原因がある場合が多いでしょうか?

以下の図は出力信号にノイズが重畳する場合の信号系モデルです。先ほどの理想的な系の場合と比べ、伝達系\(h(t)\)から出てきた純粋な出力\(v(t)\)に外乱であるノイズ\(n(t)\)が加算されて最終的な出力\(o(t)\)が得られるという所に違いがあります。

ここで、通常観測可能な物理量としては入力\(i(t)\)と出力\(o(t)\)のみです。周波数応答関数として伝達系を推定しようとしても、ノイズ\(n(t)\)に邪魔をされて純粋な出力\(v(t)\)を抽出することが困難です。

そのため、一般的にはノイズが出来るだけ入らないようにノイズを評価しながら分析していくことが重要とされています。

ノイズが重畳された系

コヒーレンス関数(関連度関数)の式

それではコヒーレンス関数(Coherence function:※関連度関数とも呼ばれます。)の式を見ていきましょう。

コヒーレンス関数\(\gamma _{io}^{2}(f) \)は式(3)で計算することが出来ます(\(f\)は周波数)。計算はクロススペクトルの絶対値の二乗を入力と出力それぞれのパワースペクトルで除するだけなので、既にFRFを計算したことがあれば簡単ですね。

\[
\gamma _{io}^{2}(f) = \frac{ |C_{io}(f)|^{2}}{P_{ii}(f) P_{oo}(f)}\tag{3}
\]

\[
0 \leq \gamma _{io}^{2}(f) \leq 1
\]

コヒーレンス関数は周波数軸で0から1の間に値をとります。1に近ければコヒーレンス、すなわち入力信号と出力信号の関連度が高いということになります。

式の導出の詳細は「モード解析入門」のP224辺りに載っています。意味合いとしては、周波数領域で定義された相関係数の二乗となるそうです。

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

コヒーレンス関数は平均化しないと意味がない

式(3)でコヒーレンス関数の計算方法はわかったので、後は実装すれば良いだけですが、1つ注意点があります。

それは、コヒーレンス関数はクロススペクトルとパワースペクトルの平均化をしないと意味がないということです(上記参考書にはP322に記載がありますね)。

この後紹介するPythonコードで計算を実行してみればわかりますが、このコヒーレンス関数は2つの関数がどれくらい似ているのかではなく、再現性があるかどうかを評価する関数です。そのため、複数回ノイズ成分を含んだ計測等を行い、平均化させる必要があります。

Pythonでコヒーレンス関数を計算するコード

Advertisements

サンプルの信号を別.pyファイルで準備する

周波数応答関数の時と同じですが、入力と出力の信号は「simulation.py」を別途作成し、シミュレーションによって用意します。

前回はランダムノイズを加振力として設定していましたが、今回はインパルス波形にしているため再掲します。

メインコード

import文

続いてメインコードです。import文は定番のnumpy, フーリエ変換用のscipy, グラフ描画用のmatplotlibの他に、先ほど作成したsimulation.pyを読むためのsimulationを記述します。

フーリエ変換する関数

今回は平均化処理を行うため、「Pythonで伝達関数!2つの信号の周波数応答関数を計算する方法」とは異なるコード構成をとっています。以下のようにフーリエ変換をするコードは別途作成しておきます。

この関数はスペクトルを得るためのinput, outputの時間波形と周波数軸を計算するための時間刻みdtを入力して使います。

FRFとコヒーレンス関数を計算する関数

周波数応答関数(FRF)とコヒーレンス関数は同じ関数の中で計算します。
その理由は、平均化を考慮しているからです。

この関数は複数回分のクロススペクトルとパワースペクトルを関数の外でも中でも保持できるように、入力と出力の両方に「〜_list」という変数を設定しています(シフトレジスタ的にしたかった)。

ボード線図描画のための関数

ボード線図のための振幅と位相の計算は関数で用意しておいた方が、今回の平均化FRFの場合は便利だったので、以下のコードのように準備しておきます。

各種関数の実行とグラフ描画

最後は各種関数を実行するための文とグラフ描画のための文を書いておしまいです。forループを使ってシミュレーションを3回計算し、毎回新しいノイズを重畳してFRFを計算していくという内容です。

もしこのプログラムを測定用に使う場合は、simulation()の部分をデータ収集コードに換えることで実現できると思います。

実行結果

このコードを実行することで、以下の図のような結果を得ます。左から時間波形、FRF(ボード線図)、コヒーレンス関数です。

実行結果

ノイズ違いで比較してみる

ガウシアンノイズの振幅を調整し、ノイズの量でコヒーレンス関数の値がどう変わるのかを下図で比較してみました。

今回は2自由度の振動シミュレーションを信号として使っているので、FRFには振動モードが存在することが特徴です。

(a)ノイズが多い場合は(b)ノイズが少ない場合に比べて全体的にコヒーレンス値が低い傾向にありますが、モード部分は比較的高い値を維持しようとしている特徴がありました。

コヒーレンス関数とノイズの関係

実験であれば測定毎に毎回コヒーレンスを評価し、著しく値が低下しないかどうかを確認することが良いとされています。
(モード周波数の部分がよければ後はOKと妥協することもありますが…。)

メインコード全文(コピペ用)

コピペしやすいようにメインコード部分全コードを以下に示します。
是非信号や平均化回数と条件を変えて遊んでみて下さい。

まとめ

本記事では周波数応答関数(FRF)の評価に必要不可欠なコヒーレンス関数の概要を説明し、Pythonによる実装方法を紹介しました。

コヒーレンス関数は平均化処理が必須であり、信号同士の再現性を分析、ノイズ成分の大小を評価することができます。

ノイズ評価は様々な工学的場面で重要となりますので、是非活用してみて下さい。

FRFとともにコヒーレンス関数のPythonコードを習得しました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*