# 贝塞尔滤波器

## 传递函数

${\displaystyle H(s)={\frac {\theta _{n}(0)}{\theta _{n}(s/\omega _{0})}}\,}$

## 简单例子

${\displaystyle H(s)={\frac {15}{s^{3}+6s^{2}+15s+15}}.}$

gain值为

${\displaystyle G(\omega )=|H(j\omega )|={\frac {15}{\sqrt {\omega ^{6}+6\omega ^{4}+45\omega ^{2}+225}}}.\,}$

${\displaystyle \phi (\omega )=-\arg(H(j\omega ))=-\arctan \left({\frac {15\omega -\omega ^{3}}{15-6\omega ^{2}}}\right).\,}$

${\displaystyle D(\omega )=-{\frac {d\phi }{d\omega }}={\frac {6\omega ^{4}+45\omega ^{2}+225}{\omega ^{6}+6\omega ^{4}+45\omega ^{2}+225}}.\,}$

${\displaystyle D(\omega )=1-{\frac {\omega ^{6}}{225}}+{\frac {\omega ^{8}}{1125}}+\cdots .}$

## Pseudo code

n ← 0 TO 1 BY 0.01
for i (0 to 1):
IF i ←← 0 THEN
pos ← 1
N ← 5
ENDIF
IF i ←← 1 THEN
pos ← 3
N ← 10
ENDIF
z, p, k ← besselap(N)
b, a ← zpk2tf(z, p, k)
h, w ← freqs(b, a, n)
magh2 ← abs(h) ** 2
phah ← unwrap(angle(h))
plt.subplot(2, 2, pos)
plt.xlim(0, 2)
plt.xlabel("w/wc")
plt.ylabel("Bessel H(jw)^2")
plt.title("N =" + INT_TO_STRING(N))
plt.grid()
plt.plot(w, magh2)
plt.subplot(2, 2, pos + 1)
plt.xlim(0, 2)
plt.xlabel("w/wc")
plt.ylabel("Bessel Ph(jw)^2")
plt.title("N =" + INT_TO_STRING(N))
plt.grid()
plt.plot(w, phah)
ENDFOR
plt.show()