python实现IIR高通低通，带通，带阻滤波器详解及应用案例

Python实现数字滤波器文章目录Python实现数字滤波器1、IIR低通、高通、带通和带阻滤波器的设计1.1、设计滤波器的函数1.2、将滤波器应用于语音由语音的产生和感知可知，基音频率的范围是60到450赫兹之间，故语音信号采集需要提取基音时，需要采用低通滤波器来获取低频基音信号，在采用计算机采集语音信号时，语音常混有50赫兹交流混音，也需要采用高通滤波器将其去除，此篇设计数字滤波器，以及实现他

树盟

14733人浏览 · 2021-03-04 20:48:10

树盟 · 2021-03-04 20:48:10 发布

Python实现数字滤波器

文章目录

Python实现数字滤波器
- 1、IIR低通、高通、带通和带阻滤波器的设计
- - 1.1、设计滤波器的函数
  - 1.2、将滤波器应用于语音

由语音的产生和感知可知，基音频率的范围是60到450赫兹之间，故语音信号采集需要提取基音时，需要采用低通滤波器来获取低频基音信号，在采用计算机采集语音信号时，语音常混有50赫兹交流混音，也需要采用高通滤波器将其去除，此篇设计数字滤波器，以及实现他们在语言中的简单应用。滤波器传递函数如下：

给定数字滤波器的M阶分子b和N阶分母a：

            jw                 -jw              -jwM
   jw    B(e  )    b[0] + b[1]e    + ... + b[M]e
H(e  ) = ------ = -----------------------------------
            jw                 -jw              -jwN
         A(e  )    a[0] + a[1]e    + ... + a[N]e

给定模拟滤波器的M阶分子b和N阶分母a：

        b[0]*(jw)**M + b[1]*(jw)**(M-1) + ... + b[M]
H(w) = ----------------------------------------------
        a[0]*(jw)**N + a[1]*(jw)**(N-1) + ... + a[N]

1、IIR低通、高通、带通和带阻滤波器的设计

1.1、设计滤波器的函数

在python的scipy.signal库中给出了很多数字信号处理的函数参考：https://docs.scipy.org/doc/scipy/reference/signal.html#module-scipy.signal 可采用以下函数来设计数字滤波器。

（1）巴特沃斯滤波器

巴特沃斯数字和模拟滤波器设计函数：

butter(N, Wn[, btype, analog, output, fs])

功能：设计一个N阶数字或模拟巴特沃斯滤波器，然后返回滤波器系数。

调用格式：

b,a=scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘ba’, fs=None)
z,p,k=scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘zpk’, fs=None)
sos=scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘sos’, fs=None)

参数：

N（int）（滤波器阶次）
Wn（array、int）（临界频率对于低通和高通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于巴特沃斯滤波器，这是增益降至通带的1/sqrt(2)的点“-3dB点”，对于数字滤波器，Wn与fs的单位相同，对于模拟滤波器，Wn是角频率rad / s）
btype：（‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’optional）（滤波类型，默认值是低通其中低通、高通、带通、带阻分别对应‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

b（ndarray）（传递函数分子）
a（ndarray）（传递函数分母）
z（ndarray）（零点）
p（ndarray）（极点）
k（float）（系统增益）
sos（ndarray）（IIR滤波器二阶部分表示）
巴特沃斯滤波器在通带中具有最大平坦的频率响应，如果要求传递函数形式，则由于根与多项式系数之间的转换是数值敏感的运算，即使对于N> = 4，也会出现数值问题。建议使用SOS表示形式。

绘制模拟滤波器频率响应：

计算频率响应的语句：

w, h = signal.freqs(b, a)#计算模拟滤波器的频率响应，返回h为频率响应，w为h的角频率
w, h = signal.freqz(b, a, worN=512, whole=False, plot=None, fs=6.283185307179586, include_nyquist=False)#计算数字滤波器的频率响应，返回h为频率响应（复数表示）w为h的角频率单位与fs相同默认情况下，w归一化为范围[0，pi）（弧度/样本），include_nyquist包含最后一个频率（奈奎斯特频率）
filtered = signal.sosfilt(sos, sig)#使用级联的二阶节沿一维过滤数据

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
b, a = signal.butter(4, 100, 'low', analog=True)  # 4阶低通临界频率为100Hz
w, h = signal.freqs(b, a)  # h为频率响应,w为频率
plt.figure(1)
plt.semilogx(w, 20 * np.log10(abs(h)))  # 用于绘制折线图，两个函数的 x 轴、y 轴分别是指数型的，并且转化为分贝
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)  # 去除画布四周白边
plt.grid(which='both', axis='both')  # 生成网格，matplotlin.pyplot.grid(b, which, axis, color, linestyle, linewidth， **kwargs)， which : 取值为'major', 'minor'， 'both'
plt.axvline(100, color='green')  # 绘制竖线
plt.show()
t = np.linspace(0, 1, 1000, False)  # 1秒，1000赫兹刻度
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)  # 合成信号
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)  # 2行1列的图
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])  # 坐标范围

在这里插入图片描述

生成由10 Hz和20 Hz组成的信号，以1 kHz采样

t = np.linspace(0, 1, 1000, False)  # 1秒，1000赫兹刻度
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)  # 合成信号
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)  # 2行1列的图
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])  # 坐标范围

设计一个15 Hz的数字高通滤波器以消除10 Hz的音调，并将其应用于信号。（建议在过滤时使用二阶节格式，以避免传递函数（ba）格式出现数值错误）：

#sos = signal.butter(10, 15, 'hp', fs=1000, output='sos')  #10阶，15赫兹
sos = signal.butter(10, 15, btype='low', analog=False, output='sos', fs=1000)
filtered = signal.sosfilt(sos, sig)  # 滤波
ax2.plot(t, filtered)
ax2.set_title('After 15 Hz high-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()

在这里插入图片描述

巴特沃斯滤波器阶次选择函数：

buttord(wp, ws, gpass, gstop[, analog, fs])

功能：返回最低阶数字或模拟巴特沃思滤波器的阶次，该滤波器在通带中的损失不超过gpass dB，在阻带中具有至少 gstop dB的衰减。

调用格式：

N, Wn = signal.buttord(wp, ws, gpass, gstop, analog=False, fs=None)

参数：

Wp（float）（通带边缘频率）
Ws（float）（阻带边缘频率）
对于数字滤波器，它们的单位与fs相同。默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。（因此，wp和ws在半周期/样本中。）例如：

低通：wp = 0.2，ws = 0.3
高通：wp = 0.3，ws = 0.2
带通：wp = [0.2，0.5]，ws = [0.1，0.6]
带阻：wp = [0.1，0.6]，ws = [0.2，0.5]

对于模拟滤波器，wp和ws是角频率（例如rad / s）。
gpass（float）（通带中的最大损耗dB）
gstop（float）（阻带中的最小衰减dB）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

N（int）（滤波器阶次）
Wn（array、int）（临界频率对于低通和高通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于巴特沃斯滤波器，这是增益降至通带的1/sqrt(2)的点“-3dB点”，对于数字滤波器，Wn与fs的单位相同，对于模拟滤波器，Wn是角频率rad / s）

设计一个模拟带通滤波器，其通带范围在20到50 rad / s之间的3 dB之内，而在14 rad以下和60 rad / s以上的噪声至少抑制-40 dB。绘制其频率响应，以灰色显示通带和阻带约束。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.buttord([20, 50], [14, 60], 3, 40, True)
b, a = signal.butter(N, Wn, 'band', True)
w, h = signal.freqs(b, a, np.logspace(1, 2, 500))
plt.semilogx(w, 20 * np.log10(abs(h)))  # 用于绘制折线图，两个函数的 x 轴、y 轴分别是指数型的，并且转化为分贝
plt.title('Butterworth bandpass filter fit to constraints')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')  # 显示网格
plt.fill([1,  14,  14,   1], [-40, -40, 99, 99], '0.9', lw=0)  # 区域颜色填充阻带
plt.fill([20, 20,  50,  50], [-99, -3, -3, -99], facecolor='gray', alpha=0.2) # 通带，[x_左下, x_左上,  x_右上, x_右下(左下角顺时针)], [y_左下, y_左上,  y_右上, y_右下]
plt.fill([60, 60, 1e9, 1e9], [99, -40, -40, 99], '0.9', lw=0)  # 阻带
plt.axis([10, 100, -60, 3])  # 画布大小
plt.show()

在这里插入图片描述

（2）切比雪夫I型滤波器

切比雪夫I型数字和模拟滤波器设计函数：

cheby1(N, rp, Wn[, btype, analog, output, fs])

功能：设计一个N阶数字或模拟切比雪夫I型滤波器，然后返回滤波器系数。

调用格式：

b,a=scipy.signal.cheby1(N, rp, Wn, btype=‘low’, analog=False, output=‘ba’, fs=None)
z,p,k=scipy.signal.cheby1(N, rp, Wn, btype=‘low’, analog=False, output=‘zpa’, fs=None)
sos=scipy.signal.cheby1(N, rp, Wn, btype=‘low’, analog=False, output=‘sos’, fs=None)

参数：

N（int）（滤波器阶次）
rp（float）（通带内允许的最大纹波低于单位增益。以分贝指定，为正数。）
Wn（array、int）（临界频率对于低通和高通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于I型滤波器，这是过渡带中增益首先下降到-rp以下的点，对于数字滤波器，Wn与fs的单位相同，对于模拟滤波器，Wn是角频率rad / s）
btype：（‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’optional）（滤波类型，默认值是低通其中低通、高通、带通、带阻分别对应‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’默认值为低通）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

b（ndarray）（传递函数分子）
a（ndarray）（传递函数分母）
z（ndarray）（零点）
p（ndarray）（极点）
k（float）（系统增益）
sos（ndarray）（IIR滤波器二阶部分表示）

Chebyshev I型滤波器最大程度地提高了频率响应的通带和阻带之间的截止速率，但代价是通带中的纹波和阶跃响应中的振铃增加；I型滤波器比II型（cheby2）滚降的速度更快，但II型滤波器的通带中没有任何波动；等纹波通带具有N个最大值或最小值。

设计一个模拟滤波器并绘制其频率响应，显示关键点：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
b, a = signal.cheby1(4, 5, 100, 'low', analog=True)
w, h = signal.freqs(b, a)
plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Chebyshev Type I frequency response (rp=5)')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)  # 去除画布四周白边
plt.grid(which='both', axis='both')
plt.axvline(100, color='green')  # # 绘制竖线，低通临界频率为100Hz
plt.axhline(-5, color='green')  # 绘制横线-rp
plt.show()

在这里插入图片描述

生成由10 Hz和20 Hz组成的信号，以1 kHz采样

t = np.linspace(0, 1, 1000, False)  
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])

设计一个15 Hz的数字高通滤波器以消除10 Hz的音调，并将其应用于信号。（建议在过滤时使用二阶节格式，以避免传递函数（ba）格式出现数值错误）：

sos = signal.cheby1(10, 1, 15, 'hp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 15 Hz high-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()

在这里插入图片描述

切比雪夫I型滤波器阶次选择函数：

cheb1ord(wp, ws, gpass, gstop[, analog, fs])

调用格式：

N, Wn = scipy.signal.cheb1ord(wp, ws, gpass, gstop, analog=False, fs=None)

参数：

Wp（float）（通带边缘频率）
Ws（float）（阻带边缘频率）
对于数字滤波器，它们的单位与fs相同。默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。（因此，wp和ws在半周期/样本中。）例如：

低通：wp = 0.2，ws = 0.3
高通：wp = 0.3，ws = 0.2
带通：wp = [0.2，0.5]，ws = [0.1，0.6]
带阻：wp = [0.1，0.6]，ws = [0.2，0.5]

对于模拟滤波器，wp和ws是角频率（例如rad / s）。
gpass（float）（通带中的最大损耗dB）
gstop（float）（阻带中的最小衰减dB）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

N（int）（滤波器阶次）
Wn（array、int）（切比雪夫固有频率（“ 3dB频率”）用于 cheby1给出滤波器结果。如果指定了fs，则其单位相同，并且还必须将fs传递给cheby1）

设计一个数字低通滤波器，使通带在3 dB以内，高达0.2 *（fs / 2），而在至少0.3dB（fs / 2）以上的频率下，至少要抑制-40 dB。绘制其频率响应，以灰色显示通带和阻带约束。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.cheb1ord(0.2, 0.3, 3, 40)
b, a = signal.cheby1(N, 3, Wn, 'low')
w, h = signal.freqz(b, a)
plt.semilogx(w / np.pi, 20 * np.log10(abs(h)))
plt.title('Chebyshev I lowpass filter fit to constraints')
plt.xlabel('Normalized frequency')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.fill([.01, 0.2, 0.2, .01], [-3, -3, -99, -99], facecolor='gray', alpha=0.2) # 通带，[x_左下, x_左上,  x_右上, x_右下(左下角顺时针)], [y_左下, y_左上,  y_右上, y_右下]
plt.fill([0.3, 0.3,   2,   2], [ 9, -40, -40,  9], '0.9', lw=0)
plt.axis([0.08, 1, -60, 3])
plt.show()

在这里插入图片描述

（3）切比雪夫II型滤波器

切比雪夫II型数字和模拟滤波器设计函数：

cheby2(N, rs, Wn[, btype, analog, output, fs])

功能：设计一个N阶数字或模拟切比雪夫II型滤波器，然后返回滤波器系数。

调用格式：

b,a=scipy.signal.cheby2(N, rs, Wn, btype=‘low’, analog=False, output=‘ba’, fs=None)
z,p,k=scipy.signal.cheby2(N, rs, Wn, btype=‘low’, analog=False, output=‘zpa’, fs=None)
sos=scipy.signal.cheby2(N, rs, Wn, btype=‘low’, analog=False, output=‘sos’, fs=None)

参数：

N（int）（滤波器阶次）
rs（float）（阻带所需的最小衰减。以分贝指定，为正数。）
Wn（array、int）（临界频率对于低通和高通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于II型滤波器，这是过渡带中增益首先达到-rs的点，对于数字滤波器，Wn与fs的单位相同，默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。（因此，Wn为半周期/样本。）对于模拟滤波器，Wn是角频率rad / s）
btype：（‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’optional）（滤波类型，默认值是低通其中低通、高通、带通、带阻分别对应‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’默认值为低通）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

b（ndarray）（传递函数分子）
a（ndarray）（传递函数分母）
z（ndarray）（零点）
p（ndarray）（极点）
k（float）（系统增益）
sos（ndarray）（IIR滤波器二阶部分表示）

Chebyshev II型滤波器最大程度地提高了频率响应的通带和阻带之间的截止速率，但以阻带中的纹波和阶跃响应中的振铃为代价，II型滤波器的滚降速度不如I（cheby1）

设计一个模拟滤波器并绘制其频率响应，显示关键点：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
b, a = signal.cheby2(4, 40, 100, 'low', analog=True)
w, h = signal.freqs(b, a)
plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Chebyshev Type II frequency response (rs=40)')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green')  # cutoff frequency
plt.axhline(-40, color='green')  # -rs
plt.show()

在这里插入图片描述

设计一个17 Hz的数字高通滤波器，以消除10 Hz的音调，并将其应用于信号。（建议在过滤时使用二阶节格式，以避免传递函数（ba）格式出现数值错误）：

t = np.linspace(0, 1, 1000, False)
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])

sos = signal.cheby2(12, 20, 17, 'hp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 17 Hz high-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.show()

在这里插入图片描述

切比雪夫II型滤波器阶次选择函数：

cheb2ord(wp, ws, gpass, gstop[, analog, fs])

调用格式：

N, Wn = scipy.signal.cheb2ord(wp, ws, gpass, gstop, analog=False, fs=None)

功能：

返回最低阶数字或模拟Chebyshev Type II滤波器的阶次，该滤波器在通带中损耗不超过gpass dB，在阻带中至少具有 gstop dB衰减。

参数：

与切比雪夫I大致一致

返回值：

与切比雪夫I大致一致

设计一个数字带阻滤波器，该滤波器能抑制从0.2 *（fs / 2）到0.5 *（fs / 2）的-60 dB，同时保持在0.1 *（fs / 2）以下或0.6 *（fs / 2）以下的3 dB之内。绘制其频率响应，以灰色显示通带和阻带约束。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.cheb2ord([0.1, 0.6], [0.2, 0.5], 3, 60)
b, a = signal.cheby2(N, 60, Wn, 'stop')
w, h = signal.freqz(b, a)
plt.semilogx(w / np.pi, 20 * np.log10(abs(h)))
plt.title('Chebyshev II bandstop filter fit to constraints')
plt.xlabel('Normalized frequency')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.fill([.01, .1, .1, .01], [-3,  -3, -99, -99], '0.9', lw=0)
plt.fill([.2,  .2, .5,  .5], [ 9, -60, -60,   9], '0.9', lw=0)
plt.fill([.6,  .6,  2,   2], [-99, -3,  -3, -99], '0.9', lw=0)
plt.axis([0.06, 1, -80, 3])
plt.show()

在这里插入图片描述

（4）椭圆滤波器

椭圆数字和模拟滤波器设计函数：

ellip(N, rp, rs, Wn[, btype, analog, output, fs])

功能：设计一个N阶数字或模拟切比雪夫II型滤波器，然后返回滤波器系数。

调用格式：

b,a=scipy.signal.ellip(N, rp, rs, Wn, btype=‘low’, analog=False, output=‘ba’, fs=None)
z,p,k=scipy.signal.ellip(N, rp, rs, Wn, btype=‘low’, analog=False, output=‘zpk’, fs=None)
sos=scipy.signal.ellip(N, rp, rs, Wn, btype=‘low’, analog=False, output=‘sos’, fs=None)

参数：

N（int）（滤波器阶次）
rp（int）（通带内允许的最大纹波低于单位增益。以分贝指定，为正数。）
rs（float）（阻带所需的最小衰减。以分贝指定，为正数。）
Wn（array、int）（临界频率对于低通和高通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于椭圆滤波器，这是过渡带中的增益首先下降到*-rp以下的点*，对于数字滤波器，Wn与fs的单位相同。默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。（因此，Wn为半周期/样本。）

对于模拟滤波器，Wn是角频率（例如rad / s）。）
btype：（‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’optional）（滤波类型，默认值是低通其中低通、高通、带通、带阻分别对应‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’默认值为低通）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

b（ndarray）（传递函数分子）
a（ndarray）（传递函数分母）
z（ndarray）（零点）
p（ndarray）（极点）
k（float）（系统增益）
sos（ndarray）（IIR滤波器二阶部分表示）

椭圆滤波器也称为Cauer或Zolotarev滤波器，它以频率波动的通带和阻带之间的跃迁速率最大化，但二者均以脉动为代价，并且阶跃响应中的振铃增加；当rp接近0时，椭圆滤波器变成Chebyshev II型滤波器（cheby2）；当rs接近0时，它成为Chebyshev I型滤波器（cheby1）；当两者都接近0时，它就变成了Butterworth过滤器（butter）；等纹波通带具有N个最大值或最小值（例如，一个5阶滤波器具有3个最大值和2个最小值）；因此，对于奇数阶滤波器，DC增益为单位；对于偶数阶滤波器，DC增益为-rp dB。

设计一个模拟滤波器并绘制其频率响应，显示关键点：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt>>> b, a = signal.ellip(4, 5, 40, 100, 'low', analog=True)
w, h = signal.freqs(b, a)
plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Elliptic filter frequency response (rp=5, rs=40)')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green') # cutoff frequency
plt.axhline(-40, color='green') # rs
plt.axhline(-5, color='green') # rp
plt.show()

在这里插入图片描述

生成由10 Hz和20 Hz组成的信号，以1 kHz采样

t = np.linspace(0, 1, 1000, False)  # 1 second
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])

设计一个17 Hz的数字高通滤波器，以消除10 Hz的音调，并将其应用于信号。（建议在过滤时使用二阶节格式，以避免传递函数（ba）格式出现数值错误）：

sos = signal.ellip(8, 1, 100, 17, 'hp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 17 Hz high-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()

在这里插入图片描述

椭圆滤波器阶次选择函数：

ellipord(wp, ws, gpass, gstop[, analog, fs])

功能：

返回最低阶数字或模拟椭圆滤波器的阶数，该阶数在通带中损失不超过gpass dB，在阻带中至少具有 gstop dB衰减。

调用格式：

N, Wn =scipy.signal.ellipord(wp, ws, gpass, gstop, analog=False, fs=None)

参数：

Wp（float）（通带边缘频率）
Ws（float）（阻带边缘频率）
对于数字滤波器，它们的单位与fs相同。默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。（因此，wp和ws在半周期/样本中。）例如：

低通：wp = 0.2，ws = 0.3
高通：wp = 0.3，ws = 0.2
带通：wp = [0.2，0.5]，ws = [0.1，0.6]
带阻：wp = [0.1，0.6]，ws = [0.2，0.5]

对于模拟滤波器，wp和ws是角频率（例如rad / s）。
gpass（float）（通带中的最大损耗dB）
gstop（float）（阻带中的最小衰减dB）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采用频率）

返回值：

N（int）（滤波器阶次）
Wn（array、int）（切比雪夫固有频率（“ 3dB频率”）用于 ellip给出滤波器结果。如果指定了fs，则其单位相同，并且还必须将fs传递给ellip）

设计一个模拟高通滤波器，使通带在30 rad / s以上的3 dB之内，而在10 rad / s时拒绝-60 dB。绘制其频率响应，以灰色显示通带和阻带约束。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.ellipord(30, 10, 3, 60, True)
b, a = signal.ellip(N, 3, 60, Wn, 'high', True)
w, h = signal.freqs(b, a, np.logspace(0, 3, 500))
plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Elliptical highpass filter fit to constraints')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.fill([.1, 10,  10,  .1], [1e4, 1e4, -60, -60], '0.9', lw=0) # stop
plt.fill([30, 30, 1e9, 1e9], [-99,  -3,  -3, -99], '0.9', lw=0) # pass
plt.axis([1, 300, -80, 3])
plt.show()

在这里插入图片描述

1.2、将滤波器应用于语音

设计一个切比雪夫II型滤波器，他的通带频率Wp为500Hz，阻带频率Ws为750Hz，采样频率为8000Hz，通带最大损耗Rp为3dB，阻带最小衰减为50dB,对bluesky31.wav数据进行滤波。

import numpy as np
import librosa, wave
from scipy import signal
import math
# 导入音频及绘图显示包
import librosa.display
# 导入绘图工作的函数集合
import matplotlib.pyplot as plt

times = librosa.get_duration(filename='./bluesky31.wav')  # 获取音频时长单位为秒
y0, sr = librosa.load('bluesky31.wav', sr=8000, offset=0.0, duration=None)  # 返回音频采样数组及采样率
print("sr=", sr)
print("times=", times)
# 绘图
PointNumbers = int(times * sr) + 1
x1 = np.arange(0, PointNumbers, 1)  # 采样点刻度
x2 = np.arange(0, times, 1 / sr)  # 时间刻度
plt.figure()
plt.xlabel("times")
plt.ylabel("amplitude")
plt.title('bluesky31.wav', fontsize=12, color='black')
plt.plot(x2, y0)
plt.show()


# 去除趋势化信息
def polydetrend(x, fs, m):
    N = len(x)  # 信号长度
    xtrend = np.zeros(len(x))  # 创建0数组
    i = np.arange(0, N, 1)  # 时间刻度
    k = i / fs  # 按信号长度取时间刻度
    a = np.polyfit(k, x, m)  # 在此m为逼近多项式系数的次数，返回逼近后的信号序列多项式拟合系数值
    xtrend[i] = np.polyval(a, k)  # 用系数构成趋势项
    y = x - xtrend  # 从需要去趋势信息的信号中减去趋势项
    return (y, xtrend)


y0, xtrend = polydetrend(y0, sr, 2)
plt.figure()
plt.xlabel("times")
plt.ylabel("amplitude")
plt.title('bluesky31.wav_Elimination_trend', fontsize=12, color='black')
plt.plot(x2, y0)
plt.show()
y1 = y0 - np.mean(y0)  # 消除直流分量
y1 = y1 / np.max(np.abs(y1))  # 幅值归一化
plt.figure()
plt.xlabel("times")
plt.ylabel("amplitude")
plt.title('The normalized bluesky31.wav', fontsize=12, color='black')
plt.plot(x2, y1)
plt.show()

sr= 8000
times= 2.375

在这里插入图片描述

# 设计切比雪夫II型滤波器
fp = 500  # 滤波器通带频率
fs1 = 750  # 滤波器阻带频率
sr2 = sr/2  # 采样频率半周期每样本
Wp = fp/sr2  # 通带频率归一化
Ws = fs1/sr2  # 阻带频率归一化
Rp = 3  # 通带最大损耗
Rs = 50  # 阻带最小衰减
n, Wn = signal.cheb2ord(Wp, Ws, Rp, Rs, analog=False, fs=None)  #当参数未归一化时，fs=sr
b, a = signal.cheby2(n, Rs, Wn, btype='lowpass', analog=False, output='ba', fs=None)
sos = signal.cheby2(n, Rs, Wn, btype='lowpass', analog=False, output='sos', fs=None)
w, h = signal.freqz(b, a ,fs=sr)  # 返回的w单位与参数fs相同
plt.plot(w, 20 * np.log10(abs(h)))
plt.title('Chebyshev Type II frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)  # 去除画布四周白边
plt.grid(which='both', axis='both')#网格
plt.axvline(Wn*sr2, color='green')  # 绘制竖线，低通截止频率(取归一化)
plt.axhline(-Rs, color='green')  # 绘制横线，阻带衰减
plt.fill([fs1,  fs1,  sr2,   sr2], [-Rs, 20,  20, -Rs], '0.9', lw=0)  # 阻带约束
plt.fill([0,  0, fp, fp], [ -100, -Rp, -Rp,   -100], '0.9', lw=0)  # 通带约束
plt.show()

在这里插入图片描述

# 对语音信号滤波
filtered = signal.sosfilt(sos, y1)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(x2, y1)
ax1.set_title('The normalized bluesky31.wav')
ax1.axis([0, times, -1, 1])
ax1.set_ylabel('Amplitude')
ax2.plot(x2, filtered)
ax2.set_title('After 500 Hz cheby2 low-pass filter')
ax2.axis([0, times, -1, 1])
ax2.set_xlabel('Time [seconds]')
ax2.set_ylabel('Amplitude')
plt.show()

在这里插入图片描述

# 信号的语谱
LenFrame=200#每帧长
PatFrame=80#帧移
nfft=512
window = signal.windows.hann(LenFrame)
f1, t1, Zxx1 = signal.stft(y1, fs=sr,window=window,nperseg=LenFrame,noverlap=PatFrame,nfft=nfft)  #返回频率刻度序列（1D张量），时间刻度序列（1D张量），以及filtered的短时傅里叶变换（2D张量）
f2, t2, Zxx2 = signal.stft(filtered, fs=sr,window=window,nperseg=LenFrame,noverlap=80,nfft=nfft) 


fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.pcolormesh(t1, f1, np.abs(Zxx1))  #绘制背景图
ax1.set_title('The normalized bluesky31.wav of Language spectra')
ax1.set_ylabel('Frequncy [Hz]')

ax2.pcolormesh(t2, f2, np.abs(Zxx2))
ax2.set_title('After 500 Hz cheby2 low-pass filter of Language spectra')
ax2.set_xlabel('Time [seconds]')
ax2.set_ylabel('Frequncy [Hz]')
plt.show()