【matlab 代码的python复现】 Matlab实现的滤波器设计实现与Python 的库函数相同实现Scipy

实现一个IIR滤波器的设计

背景

Matlab 设计的滤波器通常封装过于完整,虽然在DSP中能够实现更多功能的滤波器设计但是很难实现Python端口的实现。
我们以一段原始的生物电信号EEG信号进行处理。

EEG信号

1.信号获取

EEG信号通常通过头皮电极,经过多通道采样芯片采样,将获取到的微弱的头皮模拟电信号,转化为数字信号。获取到的数字信号在处理,通常需要数字信号处理技术与脑机接口技术。例如,滤波、通道选择、频谱分析等。我们通过电极帽获取到的数字信号,按照电路的放大比例进行的放大。因此获取到的信号可能是存在基线漂移的信号,或者非零点的信号。这些都是相对于设备采集来说的。但是信号中存在的特征往往存在是主要的。我们复现Matlab的高质量设计的滤波器,python实现更有助于工程实现。

2.Matlab滤波器设计函数

打开Matlab的信号分析器。这是一个功能强大的信号可视化APP,可以完成简单的信号的预处理的工作。
在这里插入图片描述
导入的是一个1000个采样点的EEg原始信号,这个原始信号设置采样频率为250Hz。经过分析我们得到该信号幅值位于-900,并非处于基线附近因此我们认识到,这个信号存在基线漂移问题。我们打算通过高通或者带通滤波器实现解决。
![在这里插入图片描述](https://img-blog.csdnimg.cn/direct/21b8045c82f6468b833cdb0af601be11.png

1点击显示按钮中的频谱和时频率图像,有助于我们直接观察信号的原始成分。

频谱显示的是功率谱图像,区别是频谱图他更加的光滑。具体区别见这里功率谱,功率、功率密度http://t.csdnimg.cn/Qzldh

能看到低频能量特别高,还有50Hz、100Hz的中间的工频噪声,以及高频噪声。

2点击分析器,设置带通滤波器
在这里插入图片描述
3设置滤波范围6,48Hz
在这里插入图片描述

此时此刻我们观察到这个波形是去除基线漂移并且能过去除高频的波形。位于0基线附近。

点击分析器生成函数,matlab 在工作台生成一个预处理的函数。
在这里插入图片描述
在这里插入图片描述

data =load("output.mat")
data2 =squeeze(data.data(1,1,1,:));

% 定义采样频率 (Hz)
fs = 250;

% 定义采样点数
num_samples = 1000;

% 计算采样周期 (seconds)
Ts = 1/fs;

% 创建时间向量tx,范围从0到( num_samples - 1 ) * Ts
tx = (0 : num_samples - 1) * Ts;
y = preprocess(data2,tx);
plot(y)
hold on
plot(data2)

运行结果显示出来/在这里插入图片描述
这是一个我认为完美的信号。

通过了解Matlab使用的是IIR auto设计的滤波系数。因此我尝试用Python实现。

3 基于Python 的Scipy.io实现滤波

def preprocess_input(x, tx):
    """
    预处理输入信号 x,根据时间向量 tx 进行带通滤波。

    参数:
        x (array_like): 输入信号数据.
        tx (array_like): 时间向量(单位:秒).

    返回:
        array_like: 经过带通滤波处理后的信号数据.
    """
    # 计算平均采样率
    import numpy as np
    from scipy.signal import ellip, butter, lfilter, freqz, iirdesign, cheby2, buttord, cheb2ord

    Wpass1, Wpass2 = (0.0480, 0.3840)
    Wstop1, Wstop2 = (0.0405, 0.3915)
    Apass, Astop = 0.1, 60

    b, a = iirdesign(wp=[Wpass1, Wpass2], ws=[Wstop1, Wstop2],
                     gpass=Apass, gstop=Astop,
                     analog=False, ftype='elliptic')
    # 设计滤波器
    y = filtfilt(b, a, x)

    return y


# 示例使用参数
import numpy as np
data =np.load(r"2024-04-09 14_36_5_eegdata.npy")
print(data.shape

最近更新

  1. TCP协议是安全的吗?

    2024-04-23 00:06:03       16 阅读
  2. 阿里云服务器执行yum,一直下载docker-ce-stable失败

    2024-04-23 00:06:03       16 阅读
  3. 【Python教程】压缩PDF文件大小

    2024-04-23 00:06:03       15 阅读
  4. 通过文章id递归查询所有评论(xml)

    2024-04-23 00:06:03       18 阅读

热门阅读

  1. Vue 组件通信的几种方式

    2024-04-23 00:06:03       12 阅读
  2. C++:异常处理

    2024-04-23 00:06:03       13 阅读
  3. 计算机网络——应用层(3)电子邮件

    2024-04-23 00:06:03       13 阅读
  4. .net core8 自定义一个中间件

    2024-04-23 00:06:03       13 阅读
  5. node.js 什么是模板引擎?(具体介绍underscore)

    2024-04-23 00:06:03       13 阅读
  6. Android R framework修改低电量关机值为2%

    2024-04-23 00:06:03       14 阅读
  7. 信息物理系统技术概述_1.概念和实现

    2024-04-23 00:06:03       40 阅读
  8. MongoDB 与MySQL的区别?优势?

    2024-04-23 00:06:03       11 阅读
  9. Flume

    Flume

    2024-04-23 00:06:03      39 阅读
  10. HCIP-Datacom-ARST必选题库_36_加密算法【1道题】

    2024-04-23 00:06:03       16 阅读
  11. 【centso】sqlite3.7.17升级到更新的版本

    2024-04-23 00:06:03       22 阅读