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.show()
伝達関数
設計においては系の方程式から伝達関数を求め、周波数特性を調べることも重要です。 ここでは一自由度粘性減衰系に外から加振力が作用するときの応答を例に、Python で伝達関数のボード線図を書く方法を説明します。
一自由度粘性減衰系の定常応答
簡単に式変形を載せますが、詳細は振動工学の教科書を参照してください。以下の運動方程式において
\begin{equation} m\ddot{x}+c\dot{x}+kx=f(t) \nonumber \end{equation}励振力を $f(t)=Fe^{jω t}$、これを満足する解を \(x(t)=Xe{j\omega t}\) と仮定して代入し、
\begin{equation} (-\omega^2 m+j\omega c+k)Xe^{j\omega t}=Fe^{j\omega t} \nonumber \end{equation}不減衰固有角振動数 $Ω$、臨界減衰比 $c_c$、減衰比 \(\zeta\) を用いて
\begin{equation} \frac{m}{k}=\frac{1}{\Omega},\ \frac{c}{k}=\frac{2c}{2\sqrt{mk}}\sqrt{\frac{m}{k}}=2\frac{c}{c_c}\frac{1}{\Omega}=2\zeta\frac{1}{\Omega} \nonumber \end{equation}式を整理すると
\begin{equation} \Bigl(-\left(\frac{\omega}{\Omega}\right)^2 + 2j\zeta\frac{\omega}{\Omega} + 1\Bigr)X=\frac{F}{k} \nonumber \end{equation}となり、\(\beta=\frac{\omega}{\Omega}\) として伝達関数は以下のように求まります。
\begin{equation} \frac{X}{X_{st}}=\frac{1}{1-\beta^2 + 2j\zeta\beta} \nonumber \end{equation}これを Python でプロットすると以下のようになります。
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.show()
ラプラス変換
一自由度系なら簡単ですが、もっと複雑な方程式で伝達関数を求めるのは大変です。そこで便利なのがラプラス変換です。
同じ一自由度粘性減衰系を例にします。運動方程式の各項はラプラス変換すると以下のようになります。なお、\(x_0=\dot{x}_0=0\) とします。
\begin{align} L[m\ddot{x}]&=ms^{2}X(s)-msx_{0}-m\dot{x}_{0}=ms^{2}X(s) \nonumber\\ L[c\dot{x}]&=c(sX(s)-x_0)=csX(s) \nonumber\\ L[kx]&=kX(s) \nonumber\\ L[f(t)]&=F(s) \nonumber \end{align}よって方程式は \((ms^{2}+cs+k)X(s)=F(s)\) となり、\(X_{st}(s)=\frac{F(s)}{k}\) と置けばラプラス変換後の伝達関数は
\begin{eqnarray} \frac{X(s)}{X_{st}(s)}=\frac{k}{ms^{2}+cs+k} \nonumber \end{eqnarray}と求まります。
SciPy の lti, bode 関数を使えば、伝達関数の分子、分母の係数を入力するだけで以下のようにボード線図が書けます。
bode 関数は Matlab だと Control toolbox を購入しないと使えませんが、Python ならフリーで同様のことができてしまいます。
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.show()
