Pythonでフィルタ(6)

\[ \frac{x_n+x_{n-1}}{2}=RC\frac{dy_n}{dt}+y_n \tag{1} \]

\[ s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}} \tag{2} \]

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

ここで、 $y_n=y$, $y_{n-1}=z^{-1}y$, $x_n=x$, $x_{n-1}=z^{-1}x$ という書き方にすると、

\[ y=b_0x+b_1 z^{-1} x-a_1 z^{-1} y \tag{4} \]

\[\begin{eqnarray} H\left(z \right)&=&\frac{y}{x} \tag{5}\\ H\left(z \right)&=&\frac{b_0+b_1 z^{-1}}{1+a_1 z^{-1}} \tag{6} \end{eqnarray}\]

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

こいつで、$s=j\omega$ってするらしい。なんでかわからんが決まりらしい。すると、

\[ G\left(s \right)=\frac{1}{1+CRs} \tag{8} \]

こいつに式(2)を適用すると、あらふしぎ、$G\left(s \right)$ではなく、$H\left(z \right)$というものに変わるという魔法(理論は無視するので掘りません)。

\[ H\left(z \right)=\frac{1}{1+CR\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}} \tag{9} \]

\[ H\left(z \right)= \frac { \frac{1}{1+\frac{2CR}{T}}+\frac{1}{1+\frac{2CR}{T}}z^{-1} } { 1+\frac{1-\frac{2CR}{T}}{1+\frac{2CR}{T}}z^{-1} } \tag{10} \]

式(6)と式(10)は$y_n$, $y_{n-1}$, $x_n$, $x_{n-1}$の恒等式でなければならない、すなわち$z^{-1}$の恒等式でなければならないので、

\[ b_0=\frac{1}{1+\frac{2CR}{T}} \tag{11} \]

\[ b_1=\frac{1}{1+\frac{2CR}{T}} \tag{12} \]

\[ a_1=\frac{1-\frac{2CR}{T}}{1+\frac{2CR}{T}} \tag{11} \]

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)
two_CR_div_T=2*CR_value/t_sampling_simulation

b_0=1/(1+two_CR_div_T)
b_1=1/(1+two_CR_div_T)
a_0=1
a_1=(1-two_CR_div_T)/(1+two_CR_div_T)

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()