#騒音計で収録した音声から、騒音レベルを計算したい
騒音計の時間重み付け特性は、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])$ でよし。
でできたのが、これ。
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 の戻り値の順番が違うのがいやらしい。。。
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はすごいね。
数式がプログラムで書けて、しかも数式が画像できれいに出る。
しかも微分方程式解いて、シミュレーションもできる。昔と段違いだ。。。
(昔なら自分で方程式を解いた後に、差分方程式に変換してプログラムを組んだけど、今はそんな必要もないとはね。。。）
今日見たサイト。メモとして列挙。