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

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

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

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

![RC回路](https://zawa2.com/ZZsim/images/doc_rc.gif) 

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

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

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

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

これを再現するメソッドが[scipy.signal.lsim](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lsim.html) や [control.matlab.lsim](https://python-control.readthedocs.io/en/0.8.3/generated/control.matlab.lsim.html) 。

線形時不変システムを作って、任意の信号を入力して出力を得られる。
scipyは[このサイト](https://gsmcustomeffects.hatenablog.com/entry/2019/01/13/112150)を参考にした。
例えば、$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])$ でよし。

でできたのが、これ。

```python: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()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455423/36c9de1d-8382-f615-ab74-45c0a2191cb6.png)

control.matlabは、[control.tf](https://python-control.readthedocs.io/en/0.8.3/generated/control.tf.html#control.tf) でシステムを作る。
さっきと似たような感じだけどこんな感じになる。
lsim の戻り値の順番が違うのがいやらしい。。。

```python: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
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455423/23e8820b-89a7-eed2-7e1d-f8094b8eaea5.png)


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

```python
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 で描画](https://qiita.com/code0327/items/60cdff88c7db1750a8cd)
 - [PythonControlで矩形波に対する応答を求める](https://qiita.com/nnn112358/items/dc6be19eae4e0a16ae61)
 - [RLC直列回路の伝達関数・状態空間モデルとPythonによるシミュレーション](https://qiita.com/code0327/items/423b0f0380e8c64f3580)
 - [Pythonで学ぶ制御工学 Part1:Pythonの基礎から制御系解析まで](https://y373.sakura.ne.jp/data/sice_visit_study2020/part1.pdf)
 - [Pythonで学ぶ制御工学 Part2:伝達関数モデルを用いた制御系設計](https://y373.sakura.ne.jp/data/sice_visit_study2020/part2.pdf)