Python NumPy SciPy : 周波数応答と伝達関数

何回かに渡って FFT 処理の基本をまとめてきました。
今回は周波数応答と伝達関数を求めてボード線図を書く基本的な方法について説明します。

測定データから周波数応答を求める

電気や機械の設計者は、信号を FFT 処理して周波数特性を解析するだけでなく、入力信号に対する出力信号の応答を解析することがあると思います。
機械屋さんだと振動試験で加振点近傍の加速度と任意点の加速度を測定し、周波数応答から構造物の固有振動数や応答倍率を調べたりします。
周波数応答はデータロガーで計算してくれる場合も多いですが、自分でいろいろと信号処理して周波数特性を調べたいことがあります。
周波数応答は以下のように計測データをそれぞれ FFT 処理し、 差分をとる 応答/入力で求めることができます。
ただ、実際にはパワースペクトルとクロススペクトルから求めることが多いです(こちら)
この辺りは 小野測器さんのページ などが参考になります。また、フィルタ処理や窓処理は別途必要ですし、要すれば平均化処理なども行なう必要があります。

# coding: utf-8
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
from scipy import signal

# # 時系列のサンプルデータ作成
n = 2048                         # データ数
dt = 0.005                       # サンプリング間隔
t = np.linspace(1, n, n)*dt-dt

# パラメータ設定
m = 1
c = 0.01
k = 400

# デルタ関数
yIn = np.zeros_like(t)
yIn[0] = k

# インパルス応答
num = [k]
den = [m, c, k]
s1 = signal.lti(num, den)
T, yOut = signal.impulse(s1, T=t)

# FFT 処理と周波数スケールの作成
yfIn = fftpack.fft(yIn)/(n/2)
yfOut = fftpack.fft(yOut)/(n/2)
freq = fftpack.fftfreq(n, dt)

# 周波数応答
# FRF = yfOut/yfIn
FRF = (yfOut*np.conj(yfIn))/(yfIn*np.conj(yfIn))
# FRF = (yfOut*np.conj(yfOut))/(yfIn*np.conj(yfOut))

# プロット
plt.figure()
plt.subplot(211)
plt.loglog(freq[1:n/2], np.abs(FRF[1:n/2]))
plt.ylabel("Amplitude")
plt.axis("tight")
plt.subplot(212)
plt.semilogx(freq[1:n/2], np.degrees(np.angle(FRF[1:n/2])))
plt.xlabel("Frequency[Hz]")
plt.ylabel("Phase[deg]")
plt.axis("tight")
plt.ylim(-180, 180)
plt.savefig('../files/150612FRF01.svg')
# plt.show()

Sorry, your browser does not support SVG.

伝達関数

設計においては系の方程式から伝達関数を求め、周波数特性を調べることも重要です。
ここでは一自由度粘性減衰系に外から加振力が作用するときの応答を例に、Python で伝達関数のボード線図を書く方法を説明します。

一自由度粘性減衰系の定常応答

簡単に式変形を載せますが、詳細は振動工学の教科書を参照してください。以下の運動方程式において

15061201.png

励振力を $f(t)=Fe^{jω t}$、これを満足する解を \(x(t)=Xe{j\omega t}\) と仮定して代入し、

15061202.png

不減衰固有角振動数 $Ω$、臨界減衰比 $c_c$、減衰比 \(\zeta\) を用いて

15061203.png

式を整理すると

15061204.png

となり、\(\beta=\frac{\omega}{\Omega}\) として伝達関数は以下のように求まります。

15061205.png

これを Python でプロットすると以下のようになります。

# coding: utf-8
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt

# 時系列のサンプルデータ作成
n = 1024                         # データ数
dt = 0.001                       # サンプリング間隔

# パラメータ設定
m = 1
c = 1
k = 400

freq = fftpack.fftfreq(n, dt)
zeta = c/(2*np.sqrt(m*k))
omega = np.sqrt(k/m)
beta = freq/omega

# 伝達関数
TF = 1/(1-beta**2+2*1j*zeta*beta)

# プロット
plt.figure(1)
plt.subplot(211)
plt.loglog(freq[1:n/2], np.abs(TF[1:n/2]))
plt.ylabel("Amplitude")
plt.axis("tight")
plt.subplot(212)
plt.semilogx(freq[1:n/2], np.angle(TF[1:n/2])*180/np.pi)
plt.xlabel("Frequency[Hz]")
plt.ylabel("Phase[deg]")
plt.axis("tight")
plt.ylim(-180, 180)
plt.savefig('../files/150612TF01.svg')
# plt.show()

Sorry, your browser does not support SVG.

ラプラス変換

一自由度系なら簡単ですが、もっと複雑な方程式で伝達関数を求めるのは大変です。そこで便利なのがラプラス変換です。

同じ一自由度粘性減衰系を例にします。運動方程式の各項はラプラス変換すると以下のようになります。なお、\(x_0=\dot{x}_0=0\) とします。

15061206.png

よって方程式は \((ms^{2}+cs+k)X(s)=F(s)\) となり、\(X_{st}(s)=\frac{F(s)}{k}\) と置けばラプラス変換後の伝達関数は

15061207.png

と求まります。

SciPy の lti, bode 関数を使えば、伝達関数の分子、分母の係数を入力するだけで以下のようにボード線図が書けます。

bode 関数は Matlab だと Control toolbox を購入しないと使えませんが、Python ならフリーで同様のことができてしまいます。

# coding: utf-8
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# パラメータ設定
m = 1
c = 1
k = 400

num = [k]
den = [m, c, k]
s1 = signal.lti(num, den)
w, mag, phase = signal.bode(s1, np.arange(1, 500, 1))

# プロット
plt.figure(1)
plt.subplot(211)
plt.loglog(w, 10**(mag/20))
plt.ylabel("Amplitude")
plt.axis("tight")
plt.subplot(212)
plt.semilogx(w, phase)
plt.xlabel("Frequency[Hz]")
plt.ylabel("Phase[deg]")
plt.axis("tight")
plt.ylim(-180, 180)
plt.savefig('../files/150612TF02.svg')
# plt.show()

Sorry, your browser does not support SVG.

コメント

Comments powered by Disqus
書籍更新情報
2016-10-21
Pythonによる科学技術計算 基礎編
PDF版の販売を開始しました。
販売ページはこちら

2016-09-09
Pythonによる科学技術計算 基礎編
1.2版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストア、Leanpubで販売中です
Pythonによる科学技術計算 基礎編
PDF版の販売はこちら
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share