1. 背景
PythonによるBode図表示は、scipyやcontrolといったライブラリを用いたwebサイトが多々あるが、「もう少し定量的にも把握してみたい」との動機で検討してみた。ゲイン(縦軸)を指定してその周波数(横軸)や位相(もう一つのグラフの縦軸)がいくらかを読み取り、その回路の安定性を判断することがアナログ回路設計では多い。トランジスタで構成される増幅器本体の選定や設計検証には時間を要するため、受動素子でほぼ決まる閉ループゲインを予め見積もることをpythonでも実現してみたい。
2. 目標
1.Bode図を表示し、所望のゲインに対する周波数と位相を数値で表示すること。
2.定数を変えて、波形および数値の変化が表示できること。
3.図が重ね描きできること。
を実現するスクリプトを書いてみる。
一つのライブラリに依存すると、ライブラリ更新により動作しなくなったり、更新自体がされなくなるリスクもあるため、複数のライブラリで動作確認し、特徴や可能性を評価してみた。
3. 環境
windows 10(64bit)
Jupyter Notebook 6.1.4
Python 3.8.5
control 0.9.0
scipy 1.5.2
matplotlib 3.3.2
numpy 1.19.2
4. 方針
下図の非反転増幅回路を想定し、増幅器本体は理想的とする。入力部と帰還部のインピーダンス$Z_{IN}=R_{IN}$および$Z_{F}=R_{F2}//(R_{F3}+\frac{1}{s C_{F}})$で構成される伝達関数の一定数、ここでは定数RINを配列RIN_ROWとしてfor文で3種の値99、990、9.9kΩに変更してみる。
伝達関数(導出経過省略)
1+\frac{Z_{F}}{Z_{IN}}=\frac{C_{F}(R_{F2}R_{F3}+R_{F2}R_{IN}+R_{F3}R_{IN})s+(R_{IN}R_{F2})}{C_{F}R_{IN}(R_{F2}R_{F3})s+R_{IN}}
指定するゲインは、定数GAIN＿TARGETEDと設定した。以下のスクリプトでは最大ゲインから-3dB下がったゲインを指定し、それに対応する周波数と位相を数値出力する。このゲインは出力最大振幅の半分になる周波数として、帯域幅およびゲインバンド積を評価するのによく使われる。
なお、GAIN_TARGETED = 0すればユニティゲイン周波数が一般的に求められるが、ここで採用した伝達関数では0dBで周波数軸と交差しないので、スクリプトは正常動作しない。
ゲインおよび位相曲線のデータは、それぞれ変数名mag、omegaで表される離散的な配列である。指定ゲインに対応する周波数は必ずしも配列の数値に直接対応せず、2つの配列の数値間に存在することとなり、補間法で求める。
補間法は、２点の直線補間をお手製の関数定義def myfunc_line_interpでの方法と、numpyのinterpメソッド(方法)とを利用する。いずれの方法においても、補間を行う２点の選択は、指定ゲインGAIN＿TARGETEDと、ゲインの数値列magとの大小比較することで、指定ゲインに最隣接するゲイン2点を特定する。
特定された2点のインデックスを読み取り、それらのインデックスに対応する周波数の数値、位相の数値を求める。補間で必要な２点の座標(x1,y1)および(x2,y2)は、周波数をx軸に、ゲインまたは位相をy軸とし、ゲインと位相の2種類の補間直線を決める。myfunc_line_interpでは、１次方程式の傾き(slope)と切片(offset)を求め、指定したゲインに対応する周波数(freq_myinterpまたはfreq_np_interp)を求める。
5. 結論
関数定義def myfunc_line_interpとnumpyのinterpの２方法での数値結果に誤差はあるものの、周波数で桁と小数点第１位まで一致しているので、個人的にはお手製の関数定義でも問題ないレベルと判断した。
Bode図を表示するライブラリには、以下の良し悪しがあった。matplotlibでのグラフ作成には、MATLAB形式(慣例的にpltで表記)とオブジェクト指向形式(同じくax表記)が併存しているが、正式名称がわからないので、以下、それぞれ"plt表記"および"ax表記"とする。
| ライブラリ | control | scipy + "plt表記" | scipy+"ax表記" |
|---|---|---|---|
| 長所 | グラフへの曲線の重ね描き可能。 | グラフに補助線が追加可能。 | グラフに補助線が追加可能。 |
| 短所 | グラフに補助線が追加できない？ | 非推奨の警告あり。 | グラフへの曲線の重ね描きができなかった。 |
今のところ、scipy+"plt表記"が目標には合っている。非推奨の警告があるが、for文で繰り返すインスタンス表記によるもので、現状問題なしと判断した。
6. スクリプトと評価結果
6.1. controlを用いた場合
test01.py
import numpy as np
from control.matlab import *
from matplotlib import pyplot as plt
# 伝達関数のパラメータ
CF = 16E-6
RIN = 9.9E+1
RF2 = 99E+3
RF3 = 1E+3
# 直線補間の関数定義
def myfunc_line_interp(x1,y1,x2,y2,GAIN_TARGETED):
slope = (y2-y1)/(x2-x1)
offset = -1*slope*x1+y1
return (GAIN_TARGETED-offset)/slope
# 伝達関数
num = [CF*(RF2*RF3+RF2*RIN+RF3*RIN),RIN+RF2] # 分子の係数
den = [CF*RIN*(RF2+RF3),RIN] # 分母の係数
sys = tf(num, den) # 伝達関数モデルの作成
mag, phase, omega = bode(sys, Hz=True) # ボード線図のプロット
gain_db = 20*np.log10(mag) # GainをdB単位に変換
GAIN_TARGETED = np.max(gain_db)-3 #単位dBで所望のGainを指定する。ここでは最大Gain-3dBとした。
# 指定したGain dB以上で最小のGainとその周波数の組を求める。mtt=more than target
min_gain_more_than_target = np.min(np.where(gain_db>GAIN_TARGETED, gain_db, 100))
index_min_gain_mtt = np.where(gain_db == min_gain_more_than_target)
freq_min_gain_mtt = omega[index_min_gain_mtt]/(2*np.pi)
phase_min_gain_mtt = phase[index_min_gain_mtt]/(2*np.pi)*360
# 指定したGain dB以下で最大のGainとその周波数の組を求める。ltt=less than target
max_gain_less_than_target = np.max(np.where(gain_db<GAIN_TARGETED, gain_db, -100))
index_max_gain_ltt = np.where(gain_db == max_gain_less_than_target)
freq_max_gain_ltt = omega[index_max_gain_ltt]/(2*np.pi)
phase_max_gain_ltt = phase[index_max_gain_ltt]/(2*np.pi)*360
print('指定したゲインでの周波数と位相')
# 指定したGainでの周波数を直線補間で求める。
freq_myinterp = myfunc_line_interp(freq_min_gain_mtt,min_gain_more_than_target,freq_max_gain_ltt,max_gain_less_than_target, GAIN_TARGETED)[0]
phase_myinterp = myfunc_line_interp(phase_min_gain_mtt,min_gain_more_than_target,phase_max_gain_ltt,max_gain_less_than_target, GAIN_TARGETED)[0]
print('Freqency at {:2.1f}dB (直線補間) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_myinterp))
print('Phase at {:2.1f}dB (直線補間) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_myinterp))
# numpyのinterpを用いた補間方法で求める。
freq_corner = [freq_min_gain_mtt[0],freq_max_gain_ltt[0]]
phase_corner = [phase_min_gain_mtt[0],phase_max_gain_ltt[0]]
gain_db_corner = [min_gain_more_than_target,max_gain_less_than_target]
freq_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, freq_corner)
phase_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, phase_corner)
print('Freqency at {:2.1f}dB (interp) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_np_interp))
print('Phase at {:2.1f}dB (interp) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_np_interp))
plt.show()
結果
指定したゲインでの周波数と位相
Freqency at 57.0dB (直線補間) 9.9273e-02 [Hz]
Phase at 57.0dB (直線補間) -44.3 [deg]
Freqency at 57.0dB (interp) 9.9682e-02 [Hz]
Phase at 57.0dB (interp) -44.4 [deg]
for文で繰り返し動作させてみる。上記test01.pyに対して、伝達関数の定数のうちRINを配列RIN_ROWとし、３種類の値を設定した。また、繰り返したい行の先頭にFor文を追加し、繰り返したい行をインデントさせた。以下のtest02.pyとなる。
test02.py
import numpy as np
from control.matlab import *
from matplotlib import pyplot as plt
# 伝達関数のパラメータ
CF = 16E-6
RIN_ROW = [9.9E+1, 9.9E+2, 9.9E+3]
RF2 = 99E+3
RF3 = 1E+3
# 直線補間の関数定義
def myfunc_line_interp(x1,y1,x2,y2,GAIN_TARGETED):
slope = (y2-y1)/(x2-x1)
offset = -1*slope*x1+y1
return (GAIN_TARGETED-offset)/slope
for RIN in RIN_ROW:
# 伝達関数
num = [CF*(RF2*RF3+RF2*RIN+RF3*RIN),RIN+RF2] # 分子の係数
den = [CF*RIN*(RF2+RF3),RIN] # 分母の係数
sys = tf(num, den) # 伝達関数モデルの作成
mag, phase, omega = bode(sys, Hz=True) # ボード線図のプロット
gain_db = 20*np.log10(mag) # GainをdB単位に変換
GAIN_TARGETED = np.max(gain_db)-3 #単位dBで所望のGainを指定する。ここでは最大Gain-3dBとした。
# 指定したGain dB以上で最小のGainとその周波数の組を求める。mtt=more than target
min_gain_more_than_target = np.min(np.where(gain_db>GAIN_TARGETED, gain_db, 100))
index_min_gain_mtt = np.where(gain_db == min_gain_more_than_target)
freq_min_gain_mtt = omega[index_min_gain_mtt]/(2*np.pi)
phase_min_gain_mtt = phase[index_min_gain_mtt]/(2*np.pi)*360
# 指定したGain dB以下で最大のGainとその周波数の組を求める。ltt=less than target
max_gain_less_than_target = np.max(np.where(gain_db<GAIN_TARGETED, gain_db, -100))
index_max_gain_ltt = np.where(gain_db == max_gain_less_than_target)
freq_max_gain_ltt = omega[index_max_gain_ltt]/(2*np.pi)
phase_max_gain_ltt = phase[index_max_gain_ltt]/(2*np.pi)*360
print('指定したゲインでの周波数と位相')
# 指定したGainでの周波数を直線補間で求める。
freq_myinterp = myfunc_line_interp(freq_min_gain_mtt,min_gain_more_than_target,freq_max_gain_ltt,max_gain_less_than_target, GAIN_TARGETED)[0]
phase_myinterp = myfunc_line_interp(phase_min_gain_mtt,min_gain_more_than_target,phase_max_gain_ltt,max_gain_less_than_target, GAIN_TARGETED)[0]
print('Freqency at {:2.1f}dB (直線補間) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_myinterp))
print('Phase at {:2.1f}dB (直線補間) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_myinterp))
# numpyのinterpを用いた補間方法で求める。
freq_corner = [freq_min_gain_mtt[0],freq_max_gain_ltt[0]]
phase_corner = [phase_min_gain_mtt[0],phase_max_gain_ltt[0]]
gain_db_corner = [min_gain_more_than_target,max_gain_less_than_target]
freq_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, freq_corner)
phase_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, phase_corner)
print('Freqency at {:2.1f}dB (interp) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_np_interp))
print('Phase at {:2.1f}dB (interp) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_np_interp))
plt.show()
結果
指定したゲインでの周波数と位相
Freqency at 57.0dB (直線補間) 9.9273e-02 [Hz]
Phase at 57.0dB (直線補間) -44.3 [deg]
Freqency at 57.0dB (interp) 9.9682e-02 [Hz]
Phase at 57.0dB (interp) -44.4 [deg]
指定したゲインでの周波数と位相
Freqency at 37.1dB (直線補間) 9.9300e-02 [Hz]
Phase at 37.1dB (直線補間) -43.8 [deg]
Freqency at 37.1dB (interp) 9.9682e-02 [Hz]
Phase at 37.1dB (interp) -43.9 [deg]
指定したゲインでの周波数と位相
Freqency at 17.8dB (直線補間) 1.0027e-01 [Hz]
Phase at 17.8dB (直線補間) -39.5 [deg]
Freqency at 17.8dB (interp) 1.0084e-01 [Hz]
Phase at 17.8dB (interp) -39.6 [deg]
グラフに曲線を重ね描きできたが、補助線を引くのはわからず、ここで深堀りの検討はとどめた。
6.2. scipy+ "plt表記"
以下、for文を用いたスクリプトのみ記載する。まずは"plt表記"のスクリプト。
test03.py
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
# 伝達関数のパラメータ
CF = 16E-6
RIN_ROW = [9.9E+1, 9.9E+2, 9.9E+3]
RF2 = 99E+3
RF3 = 1E+3
# def関数での直線補間
def myfunc_line_interp(x1,y1,x2,y2,GAIN_TARGETED):
slope = (y2-y1)/(x2-x1)
offset = -1*slope*x1+y1
return (GAIN_TARGETED-offset)/slope
for RIN in RIN_ROW:
# 伝達関数
num = [CF*(RF2*RF3+RF2*RIN+RF3*RIN),RIN+RF2] # 分子の係数
den = [CF*RIN*(RF2+RF3),RIN] # 分母の係数
sys = signal.lti(num, den)
# ボード線図の計算
w, mag, phase = signal.bode(sys)
GAIN_TARGETED = np.max(mag)-3 #単位dBで所望のGainを指定する。ここでは最大Gain-3dBとした。
# 指定したGain dB以上で最小のGainとその周波数の組を求める。mtt=more than target
min_gain_more_than_target = np.min(np.where(mag>GAIN_TARGETED, mag, 100))
index_min_gain_mtt = np.where(mag == min_gain_more_than_target)
freq_min_gain_mtt = w[index_min_gain_mtt]/(2*np.pi)
phase_min_gain_mtt = phase[index_min_gain_mtt]
# 指定したGain dB以下で最大のGainとその周波数の組を求める。ltt=less than target
max_gain_less_than_target = np.max(np.where(mag<GAIN_TARGETED, mag, -100))
index_max_gain_ltt = np.where(mag == max_gain_less_than_target)
freq_max_gain_ltt = w[index_max_gain_ltt]/(2*np.pi)
phase_max_gain_ltt = phase[index_max_gain_ltt]
print('指定したゲインでの周波数と位相')
# def関数での直線補間
freq_myinterp = myfunc_line_interp(freq_min_gain_mtt,min_gain_more_than_target,freq_max_gain_ltt,max_gain_less_than_target,GAIN_TARGETED)[0]
phase_myinterp = myfunc_line_interp(phase_min_gain_mtt,min_gain_more_than_target,phase_max_gain_ltt,max_gain_less_than_target,GAIN_TARGETED)[0]
print('Freqency at {:2.1f}dB (直線補間) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_myinterp))
print('Phase at {:2.1f}deg (直線補間) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_myinterp))
# numpyのinterp
freq_corner = [freq_min_gain_mtt[0],freq_max_gain_ltt[0]]
phase_corner = [phase_min_gain_mtt[0],phase_max_gain_ltt[0]]
gain_db_corner = [min_gain_more_than_target,max_gain_less_than_target]
freq_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, freq_corner)
phase_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, phase_corner)
print('Freqency at {:2.1f}dB (interp) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_np_interp))
print('Phase at {:2.1f}dB (interp) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_np_interp))
# ゲイン線図の描画
plt.subplot(2, 1, 1)
plt.semilogx(w/(2*np.pi), mag)
plt.hlines(xmin=1E-3, xmax=freq_myinterp, y=GAIN_TARGETED, linestyles='dotted')
plt.vlines(x=freq_myinterp, ymin=0, ymax=GAIN_TARGETED, linestyles='dotted')
plt.ylabel("Gain[dB]")
plt.grid(which='major',color='black',linestyle='-')
plt.grid(which='minor',color='gray',linestyle='--')
# 位相線図の描画
plt.subplot(2, 1, 2)
plt.semilogx(w/(2*np.pi), phase)
plt.hlines(xmin=1E-3, xmax=freq_np_interp, y=phase_np_interp, linestyles='dotted')
plt.vlines(x=freq_np_interp, ymin=phase_np_interp, ymax=0, linestyles='dotted')
plt.xlabel("Frequency[Hz]")
plt.ylabel("Phase[deg]")
plt.grid(which='major',color='black',linestyle='-')
plt.grid(which='minor',color='gray',linestyle='--')
plt.show()
結果
指定したゲインでの周波数と位相
Freqency at 57.0dB (直線補間) 9.9273e-02 [Hz]
Phase at 57.0deg (直線補間) -44.3 [deg]
Freqency at 57.0dB (interp) 1.0471e-01 [Hz]
Phase at 57.0dB (interp) -45.8 [deg]
指定したゲインでの周波数と位相
Freqency at 37.1dB (直線補間) 9.9301e-02 [Hz]
Phase at 37.1deg (直線補間) -43.8 [deg]
Freqency at 37.1dB (interp) 1.0471e-01 [Hz]
Phase at 37.1dB (interp) -45.3 [deg]
指定したゲインでの周波数と位相
Freqency at 17.8dB (直線補間) 1.0027e-01 [Hz]
Phase at 17.8deg (直線補間) -39.4 [deg]
Freqency at 17.8dB (interp) 1.0471e-01 [Hz]
Phase at 17.8dB (interp) -40.5 [deg]
定性的にも、定量的にも納得いく出力が得られた。なお、以下の非推奨の警告(Deprecation Warning)があった。本目的、目標では軸やその名称をfor文で繰り返すごとに変更しないため、問題なしと判断した。
<ipython-input-6-1da77193a555>:55: MatplotlibDeprecationWarning:
Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.
In a future version, a new instance will always be created and returned.
Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
plt.subplot(2, 1, 1)
.....
6.3. scipy+ "ax表記"
test04.py
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
# 伝達関数のパラメータ
CF = 16E-6
RIN_ROW = [9.9E+1, 9.9E+2, 9.9E+3]
RF2 = 99E+3
RF3 = 1E+3
# def関数での直線補間
def myfunc_line_interp(x1,y1,x2,y2,GAIN_TARGETED):
slope = (y2-y1)/(x2-x1)
offset = -1*slope*x1+y1
return (GAIN_TARGETED-offset)/slope
for RIN in RIN_ROW:
# 伝達関数
num = [CF*(RF2*RF3+RF2*RIN+RF3*RIN),RIN+RF2] # 分子の係数
den = [CF*RIN*(RF2+RF3),RIN] # 分母の係数
sys = signal.lti(num, den)
# ボード線図の計算
w, mag, phase = signal.bode(sys)
GAIN_TARGETED = np.max(mag)-3 #単位dBで所望のGainを指定する。ここでは最大Gain-3dBとした。
# 指定したGain dB以上で最小のGainとその周波数の組を求める。mtt=more than target
min_gain_more_than_target = np.min(np.where(mag>GAIN_TARGETED, mag, 100))
index_min_gain_mtt = np.where(mag == min_gain_more_than_target)
freq_min_gain_mtt = w[index_min_gain_mtt]/(2*np.pi)
phase_min_gain_mtt = phase[index_min_gain_mtt]
# 指定したGain dB以下で最大のGainとその周波数の組を求める。ltt=less than target
max_gain_less_than_target = np.max(np.where(mag<GAIN_TARGETED, mag, -100))
index_max_gain_ltt = np.where(mag == max_gain_less_than_target)
freq_max_gain_ltt = w[index_max_gain_ltt]/(2*np.pi)
phase_max_gain_ltt = phase[index_max_gain_ltt]
print('指定したゲインでの周波数と位相')
# def関数での直線補間
freq_myinterp = myfunc_line_interp(freq_min_gain_mtt,min_gain_more_than_target,freq_max_gain_ltt,max_gain_less_than_target,GAIN_TARGETED)[0]
phase_myinterp = myfunc_line_interp(phase_min_gain_mtt,min_gain_more_than_target,phase_max_gain_ltt,max_gain_less_than_target,GAIN_TARGETED)[0]
print('Freqency at {:2.1f}dB (直線補間) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_myinterp))
print('Phase at {:2.1f}deg (直線補間) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_myinterp))
# numpyのinterp
freq_corner = [freq_min_gain_mtt[0],freq_max_gain_ltt[0]]
phase_corner = [phase_min_gain_mtt[0],phase_max_gain_ltt[0]]
gain_db_corner = [min_gain_more_than_target,max_gain_less_than_target]
freq_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, freq_corner)
phase_np_interp = np.interp(GAIN_TARGETED, gain_db_corner, phase_corner)
print('Freqency at {:2.1f}dB (interp) {:2.4e} [Hz]'.format(GAIN_TARGETED, freq_np_interp))
print('Phase at {:2.1f}dB (interp) {:3.1f} [deg]'.format(GAIN_TARGETED, phase_np_interp))
fig = plt.figure()
# ゲイン線図の描画
ax1 = fig.add_subplot(2, 1, 1)
ax1.plot(w/(2*np.pi), mag)
ax1.set_xscale("log")
ax1.hlines(xmin=1E-3, xmax=freq_myinterp, y=GAIN_TARGETED, linestyles='dotted')
ax1.vlines(x=freq_myinterp, ymin=0, ymax=GAIN_TARGETED, linestyles='dotted')
ax1.set_ylabel("Gain[dB]")
ax1.grid(which='major',color='black',linestyle='-')
ax1.grid(which='minor',color='gray',linestyle='--')
# 位相線図の描画
ax2 = fig.add_subplot(2, 1, 2)
ax2.plot(w/(2*np.pi), phase)
ax2.set_xscale("log")
ax2.hlines(xmin=1E-3, xmax=freq_np_interp, y=phase_np_interp, linestyles='dotted')
ax2.vlines(x=freq_np_interp, ymin=phase_np_interp, ymax=0, linestyles='dotted')
ax2.set_xlabel("Frequency[Hz]")
ax2.set_ylabel("Phase[deg]")
ax2.grid(which='major',color='black',linestyle='-')
ax2.grid(which='minor',color='gray',linestyle='--')
plt.show()
結果
数値出力は"plt表記"と同じなので省略する。警告はなかったものの、現時点ではグラフは重ね描きができなかった。当初の目標は別途達成できたので深堀の検討はしない。
7.参考サイト
利用したライブラリ
[1]https://python-control.readthedocs.io/en/latest/index.html
[2]https://docs.scipy.org/doc/scipy/reference/index.html
利用したBode図表示メソッド
[3]https://python-control.readthedocs.io/en/latest/generated/control.matlab.bode.html?highlight=bode
[4]https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lti.html
matplotlibのグラフ表示形式
[5]https://matplotlib.org/2.0.2/faq/usage_faq.html#coding-styles
8.考察（利用範囲の限界や今後の課題）
伝達関数は他にはやっていないが、ローパスフィルタ特性を持ち、一つのゲインに対して一つの周波数、位相が対応する伝達特性なら問題ないと考えています。
ただ、補間で配列magから特定した２点と指定したゲインGAIN_TARGETEDとが、ドンピシャで等しくなるケースは条件に含んでいないので、正常動作しないかもしれません。
9.感想
今までexcelでチマチマとやっていたが、汎用性やポータビリティが低く、お盆休みを利用してpython化してみました。pythonは、目的に合うライブラリ発掘と、自力でのスクリプト作成の両輪でアプローチすることになるが、経験や情報が少ない初心者には苦労しますね。もっとスマートな方法があるかもしれません。でも、pythonは初心者にも入りやすくありがたいです。
「定量的にも」との個人的目標で動作確認をしましたが、「周波数帯域が低すぎでは？」といった疑問がその筋の専門の方には思い浮かぶかもしれません。本資料はスクリプトの動作確認のためであり、数値はテキトーですので、突っ込みはご勘弁を。。。