信息处理 - Python 和 MATLAB filtfilt 函数的区别 - 吾爱随笔录

Python 和 MATLAB filtfilt 函数的区别

信息处理 matlab 过滤器 Python

2022-01-16 06:52:48

我正在尝试将一些 MATLAB 代码移植到 Python 并遇到一些奇怪的行为。我正在实现一个 5 阶巴特沃斯带通滤波器。采样率为 30 Hz。

在 Windows 7 x64 上运行 MATLAB R2012b、Spyder 2.2.0 和 Python 2.7、SciPy 0.12.0。

在 MATLAB 中：

[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

这是原始信号：

原始信号

和滤波后的信号：

MATLAB 滤波信号

及其功率谱：

MATLAB 功率谱

这是有道理的，因为归一化的带通频率是 0.05 - 0.33。

我之前发现 SciPy 的 butter 函数没有给出与 MATLAB 相同的系数，所以我使用 hdf5write 将滤波器系数从 MATLAB 导出到 Python（参见此处：https ://stackoverflow.com/questions/7117797/export-matlab-variable-to -text-for-python-usage )

在 Python 中：

y = signal.filtfilt(b, a, input_signal, padtype = None)

输出是：

Python过滤信号

及其功率谱：

Python 功率谱

我使用了 padtype = None，因为默认情况下它是 padtype = 'odd'。但是，我尝试了所有不同的填充选项，它们看起来都差不多。

我不完全确定出了什么问题......任何帮助将不胜感激。

编辑：添加了 padtype = "odd" 和 b 的图形以及使用的滤波器系数

padtype = "odd" 的信号：

padtype = 奇数信号

似乎有一些潜在的信号，所以我不认为这是脉冲响应，尽管一开始的大瞬态很奇怪。

Python 滤波器系数（5 阶巴特沃斯滤波器）：

passband = [0.75*2/30, 5.0*2/30]
b, a = scipy.signal.butter(5, passband, 'bandpass')

b, a 是 float64 类型的数组。

b = array([  5.49209388e-03,   0.00000000e+00,  -2.74604694e-02,
    -1.97776791e-17,   5.49209388e-02,   2.47220989e-17,
    -5.49209388e-02,  -1.97776791e-17,   2.74604694e-02,
     0.00000000e+00,  -5.49209388e-03])
a = array([  1.        ,  -6.52098852,  19.50384534, -35.47189804,
    43.65758795, -38.07760914,  23.83047021, -10.55670367,
     3.16710078,  -0.58122912,   0.04945954])

MATLAB 滤波器系数

b =   0.0055, 0, -0.0275, 0, 0.0549, 0, -0.0549, 0, 0.0275, 0, -0.0055
a =   1.0000, -6.5210, 19.5038, -35.4719, 43.6576, -38.0776, 23.8305, 
      -10.5567, 3.1671, -0.5812, 0.0495

结果系数大致相同，我相信我之前遇到的问题是由于频率要高得多（在 MATLAB 和 Python 中必须以不同的方式生成系数）。

输入信号是一个 float32 列表，尽管我尝试使用 numpy.array 转换为数组，结果是一样的。

编辑：有关 padtypes 和 Python 与 MATLAB 的更多信息

Python SciPy 的 filtfilt 函数包括一个名为 padtype 的参数，它指示在信号两侧扩展的填充类型。该填充用于减少瞬变。奇数和偶数是这些扩展与端点的对称类型的描述符（http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html）。

padtypes差异的一些插图：

padtype =“偶数”

均匀填充

padtype =“奇数”

奇数填充

padtype = 无

没有填充

根据结果，MATLAB 似乎使用了奇数填充（这也是 Python 的默认值）：

MATLAB

4个回答

根据新闻组https://mail.python.org/pipermail/scipy-user/2014-April/035648.html，

不同之处在于默认填充长度。在matlab的filtfilt中为3*(max(len(a), len(b)) - 1)，在scipy的filtfilt中为3*max(len(a), len(b))。

因此，我在 Python 中调用 filtfilt 时修改了 padlen。MATLAB 和 Python 之间的结果将是相同的。

MATLAB代码

input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

带有 SciPy 代码的 Python

from scipy.signal import butter, filtfilt
input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal, padtype = 'odd', padlen=3*(max(len(b),len(a))-1))

这个对我有用：

from scipy.io import loadmat
from scipy.signal import butter, filtfilt
from matplotlib.pyplot import plot

signaldata = loadmat('signaldata.mat')

input_signal = signaldata['input_signal'][0]

passband = [0.75*2/30, 5.0*2/30]
b, a = butter(5, passband, 'bandpass')

y = filtfilt(b, a, input_signal)
plot(y)

在此处输入图像描述

我不知道你的问题是什么，但也许你可以从这里解决。

我在 Paco 的示例中得到了相同的结果，但是将输入信号更改为向量并不会在 MATLAB 和 Python（Scipy 1.0.0）中产生相同的结果：

MATLAB：

input_signal = ones(10000,1);
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal);
y(1:5)'

ans =

    1.0e-13 *

    -0.306211857054885  -0.468000282200120  -0.612592864567431  -0.736885525141862  -0.838459113089290

Python：

from numpy import ones
from scipy.signal import butter, filtfilt
input_signal = ones((10000,))
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal, padtype = 'odd', padlen=3*(max(len(b),len(a))-1))
y[:5]

Out[49]: 
array([-2.52233927e-14, -1.79580673e-14, -1.07857124e-14, -3.97094772e-15, 2.30161154e-15])

经过多次试验和错误，我发现以下两个语句是等效的：

% Matlab
rf_data_filt = filtfilt(b, a, rf_data);

# Python
rf_data_filt = filtfilt(b, a, rf_data, axis=0, padtype='odd', padlen=3*(max(len(b),len(a))-1))

似乎需要三个 kwargs：axis=0, padtype='odd', padlen=3*(max(len(b),len(a))-1)

其它你可能感兴趣的问题

上一篇计算离散信号的 n 导数的方法下一篇音频指纹问题