RC回路のシミュレーション

13197 ワード
Python scipy PythonControl 電子回路 Python テキストリンク

騒音計で収録した音声から、騒音レベルを計算したい

騒音計の時間重み付け特性は、Fastが時定数125 msとなっているらしい。
これをPythonでやるために、色々調べたのでメモ。

どうも使っている騒音計では、借りてきたこの画像でいうところの、
$V_{in}$に2乗した音圧が入力され、$V_{out}$が出力されるらしい。
回路なんて全然やっていなかったから手探りで色々調べた。

ちゃんとシミュレートできているかの確認は下記サイトの出力と見比べることとした。
（素人には、何が正しいかわからんしね。ちなみに回路図、絵だと思っていたら、操作できてびっくり！）
1次RC回路の過渡解析・AC解析

色々調べていたら、 scipy.signal, control.matlab でできるみたい。
（はじめ両方混同していて、パニックに。。。）

シミュレーションしたいRC回路の周波数特性はこう。

\begin{align}
\text{周波数特性:}\quad \frac{V_{out}(s)}{V_{in}(s)}&=\frac{1}{1+sCR}\\
\end{align}

これを再現するメソッドがscipy.signal.lsim や control.matlab.lsim 。

線形時不変システムを作って、任意の信号を入力して出力を得られる。
scipyはこのサイトを参考にした。
例えば、$lti([a_1, a_2, a_3], [b_1, b_2, b_3])$ と指定すると
周波数特性が $\frac{a_1 s^2 + a_2 s + a_3}{b_1 s^2 + b_2 s + b_3}$ となり
時数は必要な数並べると増えるみたい。
RC回路だと、$lti([1], [R * C, 1])$ でよし。

でできたのが、これ。

RC_circuit_signal.py
from matplotlib import pyplot as plt
import scipy.signal as sg
import numpy as np

R    = 1000        # 1k Ohm(s/F)
C    = 0.000_000_1 # 0.1 F 
# 1 / (R * C * s + 1)
num  = [1]
den  = [R * C, 1]
t    = np.linspace(0, 0.01, 1000) 
freq = 250   
v_in = 0.5 + 0.5 * np.sign(np.sin(2 * np.pi * freq * (t - 0.001))) 

system = sg.lti(num, den) 
t_out, v_out, x = sg.lsim(system, v_in, t)

plt.style.use('grayscale')
plt.plot(t,     v_in,  label="$V_{in}$")
plt.plot(t_out, v_out, label="$V_{out}$")
plt.xlabel('Time[s]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid()
plt.show()

control.matlabは、control.tf でシステムを作る。
さっきと似たような感じだけどこんな感じになる。
lsim の戻り値の順番が違うのがいやらしい。。。

RC_circuit_control_matlab.ipynb
from matplotlib import pyplot as plt
import control.matlab as ctrl
import numpy as np

R    = 1000        # 1k Ohm(s/F)
C    = 0.000_000_1 # 0.1 F 
# 1 / (R * C * s + 1)
num  = [1]
den  = [R * C, 1]
t    = np.linspace(0, 0.01, 1000) 
freq = 250   
v_in = 0.5 + 0.5 * np.sign(np.sin(2 * np.pi * freq * (t - 0.001))) 

system = ctrl.tf(num, den)
print(system)
v_out, t_out, x = ctrl.lsim(system, v_in, t)

plt.style.use('grayscale')
plt.plot(t,     v_in,  label="$V_{in}$")
plt.plot(t_out, v_out, label="$V_{out}$")
plt.xlabel('Time[s]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid()
plt.show()

      1
------------
0.0001 s + 1

ちなみに、作ったシステムを出力すると数式が表示される。賢い！

import control.matlab as ctrl
ctrl.tf([1,2,3], [4,5,6])

$$\frac{s^2 + 2 s + 3}{4 s^2 + 5 s + 6}$$

それにしても、今回初めて知ったけど、Python Control Systems Libraryはすごいね。
数式がプログラムで書けて、しかも数式が画像できれいに出る。
しかも微分方程式解いて、シミュレーションもできる。昔と段違いだ。。。
(昔なら自分で方程式を解いた後に、差分方程式に変換してプログラムを組んだけど、今はそんな必要もないとはね。。。）

今日見たサイト。メモとして列挙。
- RC直列回路の過渡現象を sympy で解析、matplotlib で描画
- PythonControlで矩形波に対する応答を求める
- RLC直列回路の伝達関数・状態空間モデルとPythonによるシミュレーション
- Pythonで学ぶ制御工学 Part1:Pythonの基礎から制御系解析まで
- Pythonで学ぶ制御工学 Part2:伝達関数モデルを用いた制御系設計