当前位置: 首页 > news >正文

Matlab 与 Python 基于窗函数的滤波器设计对比 之 凯瑟窗

目录

  • 凯瑟窗的参数估计
    • Matlab 中的 kaiserord
    • Python 中的 kaiserord
    • 区别
  • 基于窗函数的滤波器设计函数
    • Matlab 中的 fir1
    • Python 中的 firwin
    • 区别
  • 滤波器设计举例
    • 低通滤波器设计
    • 高通滤波器设计
    • 带通滤波器设计
    • 带阻滤波器设计

Kaiser 窗近似于扁长椭圆形窗,它使主瓣能量与旁瓣能量之比最大。对于特定长度的 Kaiser 窗,参数 β 控制相对旁瓣衰减。对于给定的 β,相对旁瓣衰减相对于窗长度是固定的。随着 β 的增加,相对旁瓣衰减降低,主瓣宽度增加。(以上内容来自: Kaiser 窗)

凯瑟窗的参数估计

MatlabPython都提供了kaiserord函数, 用于估算特定频率响应条件下的Kaiser窗的阶数(阶数为窗长度-1)和 β \beta β参数.

python中的kaiserord特指scipy.signal模块中的kaiserord函数.

Matlab 中的 kaiserord

Matlab中的函数语法形式如下:

  1. [n,Wn,beta,ftype] = kaiserord(f,a,dev)
  2. [n,Wn,beta,ftype] = kaiserord(f,a,dev,fs)
  3. c = kaiserord(f,a,dev,fs,‘cell’)

参数:

  • f : 通带边界频率阻带截止频率构成的向量, 对于低通高通滤波器, f的长度为2, 对于带通带阻滤波器, f的长度为4; 形式1中的f为数字归一化频率.
  • a: 频带幅度. 低通为[1, 0], 高通为[0, 1], 带通为[0,1,0], 带阻为[1,0,1]
  • dev: 最大允许偏差(线性值).
  • fs: 采样率

返回值:

  • n: 滤波器阶数. n = 窗长度 - 1
  • Wn: 数字归一化频率
  • bata: β \beta β参数
  • ftype: 滤波器类型. 取值: 'low' | 'bandpass' | 'high' | 'stop' | 'DC-0' | 'DC-1'

Python 中的 kaiserord

Python中的函数定义如下:

numtaps, beta = scipy.signal.kaiserord(ripple, width)

参数:

  • ripple: 最大允许偏差(dB)
  • width: 数字归一化的过渡带宽, 等于通带边界频率与阻带截止频率的差.
    返回值:
  • numtaps: 窗的长度.
  • beta: β \beta β参数

区别

因为线性相位特性为Type II(当win_len为偶数时)的滤波器无法实现高通滤波器和带阻滤波器(这里不做论证, 感兴趣的自行查找). 因此使用Matlab中的kaiserord 估算高通和带阻滤波器的阶数时, 阶数n必为偶数(win_len为奇数), 但python中的scipy.signal.kaiserord并不知道要实现何种滤波器, 因此numtaps可能为偶数 (参见:高通滤波器设计示例和带阻滤波器设计示例).

基于窗函数的滤波器设计函数

Matlab中使用fir1函数实现基于窗函数的滤波器设计, python中使用scipy.signal.firwin函数来实现. 这两个函数都能够实现Type I(当win_len为奇数时) 和 Type II(当win_len为偶数时)的滤波器设计.

Matlab 中的 fir1

Matlab中的fir1的语法形式如下:

  • b = fir1(n,Wn,ftype, window, scaleopt)

参数:

  • n: 滤波器阶数
  • Wn: 归一化数字频率
  • ftype: 滤波器类型. 取值: ‘low’ | ‘bandpass’ | ‘high’ | ‘stop’ | ‘DC-0’ | ‘DC-1’
  • window: 窗向量
  • scaleopt: 缩放选项. 取值: ‘scale’ (默认) | ‘noscale’
    • ‘scale’: 对系数进行缩放, 使得通带的幅频响应为 0 dB.
    • ‘noscale’: 不对系数进行缩放.

返回值:

b: 滤波器系数

Python 中的 firwin

python中的firwin的语法形式如下:

scipy.signal.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None)

参数:

  • numtaps: 窗长度.
  • cutoff: 截止频率. 对于低通和高通, 该参数为一个数值;对于带通和带阻, 该参数为一个二元元组.
  • width: 过渡带带宽.
  • window: 窗的名字, 或者包含窗名和参数的元组, 如: ('kaiser', beta).
  • scale: 是否缩放滤波器系数, 使得通带中心频率的幅度响应为0 dB.
  • pass_zero: 表示直流信号(频率为0)是否通过. 取值: TrueFalse, 也可直接给滤波器的类型: { ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’}. 对于低通和带阻, 该参数设为: True, 对于高通和带通, 该参数设为False.

返回值:

  • h: 滤波器系数

区别

Matlab中的fir1函数同python中的scipy.signal.firwin函数完成相同的功能. 就凯瑟窗而言, 主要区别如下:

  • Matlabfir1函数不需要指定滤波器的类型(低通, 高通, 带通或带阻), 因为已经在kaiserord已经指定了滤波器类型. 而python中的firwin函数需要通过pass_zero参数指定滤波器类型.
  • fir1函数使用滤波器的阶数n, 而firwin函数使用窗长度numtaps

滤波器设计举例

以下设计示例中, MatlabPython生成的滤波器系数完全相同.

低通滤波器设计

基于凯瑟窗, 设计一个低通滤波器, 采样频率为1000 Hz, 通带边界频率为 150 Hz, 阻带截止频率为200 Hz, 通带最大允许偏差为0.05, 阻带最大允许偏差为 0.01.

  • 使用Matlab 设计此滤波器
figure;
fs = 1000;
fcuts = [150, 200];
mags = [1, 0];
devs = [0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('低通滤波器')

在这里插入图片描述

  • 使用Python设计此滤波器
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

fs = 1000
fcuts = np.array([150, 200]) / (fs / 2)
devs = np.array([0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = fcuts[1] - fcuts[0]
numtaps, beta = signal.kaiserord(ripple, width)
taps = signal.firwin(numtaps, np.mean(fcuts), window=('kaiser', beta), 
                     scale=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

高通滤波器设计

基于凯瑟窗, 设计一个高通滤波器, 采样频率为1000 Hz, 通带边界频率为 350 Hz, 阻带截止频率为300 Hz, 阻带最大允许偏差为0.05, 通带最大允许偏差为 0.01.

  • 使用Matlab 设计此滤波器
figure;
fs = 1000;
fcuts = [300, 350];
mags = [0, 1];
devs = [0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('高通滤波器')

在这里插入图片描述

  • 使用Python 设计此滤波器
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

fs = 1000
fcuts = np.array([300, 350]) / (fs / 2)
devs = np.array([0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = fcuts[1] - fcuts[0]
numtaps, beta = signal.kaiserord(ripple, width)
if numtaps % 2 == 0:
    numtaps += 1
taps = signal.firwin(numtaps, np.mean(fcuts), window=('kaiser', beta), 
                     scale=False, pass_zero=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

带通滤波器设计

基于凯瑟窗, 设计一个带通滤波器, 采样频率为1000 Hz, 通带边界频率分别为 200 Hz300 Hz, 阻带截止频率分别为150 Hz350 Hz, 通带最大允许偏差为 0.05, 阻带最大允许偏差为 0.01.

  • 使用Matlab 设计此滤波器
figure;
fs = 1000;
fcuts = [150, 200, 300, 350];
mags = [0, 1, 0];
devs = [0.01, 0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('带通滤波器')

在这里插入图片描述

  • 使用Python 设计此滤波器
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

fs = 1000
fcuts = np.array([150, 200, 300, 350]) / (fs / 2)
devs = np.array([0.01, 0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = np.min(fcuts[1::2] - fcuts[0::2])
numtaps, beta = signal.kaiserord(ripple, width)
taps = signal.firwin(numtaps, (fcuts[1] - width / 2, fcuts[2] + width / 2), 
                     window=('kaiser', beta), scale=False, pass_zero=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

带阻滤波器设计

基于凯瑟窗, 设计一个带阻滤波器, 采样频率为1000 Hz, 通带边界频率分别为150 Hz350 Hz, 阻带截止频率分别为 200 Hz300 Hz, 阻带最大允许偏差为 0.05, 通带最大允许偏差为 0.01.

  • 使用Matlab 设计此滤波器
figure;
fs = 1000;
fcuts = [150, 200, 300, 350];
mags = [1, 0, 1];
devs = [0.01, 0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('带阻滤波器')

在这里插入图片描述

  • 使用Python 设计此滤波器
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

fs = 1000
fcuts = np.array([150, 200, 300, 350]) / (fs / 2)
devs = np.array([0.01, 0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = np.min(fcuts[1::2] - fcuts[0::2])
numtaps, beta = signal.kaiserord(ripple, width)
if numtaps % 2 == 0:
    numtaps += 1
taps = signal.firwin(numtaps, (fcuts[1] - width / 2, fcuts[2] + width / 2), 
                     window=('kaiser', beta), scale=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

以上内容为个人的经验总结. 有不当之处, 欢迎留言指正!!!

相关文章:

  • java web开发(从spring boot到spring cloud)
  • 看呆了!二面高德 Java 岗,问了一堆源码,微服务,分布式,Redis,心累
  • 2022华为杯研究生数学建模竞赛B题思路解析
  • 2022华为杯研究生数学建模竞赛E题思路解析
  • 【C语言】学生考勤管理系统
  • 常用的调试技巧(如何检测bug)
  • SpringBoot二十六课大纲和目录
  • 2022年中国研究生数学建模竞赛C题-汽车制造涂装-总装缓存调序区调度优化问题
  • 2022研究生数模A题——移动场景超分辨定位问题
  • vue2脚手架之全局事件总线
  • spring boot学生社团管理系统的设计与实现毕业设计源码151109
  • STM32CUBEMX开发GD32F303(15)----外部中断EXTI
  • 算法竞赛Java选手的语言快速熟悉指南
  • 【电商数仓】数仓搭建之服务数据(data warehouse service-- DWS)层(DWS层概述、几个系统函数和用户主题的建立与数据导入)
  • Window系统安装FFmpeg
  • $translatePartialLoader加载失败及解决方式
  • 「前端」从UglifyJSPlugin强制开启css压缩探究webpack插件运行机制
  • 10个最佳ES6特性 ES7与ES8的特性
  • 11111111
  • 4月23日世界读书日 网络营销论坛推荐《正在爆发的营销革命》
  • ES学习笔记(12)--Symbol
  • HomeBrew常规使用教程
  • JS题目及答案整理
  • Logstash 参考指南(目录)
  • MySQL常见的两种存储引擎:MyISAM与InnoDB的爱恨情仇
  • OSS Web直传 (文件图片)
  • Terraform入门 - 1. 安装Terraform
  • 阿里中间件开源组件:Sentinel 0.2.0正式发布
  • 算法-插入排序
  • 问:在指定的JSON数据中(最外层是数组)根据指定条件拿到匹配到的结果
  • 学习使用ExpressJS 4.0中的新Router
  • 由插件封装引出的一丢丢思考
  • FaaS 的简单实践
  • Java性能优化之JVM GC(垃圾回收机制)
  • 摩拜创始人胡玮炜也彻底离开了,共享单车行业还有未来吗? ...
  • 移动端高清、多屏适配方案
  • ​一些不规范的GTID使用场景
  • # Python csv、xlsx、json、二进制(MP3) 文件读写基本使用
  • (3)选择元素——(14)接触DOM元素(Accessing DOM elements)
  • (阿里云万网)-域名注册购买实名流程
  • (翻译)Quartz官方教程——第一课:Quartz入门
  • (机器学习-深度学习快速入门)第一章第一节:Python环境和数据分析
  • (力扣)1314.矩阵区域和
  • (六)vue-router+UI组件库
  • (牛客腾讯思维编程题)编码编码分组打印下标(java 版本+ C版本)
  • (七)Knockout 创建自定义绑定
  • (转)可以带来幸福的一本书
  • (转)项目管理杂谈-我所期望的新人
  • (转载)(官方)UE4--图像编程----着色器开发
  • .NET/ASP.NETMVC 深入剖析 Model元数据、HtmlHelper、自定义模板、模板的装饰者模式(二)...
  • .NET/C# 将一个命令行参数字符串转换为命令行参数数组 args
  • .Net的DataSet直接与SQL2005交互
  • .NET企业级应用架构设计系列之技术选型
  • .net下的富文本编辑器FCKeditor的配置方法
  • .NET与java的MVC模式(2):struts2核心工作流程与原理