滤波器特性分析_信号与系统作业

filter.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
# 信号与系统课程设计：滤波器特性分析
# 作者：马晓忠（Blanboom）
# http://blanboom.org
# 2014 年 06 月 17 日

# 环境：
#       Python:     2.7.5
#       NumPy:      1.6.2
#       SciPy:      0.11.0
#       matplotlib: 1.1.1

# 程序后面部分运行较慢，请耐心等待

import scipy.signal as signal
import pylab as pl
import numpy as np

pl.rcParams['font.sans-serif'] = ['Hiragino Sans GB'] # 指定默认字体
pl.rcParams['axes.unicode_minus'] = False             # 解决保存图像负号'-'显示为方块的问题


# 函数：drawfilter()
# 功能：生成指定参数的滤波器，通过 freqz() 函数得到该滤波器的频率响应曲线
# 输入：
#       wp：通带频率范围，范围为 0 ~ 1，单位为 4kHz
#       ws：阻带频率范围，范围为 0 ~ 1，单位为 4kHz
#       ftype：滤波器类型，'butter' 为巴特沃斯滤波器，'cheby1' 为I切比雪夫滤波器，'cheby2'为II型切比雪夫滤波器，'ellip' 为椭圆滤波器
#       title：生成的频率响应曲线的名称
# 注意：
#       1. 使用此函数并不直接显示图形，而是需要通过 pl.show() 显示. 也可通过 pl.subplot() 等使多张图表显示在同一窗口内
#       2. 此函数对应的滤波器通带增益最大衰减值为 2dB，阻带增益最大衰减值为 60dB
def drawfilter(wp, ws, ftype, title):
    # 生成一个滤波器
    b, a = signal.iirdesign(wp, ws, 2, 60, False, ftype)

    # 计算滤波器的频率响应
    w, h = signal.freqz(b, a)

    # 计算增益
    power = 20 * np.log10(np.clip(np.abs(h), 1e-8, 1e100))

    # 画出图像
    pl.plot(w/np.pi*4000, power, linewidth = 2.4)
    pl.title(title)
    pl.ylabel(u"增益(dB)")
    pl.xlabel(u"频率(Hz)")
    pl.ylim(-165,5)

# ---------------------------------- 低通滤波器的特性分析 ---------------------------------- #
# 通带为 < 0.2*4kHz，阻带为 > 0.4*4kHz

# 调整图表格式
pl.figure(figsize=(12, 8))
pl.subplots_adjust(left = 0.08)
pl.subplots_adjust(bottom = 0.07)
pl.subplots_adjust(right = 0.94)
pl.subplots_adjust(wspace = 0.29)
pl.subplots_adjust(hspace = 0.39)
pl.subplots_adjust(top = 0.95)

# 开始生成并显示图像
pl.subplot(2, 2, 1)
drawfilter(0.2, 0.4, 'butter', u"巴特沃斯低通滤波器")
pl.subplot(2, 2, 2)
drawfilter(0.2, 0.4, 'cheby1', u"I型切比雪夫低通滤波器")
pl.subplot(2, 2, 3)
drawfilter(0.2, 0.4, 'cheby2', u"II型切比雪夫低通滤波器")
pl.subplot(2, 2, 4)
drawfilter(0.2, 0.4, 'ellip', u"椭圆函数低通滤波器")
pl.show()
# --------------------------------------------------------------------------------------- #


# ---------------------------------- 高通滤波器的特性分析 ---------------------------------- #
# 通带为 > 0.4*4kHz，阻带为 < 0.2*4kHz

# 调整图表格式
pl.figure(figsize=(12, 8))
pl.subplots_adjust(left = 0.08)
pl.subplots_adjust(bottom = 0.07)
pl.subplots_adjust(right = 0.94)
pl.subplots_adjust(wspace = 0.29)
pl.subplots_adjust(hspace = 0.39)
pl.subplots_adjust(top = 0.95)

# 开始生成并显示图像
pl.subplot(2, 2, 1)
drawfilter(0.4, 0.2, 'butter', u"巴特沃斯高通滤波器")
pl.subplot(2, 2, 2)
drawfilter(0.4, 0.2, 'cheby1', u"I型切比雪夫高通滤波器")
pl.subplot(2, 2, 3)
drawfilter(0.4, 0.2, 'cheby2', u"II型切比雪夫高通滤波器")
pl.subplot(2, 2, 4)
drawfilter(0.4, 0.2, 'ellip', u"椭圆函数高通滤波器")
pl.show()
# --------------------------------------------------------------------------------------- #


# ---------------------------------- 带通滤波器的特性分析 ---------------------------------- #
# 通带为 0.2*4kHz ~ 0.5*4kHz，阻带为 < 0.1*4kHz  >0.6*4kHz

# 调整图表格式
pl.figure(figsize=(12, 8))
pl.subplots_adjust(left = 0.08)
pl.subplots_adjust(bottom = 0.07)
pl.subplots_adjust(right = 0.94)
pl.subplots_adjust(wspace = 0.29)
pl.subplots_adjust(hspace = 0.39)
pl.subplots_adjust(top = 0.95)

# 开始生成并显示图像
pl.subplot(2, 2, 1)
drawfilter([0.2, 0.5], [0.1, 0.6], 'butter', u"巴特沃斯带通滤波器")
pl.subplot(2, 2, 2)
drawfilter([0.2, 0.5], [0.1, 0.6], 'cheby1', u"I型切比雪夫带通滤波器")
pl.subplot(2, 2, 3)
drawfilter([0.2, 0.5], [0.1, 0.6], 'cheby2', u"II型切比雪夫带通滤波器")
pl.subplot(2, 2, 4)
drawfilter([0.2, 0.5], [0.1, 0.6], 'ellip', u"椭圆函数带通滤波器")
pl.show()
# --------------------------------------------------------------------------------------- #


# ---------------------------------- 带阻滤波器的特性分析 ---------------------------------- #
# 阻带为 0.2*4kHz ~ 0.5*4kHz，通带为 < 0.1*4kHz  >0.6*4kHz

# 调整图表格式
pl.figure(figsize=(12, 8))
pl.subplots_adjust(left = 0.08)
pl.subplots_adjust(bottom = 0.07)
pl.subplots_adjust(right = 0.94)
pl.subplots_adjust(wspace = 0.29)
pl.subplots_adjust(hspace = 0.39)
pl.subplots_adjust(top = 0.95)

# 开始生成并显示图像
pl.subplot(2, 2, 1)
drawfilter([0.1, 0.6], [0.2, 0.5], 'butter', u"巴特沃斯带阻滤波器")
pl.subplot(2, 2, 2)
drawfilter([0.1, 0.6], [0.2, 0.5], 'cheby1', u"I型切比雪夫带阻滤波器")
pl.subplot(2, 2, 3)
drawfilter([0.1, 0.6], [0.2, 0.5], 'cheby2', u"II型切比雪夫带阻滤波器")
pl.subplot(2, 2, 4)
drawfilter([0.1, 0.6], [0.2, 0.5], 'ellip', u"椭圆函数带阻滤波器")
pl.show()
# --------------------------------------------------------------------------------------- #


# -------------------------- 通过 8kHz 频率扫描波得到带阻滤波器的频谱 ------------------------ #
# 这段程序产生的图形应该与前面程序产生的基本一致

# 生成一个巴特沃斯带阻滤波器，参数与前面相同
b, a = signal.iirdesign([0.1, 0.6], [0.2, 0.5], 2, 60, False, 'butter')

# 产生 8kHz 的频率扫描信号，开始频率为 0，结束频率为 4kHz
t = np.arange(0, 2, 1/8000.0)
sweep = signal.chirp(t, f0 = 0, t1 = 2, f1 = 4000.0)

# 将频率扫描信号进行滤波
out = signal.lfilter(b, a, sweep)

# 将波形转换为能量
out = 20 * np.log10(np.abs(out))

# 找到所有局部最大值的下标
index = np.where(np.logical_and(out[1:-1] > out[:-2], out[1:-1] > out[2:]))[0] + 1

# 画出图像
pl.plot(t[index]/2.0*4000, out[index], linewidth = 1.1)
pl.title(u"通过频率扫描波得到的巴特沃斯带阻滤波器频谱")
pl.ylabel(u"增益(dB)")
pl.xlabel(u"频率(Hz)")
pl.show()
# --------------------------------------------------------------------------------------- #


# ------------------------------------- 滤波器的级联 -------------------------------------- #
# 巴特沃斯低通滤波器，通带为 < 0.6*4kHz，阻带为 > 0.7*4kHz
b1, a1 = signal.iirdesign(0.6, 0.7, 2, 60, False, 'butter')

# 巴特沃斯高通滤波器，通带为 > 0.1*4kHz，阻带为 < 0.2*4kHz
b2, a2 = signal.iirdesign(0.2, 0.1, 2, 60, False, 'butter')

# 计算滤波器的频率响应
w1, h1 = signal.freqz(b1, a1)
w2, h2 = signal.freqz(b2, a2)

# 将 Python 数组转换为 Numpy 数组，加快数据处理速度
h1_numpy = np.array(h1)
h2_numpy = np.array(h2)

# 两个滤波器的频率响应相乘，得到级联后滤波器的频率响应
h_numpy = h1_numpy * h2_numpy

# 计算增益
power1 = 20*np.log10(np.clip(np.abs(h1_numpy),      1e-8, 1e100))
power2 = 20*np.log10(np.clip(np.abs(h2_numpy),      1e-8, 1e100))
power  = 20*np.log10(np.clip(np.abs(h_numpy),       1e-8, 1e100))

# 画出图像
pl.figure(figsize=(15, 6))
pl.subplots_adjust(left = 0.06)
pl.subplots_adjust(bottom = 0.12)
pl.subplots_adjust(right = 0.96)
pl.subplots_adjust(wspace = 0.35)
pl.subplots_adjust(top = 0.92)
pl.subplot(1, 3, 1)
pl.plot(w1/np.pi*4000, power1, linewidth = 2.4)
pl.title(u"巴特沃斯低通滤波器频谱")
pl.ylabel(u"增益(dB)")
pl.xlabel(u"频率(Hz)")
pl.ylim(-165,5)
pl.subplot(1, 3, 2)
pl.plot(w2/np.pi*4000, power2, linewidth = 2.4)
pl.title(u"巴特沃斯高通滤波器频谱")
pl.ylabel(u"增益(dB)")
pl.xlabel(u"频率(Hz)")
pl.ylim(-165,5)
pl.subplot(1, 3, 3)
pl.plot(w1/np.pi*4000, power, linewidth = 2.4)
pl.title(u"两个滤波器级联后的频谱")
pl.ylabel(u"增益(dB)")
pl.xlabel(u"频率(Hz)")
pl.ylim(-165,5)
pl.show()
# --------------------------------------------------------------------------------------- #

########################################## EOF ############################################