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

一 实现一个IIR滤波器的设计

背景

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

EEG信号

1.信号获取

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

2.Matlab滤波器设计函数

打开Matlab的信号分析器。这是一个功能强大的信号可视化APP,可以完成简单的信号的预处理的工作。

导入的是一个1000个采样点的EEg原始信号，这个原始信号设置采样频率为250Hz。经过分析我们得到该信号幅值位于-900,并非处于基线附近因此我们认识到，这个信号存在基线漂移问题。我们打算通过高通或者带通滤波器实现解决。

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)
sampling_rate = 1000.0  # 假设采样率为1000 Hz
signal =data[0,0,0,:]
time_vector = np.linspace(0, 4, 1000, endpoint=False)  # 示例时间向量
filtered_signal = preprocess_input(signal, time_vector)

#将原始信号与滤波后的信号画图
import matplotlib.pyplot as plt
plt.plot(signal, label='Original Signal')
plt.plot(filtered_signal, label='Filtered Signal')
plt.legend()
plt.show()
# 现在 `filtered_signal` 是经过带通滤波处理后的信号数据

画出的结果是这样的。

为什么使用了相同的参数输出的结果为什么和matlab差别如此之大？
为什么会出现初始的信号的阶跃震荡？
Matlab自动实现的去除基线漂移的效果为什么Python不能是实现？

4.问题分析

我们使用Matlab设计零相位滤波器带通滤波器

load noisysignals x;                    % noisy waveform
x=x-500;
[b,a] = butter(1,[24/1000,48/250],'bandpass');           % IIR filter design
y = filtfilt(b,a,x);                    % zero-phase filtering
y2