python · 22日 4月 2023

Pythonでフィルタ(5)

フィルタの構成方法と適用方法について記録を残す。理論はさておき実際にフィルタを適用した波形を見ることを目的とする。ここまでMATLAB-style IIRフィルタをやってきたけど、アナログ回路での時間ドメイン計算ができればいいじゃんって考える。

RC LPFです。

Cを流れる電流を$i_c$とすると、

\[ i_c=C\frac{dv_o\left(t \right)}{dt} \tag{1}\]

$v_o$端子から流れ出る電流をゼロとすると、

\[ v_i\left(t \right)=RC\frac{dv_o\left(t \right)}{dt}+v_o\left(t \right) \tag{2}\]

サンプリング周期をTとして、$t_0=\left( n-1 \right)T, t_1=nT$とすると、

\[ \begin{eqnarray} \frac{dv_o\left(t \right)}{dt}&=&\frac{v_o\left(t_1 \right)-v_o\left(t_0 \right)}{t_1-t_0} \\ &=&\frac{v_o\left(nT \right)-v_o\left( \left(n-1 \right)T \right)}{T} \end{eqnarray} \tag{3}\]

(3)を(2)に代入して整理すると、

\[ v_o\left(nT \right)=\frac{\frac{RC}{T}}{\frac{RC}{T}+1}v_o\left(\left(n-1 \right)T \right)+\frac{1}{\frac{RC}{T}+1}v_i\left(nT \right) \tag{4} \]

いっこ前の出力と現時点の入力から現時点の出力を計算できるようになりました。

次は、周波数特性を計算します。

抵抗Rのインピーダンスを$Z_R$、コンデンサCのインピーダンスを$Z_C$とすると、

\[ v_o\left(t \right)=\frac{Z_C}{Z_R+Z_C}v_i\left(t \right) \tag{5} \]

$Z_R=R, Z_C=\frac{1}{j\omega C}$なので、ゲイン$G\left(j\omega \right)$は、

\[ G\left(j\omega \right)=\frac{v_o\left(t \right)}{v_i\left(t \right)}=\frac{1}{1+j\omega CR} \tag{6} \]

となる。で、カットオフ周波数を計算してみる。カットオフ周波数は、

電圧が$\frac{1}{\sqrt{2}}$となる周波数なので、

\[ |G\left(j\omega_c \right)|=\frac{1}{\sqrt{2}} \tag{7} \]

\[ \frac{1}{\sqrt{1+\left(\omega_c CR \right)^2}}=\frac{1}{\sqrt{2}} \tag{8} \]

\[ \omega=\frac{1}{CR} \tag{9} \]

\[ f_c=\frac{1}{2 \pi CR} \tag{10} \]

まぁよく見るやつです。

ちょいと、計算してみます。目指す特性はこれまでと同じです。

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

# よく使う変数
pi=np.pi
deg2rad=pi/180.0
rad2deg=180.0/pi
twopi=2*pi
        
f_carrier=500e3

#フィルタ生成
q_value=10
zeta=1/(2*q_value)
fcut=zeta*f_carrier
print(fcut)

omega_c=twopi*fcut
R_value=100
C_value=1/(omega_c*R_value)
print(R_value)
print(C_value)

f_start=1
f_end=50e3
f_step=1
freqz=np.arange(f_start,f_end,f_step)
omegaz=freqz*twopi

G_jw=1/(1+1j*omegaz*C_value*R_value)
G_jw_dB=20*np.log10(np.abs(G_jw))
G_jw_ang=np.angle(G_jw)*rad2deg

#プロット
fig = plt.figure()
ax1=plt.subplot(2,1,1)
ax1.grid(True)
ax1.plot(freqz,G_jw_dB, '-')
ax1.set_xlim(0,50e3)
ax1.set_ylim(-10,0)
ax1.set_xlabel('frequency [Hz]')
ax1.set_ylabel('Amplitude [dB]')
ax1_2=plt.subplot(2,1,2)
#ax1_2=ax1.twinx()
ax1_2.grid(True)
ax1_2.plot(freqz,G_jw_ang, 'r-')
ax1_2.set_xlim(0,50e3)
ax1_2.set_ylim(-180,180)
ax1_2.set_yticks(np.arange(-180,181,step=90))
ax1_2.set_xlabel('frequency [Hz]')
ax1_2.set_ylabel('Phase [deg.]')

fig.tight_layout()
plt.show()

カットオフ周波数は25kHzになっていて、うまくできているようです。

ところでRを100Ωってしたけど、これ決めるのってなかなか難しい。実際はフィルタの入力側の出力インピーダンスやフィルタの出力側の負荷の入力インピーダンスで決めるものらしい。がしかし、今回は実際の回路は組まずに計算で波形を出すことが目的なので、適当なんです。そして、今回は具体的なCとRの値は不要で、CR(掛け算の結果)だけあればよかったんです。Rを100Ωってする必要すらなかった。

次は時間ドメインでの計算にチャレンジしたいが、その前に、これまでのコードの再利用のために、前回のブロック線図をいじってみる。

この絵から、

\[ \begin{eqnarray} p_n&=&x_n-a_1p_{n-1} \\ y_n&=&b_0p_n-b_1p_{n-1} \end{eqnarray} \]

こんな感じでフィルタを当てました。

で、この式をこねくり回すと、

\[ y_n=b_0x_n+b_1x_{n-1}-a_1y_{n-1} \tag{11} \]

たぶん(、あちこちのサイトがブロック線図の変形をさも当たり前のようにやっているので)、まずブロック線図をこねくり回してから式を導出するのが正しい教育を受けた人の手順なんだろうけど、我流なので。この式からブロック線図を書くとこうなる。

（ちなみに、最近いろんなサイトを見てるとaとbの関係がわけわかんなくなってきた、あっちこっちのサイトで必ずしも一致していないような、誤記なのか、、日曜技術者はscipyでの定義に合わせていくことにしました。）

で、この絵または式(11)と式(4)を眺めると、

\[ a_1=-\frac{\frac{RC}{T}}{\frac{RC}{T}+1} \tag{12} \] \[ b_0=\frac{1}{\frac{RC}{T}+1} \tag{13} \] \[ b_1=0 \tag{14} \]

scipy.signal.butterとは違う値になっていると思う。で、それでいいのかどうか。

こんな感じで比較する。

(ちなみにa0=1は常識らしい)

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

# よく使う変数
pi=np.pi
deg2rad=pi/180.0
rad2deg=180.0/pi
twopi=2*pi

f_sampling_simulation=5e6
f_carrier=500e3
t_sampling_simulation=1/f_sampling_simulation

q_value=10
zeta=1/(2*q_value)
fcut=0+zeta*f_carrier
print(fcut)

#フィルタ生成(butter)
wcut=fcut/(f_sampling_simulation/2)
b,a=signal.butter(1,wcut,'lowpass')
print(b)
print(a)
w,h=signal.freqz(b,a,worN=2**20)
freqs_but=w*f_sampling_simulation/twopi
h_abs_but=np.abs(h)
amps_db_but=20*np.log10(h_abs_but)
angles_but=np.angle(h)*rad2deg

#フィルタ生成(CR LPF)
CR_value=1/(twopi*fcut)
CR_div_T=CR_value/t_sampling_simulation

b_0=1/(CR_div_T+1)
b_1=0
a_0=1
a_1=-CR_div_T/(CR_div_T+1)

b=np.array([b_0,b_1])
a=np.array([a_0,a_1])
print(b)
print(a)
w,h=signal.freqz(b,a,worN=2**20)
freqs_cr=w*f_sampling_simulation/twopi
h_abs_cr=np.abs(h)
amps_db_cr=20*np.log10(h_abs_cr)
angles_cr=np.angle(h)*rad2deg

#プロット
fig = plt.figure()
ax1=plt.subplot(2,2,1)
ax1.grid(True)
ax1.plot(freqs_but,amps_db_but, '-')
ax1.plot(freqs_cr,amps_db_cr, '-')
ax1.set_xlim(0,50e3)
ax1.set_ylim(-10,0)
ax1.set_xlabel('frequency [Hz]')
ax1.set_ylabel('Amplitude [dB]')

ax2=plt.subplot(2,2,2)
ax2.grid(True)
ax2.plot(freqs_but,amps_db_but, '-')
ax2.plot(freqs_cr,amps_db_cr, '-')
ax2.set_xlim(1e2,1e7)
ax2.set_xscale('log')
ax2.set_ylim(-50,0)
ax2.set_xlabel('frequency [Hz]')
ax2.set_ylabel('Amplitude [dB]')

ax3=plt.subplot(2,2,3)
ax3.grid(True)
ax3.plot(freqs_but,angles_but, '-')
ax3.plot(freqs_cr,angles_cr, '-')
ax3.set_xlim(0,50e3)
ax3.set_ylim(-90,90)
ax3.set_xlabel('frequency [Hz]')
ax3.set_ylabel('angle [deg.]')

ax4=plt.subplot(2,2,4)
ax4.grid(True)
ax4.plot(freqs_but,angles_but, '-')
ax4.plot(freqs_cr,angles_cr, '-')
ax4.set_xlim(1e2,1e7)
ax4.set_xscale('log')
ax4.set_ylim(-90,90)
ax4.set_xlabel('frequency [Hz]')
ax4.set_ylabel('angle [deg.]')


fig.tight_layout()
plt.show()

コンソール出力はこうなる。

25000.0

[0.01546629 0.01546629]

[ 1. -0.96906742]

[0.03045903 0. ]

[ 1. -0.96954097]

aはほぼ同じだけど、bが異なる。

で、周波数特性はこうなる。

カットオフ周波数付近ではほぼ同じ特性。高周波側でかなり違いが出てくる。これが時間ドメイン波形にどう影響するのか、、、

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

# よく使う変数
pi=np.pi
deg2rad=pi/180.0
rad2deg=180.0/pi
twopi=2*pi
 
# シンボル列を引き延ばす
def extend_symbol_to_simulation_sampling_rate(symbols,t_symbol,t_simulation_sampling):
    # シンボル数
    size_of_symbols=symbols.size
    # 最後のシンボルが終わる時間
    end_time=size_of_symbols*t_symbol
    # 計算タイミング
    t=np.arange(0,end_time,t_simulation_sampling)
    # 各シンボルの開始（終了）時間
    t_periods_of_symbols=np.append(0,np.arange(1,size_of_symbols,1)*t_symbol)
    # 結果データ領域確保
    ex_symbols=np.zeros_like(t)
    # 先頭データは入れておく
    ex_symbols[0]=symbols[0]
    # 各シンボルについて、その期間のインデックスを取得し、データを入れ込む
    for i in range(0,size_of_symbols):
        ex_symbols[np.where(t>t_periods_of_symbols[i])]=symbols[i]
    return ex_symbols,end_time
        
f_sampling_simulation=5e6
f_carrier=500e3

#generate baseband waveform
f_symbol=10e3                 # シンボルレート in Hz
t_sampling_simulation=1/f_sampling_simulation
symbols_in_base=np.array([1,0,1,0,1,0,1,0,1,0,1,0])
baseband_t,t_end=extend_symbol_to_simulation_sampling_rate(symbols_in_base,1/f_symbol,t_sampling_simulation)
t=np.arange(0,t_end,t_sampling_simulation)
number_of_samples=t.size

#フィルタ生成(butter)
q_value=10
zeta=1/(2*q_value)
fcut=0+zeta*f_carrier
print(fcut)
wcut=fcut/(f_sampling_simulation/2)
b,a=signal.butter(1,wcut,'lowpass')
print(b)
print(a)
#w,h=signal.freqz(b,a,worN=2**20)
#freqs=w*f_sampling_simulation/twopi
#h_abs=np.abs(h)
#amps_db=20*np.log10(h_abs)
#angles=np.angle(h)*rad2deg

#フィルタ適用
baseband_filtered_t=signal.lfilter(b,a,baseband_t)
#保存
np.savez('o_t_baseband_filtered_t',t,baseband_filtered_t)

#フィルタ生成(CR LPF)
CR_value=1/(twopi*fcut)
CR_div_T=CR_value/t_sampling_simulation

b_0=1/(CR_div_T+1)
b_1=0
a_0=1
a_1=-CR_div_T/(CR_div_T+1)

b=np.array([b_0,b_1])
a=np.array([a_0,a_1])
print(b)
print(a)
w,h=signal.freqz(b,a,worN=2**20)
freqs=w*f_sampling_simulation/twopi
h_abs=np.abs(h)
amps_db=20*np.log10(h_abs)
angles=np.angle(h)*rad2deg
#フィルタ適用(自力2)
baseband_filtered_t=np.zeros_like(baseband_t)
prevx=0
prevy=0
for i in range(t.size):
    presx=baseband_t[i]
    presy=b[0]*presx+b[1]*prevx-a[1]*prevy
    baseband_filtered_t[i]=presy
    prevx=presx
    prevy=presy

#butterを適用した波形をファイルからロード
npz = np.load('o_t_baseband_filtered_t.npz')
t_from_file=npz['arr_0']
amp_t_from_file=npz['arr_1']

#プロット
fig = plt.figure()
ax1=plt.subplot(3,1,1)
ax1.grid(True)
ax1.plot(freqs,amps_db, '-')
ax1.set_xlim(0,50e3)
ax1.set_ylim(-10,0)
ax1.set_xlabel('frequency [Hz]')
ax1.set_ylabel('Amplitude [dB]')
ax1_2=ax1.twinx()
ax1_2.grid(True)
ax1_2.plot(freqs,angles, 'r-')
ax1_2.set_xlim(0,50e3)
ax1_2.set_ylim(-90,90)
ax1_2.set_yticks(np.arange(-90,91,step=45))
ax1.set_xlabel('frequency [Hz]')
ax1_2.set_ylabel('Phase [deg.]')

ax2=plt.subplot(3,1,2)
ax2.plot(t,baseband_t)

ax3=plt.subplot(3,1,3)
ax3.plot(t,baseband_filtered_t)
ax3.plot(t_from_file,amp_t_from_file)

fig.tight_layout()
plt.show()

相変わらず、説明なしのいきなりのコード。２つの計算結果を混同してしまわないように(絶対に混同していないことを担保するために)butterでの計算結果は一旦ファイルに放り込んでいます。

ほぼ一緒。拡大してみる。

ほぼ一緒。

ということで、scipy.signal.butterとアナログCR LPFはほぼ同じってことが確認できました。フィルタ係数のbを見比べると、butterはアナログ回路のいっこ前の入力を使っているのではなく0.5前と0.5後を使ったものと一致させることを狙っているんじゃなかろうかと思う。確認したいけど、しんどいのでいつか気が向いたらやる。CR LPFでサンプリング周波数付近で減衰量が悪化するのも原理的にはあってない気がする。

butterの計算を掘り進めるのがしんどい(歳とると腰が引けてくるんです)ので、CR LPFを使っていこうかとちょっと思ったが、前述のような理由からやっぱりbutterかなー。

さて、理論はさておき、とか言いながら、今回、式書いてるやん、、、って言わんでください。理論ではない当たり前のことをメモとして書き残しただけなんです(理由のない上から目線の言い訳)；

ああーあれもこれもせんといかんー。こんなことやっている場合じゃないー｡ﾟ(ﾟ´Д｀ﾟ)ﾟ｡

カテゴリ：非現実

コメントをお書きください

コメント: 0