模拟滤波器的基础知识和设计

原创

已于 2023-10-07 22:42:47 修改 · 3.3k 阅读

25 ·

CC 4.0 BY-SA版权

文章标签：

#信号处理 #python #scipy

于 2023-10-07 22:42:18 首次发布

该文章已生成可运行项目，

信号处理工作中滤波器的应用是非常广泛的，可以分成模拟滤波器和数字滤波器两种，数字滤波器主要包括两种，IIR和FIR，这两种滤波器后面统一说，今天先来说一说模拟滤波器（主要是我先用Python实现了Matlab书里面模拟滤波器的一些内容）。

首先，什么是滤波器，什么又是模拟滤波器？

滤波器：具有频率选择作用的电路或者运算处理系统，具有滤除噪声、分离不同信号的功能，今天主要写的是，1、巴特沃斯滤波器、2、切比雪夫滤波器，3、椭圆滤波器，4、低通到低通的频带转换

模拟滤波器：更具一组设计规范来设计模拟系统函数 $H_a(s)$ ,使其逼近某个理想滤波器的特性。

各种模拟滤波器的设计过程都是先设计出低通滤波器，再通过频率变换将低通滤波器转换成其他类型模拟滤波器。

我们考虑因果系统：

$H_a(j\Omega)=\int_{0}^{\inf}h_a(t)e^{-j\omega t}dt$

其中， $h_a(s)$ 是系统的单位脉冲响应，是实函数，那么就有：

$H_a(j\Omega)=\int_{0}^{\inf}h_a(t)(\cos \Omega t-j\sin\Omega t)dt$

实际上有： $H_a(-j\Omega)=H_a^*(j\Omega)$

定义模拟滤波器的振幅平方函数为：

$A(\Omega^2)=|H_a(j\Omega)|^2=H_a(j\Omega)H_a^*(j\Omega)$

令： $s=j\Omega$

如果要系统稳定，那么 $A(\Omega^2)=A(-s^2)$

如果我们要让系统函数稳定，就应该选用 $A(-s^2)$ 在s剖面的左半平面的极点作为 $H_a(s)$ 的极点。

来看看今天的内容：

首先，和Jupyter笔记本一样，先导入我们需要的包：

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

1、巴特沃斯滤波器

其振幅平方函数为：

$A(\Omega)=|H_a(j\Omega)|^2=\frac{1}{1+(\frac{j\Omega}{j\Omega_c})^{2N}}=\frac{1}{1+(\Omega/\Omega_c)^{2N}}$

其中，N是滤波器的阶数，N越大，带通和傣族的近似性越好，过渡带也就越陡。

tips：之前一大段时间没有更新，一个是野外没条件，另一个原因就是懒得没有好好去读Scipy.signal的文档，所以说博客有一大段时间空下来了，其实这两天再去读文档，同时对照着Matlab书里面的函数讲解，发现很多都是一样的。

MATLAB中，buttap函数用来计算N阶巴特沃斯归一化，模拟低通原型滤波器系统函数的零点、极点、增益因子的，Python也一样，返回的都是z，p，k，分别是G(p)的极点、零点、增益

我们来看一个最简单的例子：产生一个20阶低通模拟滤波器原型，表示为零极点增益形式：

[z,p,k]=signal.buttap(20)
[n,den]=signal.zpk2tf(z,p,k)
[h,w]=signal.freqs(n,den)
plt.subplot(211)
plt.plot(np.abs(h))
plt.grid(True)
plt.subplot(212)
plt.plot(w)
plt.grid(True)

来看看结果：chule

说实话，除了这张图以外，其他的我都能和Matlab的对的上。

那就来看一看不同阶数下的巴特沃斯滤波器的幅频响应曲线：

n=np.linspace(0,2,200,dtype='float')
[z1,p1,k1]=signal.buttap(1)
[num1,den1]=signal.zpk2tf(z1,p1,k1)
[w1,h1]=signal.freqs(num1,den1)
magh1=abs(h1)
[z2,p2,k2]=signal.buttap(3)
[num2,den2]=signal.zpk2tf(z2,p2,k2)
[w2,h2]=signal.freqs(num2,den2)
magh2=abs(h2)
[z3,p3,k3]=signal.buttap(8)
[num3,den3]=signal.zpk2tf(z3,p3,k3)
[w3,h3]=signal.freqs(num3,den3)
magh3=abs(h3)
[z,p,k]=signal.buttap(12)
[num,den]=signal.zpk2tf(z,p,k)
[w,h]=signal.freqs(num,den)
magh=abs(h)
plt.subplot(2,2,1)
plt.plot(magh1)
plt.grid(True)
plt.subplot(2,2,2)
plt.plot(magh2)
plt.grid(True)
plt.subplot(2,2,3)
plt.plot(magh3)
plt.grid(True)
plt.subplot(2,2,4)
plt.plot(magh)
plt.grid(True)

在已知设计参数 $\omega_p,\omega_s,R_p,R_s$ 之后，利用buttord函数可以求出所需要的滤波器的阶数和3dB截止频率：

[n,Wn]=signal.buttord(Wp,Ws,Rp,Rs)

其中：

Wp:带通截止频率

Ws:带阻起始频率

Rp:通带内波动

Rs:阻带内最小衰减

1、低通滤波器：

# 采样速率为10000H在，设计一个低通滤波器，fp=2000Hz，fs=3000H在，Rp=4dB，Rs=30dB
fn=10000
fp=900
fs=600
Rp=3
Rs=20
Wp=(fp/(fn/2))
Ws=fs/(fn/2)
[n,Wn]=signal.buttord(Wp,Ws,Rp,Rs)
[b,a]=signal.butter(n,Wn)
[H,F]=signal.freqz(b,a,1000,8000)
plt.subplot(211)
plt.plot(H,20*np.log10(abs(F)))
plt.xlabel("frequency")
plt.ylabel('altitude')
plt.title("LowPass")
plt.grid(True)
pha=np.angle(F)*180/np.pi
plt.subplot(212)
plt.plot(H,pha)
plt.xlabel("frequency")
plt.ylabel('angle')
plt.grid(True)

这里有个地方注意一下：Matlab和Python的signal.freqs/z两个函数的输出顺序是不同的，Matlab的输出的H和W和Python输出的H的W两者刚好调换了位置。sheshe

本文章已经生成可运行项目