MATLAB 的 sgolay(k, f) 的 Python 等效项是什么？

Question

我在 MATLAB 中有一个函数

[b,g] = sgolay(k, f);

它输出一个 f x f 矩阵。

当我在 Python 中为相同的 k 和 f 值运行相同的值时，使用：

scipy.signal.savgol_coeffs(f, k)

它输出一个完全不同的数组，只有 f 个元素。

考虑的值是：

k = 4，f = 21

savol_filter() 接受三个参数，包括数组，而sgolay() 只接受两个。此外，savo_coeffs 没有生成所需的矩阵。

在matlab中获取sgolay(k, f)生成的矩阵的Python等价物是什么？

Answer 1

如果您检查 Matlab 的 sgolay 函数返回的矩阵 b，您会发现 中心行 与 SciPy 的 savgol_coeffs 返回的一维数组相同。 b 的上半部分和下半部分，每个部分都有(framelen - 1)/2 行，是要应用于信号末端的 Savitzky-Golay 滤波器的系数，其中滤波器不对称。也就是说，对于信号每一端的(framelen - 1)/2 值，每个滤波后的值都是使用一组不同的系数计算的。

您可以使用savgol_coeffs 通过迭代pos 参数来生成b。以下 ipython 会话显示了一个示例。

In [74]: import numpy as np

In [75]: from scipy.signal import savgol_coeffs

In [76]: np.set_printoptions(precision=11, linewidth=90)

In [77]: order = 3

In [78]: windowlen = 5

这些是对称（即居中）Savitzky-Golay 滤波器的系数。一维数组应该匹配sgolay返回的矩阵的中心行：

In [79]: savgol_coeffs(windowlen, order)
Out[79]: array([-0.08571428571,  0.34285714286,  0.48571428571,  0.34285714286, -0.08571428571])

如果我们设置pos=windowlen-1，我们会得到设计用于评估窗口一端的滤波器的系数。这些应该匹配sgolay返回的数组的第一行：

In [80]: savgol_coeffs(windowlen, order, pos=windowlen-1)
Out[80]: array([ 0.98571428571,  0.05714285714, -0.08571428571,  0.05714285714, -0.01428571429])

同样，pos=0 给出了窗口另一端的系数。这些应该匹配sgolay返回的矩阵的最后一行：

In [81]: savgol_coeffs(windowlen, order, pos=0)
Out[81]: array([-0.01428571429,  0.05714285714, -0.08571428571,  0.05714285714,  0.98571428571])

这是匹配Matlab的sgolay返回值的完整数组：

In [82]: b = np.array([savgol_coeffs(windowlen, order, pos=p) for p in range(windowlen-1, -1, -1)])

In [83]: b
Out[83]: 
array([[ 0.98571428571,  0.05714285714, -0.08571428571,  0.05714285714, -0.01428571429],
       [ 0.05714285714,  0.77142857143,  0.34285714286, -0.22857142857,  0.05714285714],
       [-0.08571428571,  0.34285714286,  0.48571428571,  0.34285714286, -0.08571428571],
       [ 0.05714285714, -0.22857142857,  0.34285714286,  0.77142857143,  0.05714285714],
       [-0.01428571429,  0.05714285714, -0.08571428571,  0.05714285714,  0.98571428571]])

如果您将此与 Matlab 中 b = sgolay(3, 5) 的结果进行比较，您会发现它们是相同的。

要获得sgolay 返回的g 矩阵，您必须调用savgol_coeffs 并将deriv 设置为range(order+1) 中的值，反转和转置数组，并按阶乘缩放导数顺序。要反转系数，您可以使用::-1 形式的切片，也可以使用savgol_coeffs 的use 选项。

这是使用savgol_coeffs 生成带有order=3 和windowlen=5 的g 矩阵的一种方法：

In [12]: import numpy as np

In [13]: from scipy.signal import savgol_coeffs

In [14]: from scipy.special import factorial

In [15]: np.set_printoptions(precision=11, linewidth=90, suppress=True)

In [16]: order = 3

In [17]: windowlen = 5

In [18]: g = np.array([savgol_coeffs(windowlen, order, deriv=d, use='dot') for d in range(order+1)]).T / factorial(np.arange(order+1))

In [19]: g
Out[19]: 
array([[-0.08571428571,  0.08333333333,  0.14285714286, -0.08333333333],
       [ 0.34285714286, -0.66666666667, -0.07142857143,  0.16666666667],
       [ 0.48571428571,  0.           , -0.14285714286,  0.           ],
       [ 0.34285714286,  0.66666666667, -0.07142857143, -0.16666666667],
       [-0.08571428571, -0.08333333333,  0.14285714286,  0.08333333333]]) 

您没有说明为什么要在 Python 中使用完整的 windowlen x windowlen 数组。你不需要它来使用savgol_filter。