如何实现像 scipy.signal.lfilter 这样的过滤器答案

【问题标题】：How to implement a filter like scipy.signal.lfilter如何实现像 scipy.signal.lfilter 这样的过滤器
【发布时间】：2014-01-21 21:53:25
【问题描述】：

我在 python 中制作了一个原型，并将其转换为 iOS 应用程序。不幸的是，scipy 和 numpy 的所有不错的功能在 Objective-C 中都不可用。所以，显然我需要从头开始在目标 C 中实现一个过滤器。作为第一步，我试图在 python 中从头开始实现 IIR。如果我能理解如何在 python 中做到这一点，我将能够用 C 编写代码。

作为旁注，我很感激有关在 iOS 中进行过滤的资源的任何建议。作为一个习惯于 matlab 和 python 的 Objective C 新手，我很震惊像 Audio Toolboxes 和 Accelerate Frameworks 和 Amazing Audio Engines 这样的东西没有相当于 scipy.signal.filtfilt 的东西，也没有像 scipy 这样的过滤器设计功能。 signal.butter 等。

因此，在以下代码中，我以五种方式实现过滤器。 1）scipy.signal.lfilter（用于比较），2）使用由Matlab的黄油函数计算的A，B，C，D矩阵的状态空间形式。 3) 使用由 scipy.signal.tf2ss 计算的 A、B、C、D 矩阵的状态空间形式。 4) 直接形式 I，5) 直接形式 II。

如您所见，使用 Matlab 矩阵的状态空间形式运行良好，足以让我在我的应用程序中使用它。但是，我仍在寻求了解为什么其他人不能很好地工作。

import numpy as np
from scipy.signal  import butter, lfilter, tf2ss

# building the test signal, a sum of two sines;
N = 32 
x = np.sin(np.arange(N)/6. * 2 * np.pi)+\
    np.sin(np.arange(N)/32. * 2 * np.pi)
x = np.append([0 for i in range(6)], x)

# getting filter coefficients from scipy 
b,a = butter(N=6, Wn=0.5)

# getting matrices for the state-space form of the filter from scipy.
A_spy, B_spy, C_spy, D_spy = tf2ss(b,a)

# matrices for the state-space form as generated by matlab (different to scipy's)
A_mlb = np.array([[-0.4913, -0.5087, 0, 0, 0, 0],
        [0.5087, 0.4913, 0, 0, 0, 0],
        [0.1490, 0.4368, -0.4142, -0.5858, 0, 0],
        [0.1490, 0.4368, 0.5858, 0.4142, 0, 0],
        [0.0592, 0.1735, 0.2327, 0.5617, -0.2056, -0.7944],
        [0.0592, 0.1735, 0.2327, 0.5617, 0.7944, 0.2056]])

B_mlb = np.array([0.7194, 0.7194, 0.2107, 0.2107, 0.0837, 0.0837])

C_mlb = np.array([0.0209, 0.0613, 0.0823, 0.1986, 0.2809, 0.4262])

D_mlb = 0.0296

# getting results of scipy lfilter to test my implementation against
y_lfilter = lfilter(b, a, x)

# initializing y_df1, the result of the Direct Form I method.
y_df1 = np.zeros(6) 

# initializing y_df2, the result of the Direct Form II method.
# g is an array also used in the calculation of Direct Form II
y_df2 = np.array([])
g = np.zeros(6)

# initializing z and y for the state space version with scipy matrices.
z_ss_spy = np.zeros(6)
y_ss_spy = np.array([])

# initializing z and y for the state space version with matlab matrices.
z_ss_mlb = np.zeros(6)
y_ss_mlb = np.array([])


# applying the IIR filter, in it's different implementations
for n in range(N):
    # The Direct Form I
    y_df1 = np.append(y_df1, y_df1[-6:].dot(a[:0:-1]) + x[n:n+7].dot(b[::-1]))

    # The Direct Form II
    g = np.append(g, x[n] + g[-6:].dot(a[:0:-1]))
    y_df2 = np.append(y_df2, g[-7:].dot(b[::-1]))

    # State space with scipy's matrices
    y_ss_spy = np.append(y_ss_spy, C_spy.dot(z_ss_spy) + D_spy * x[n+6])
    z_ss_spy = A_spy.dot(z_ss_spy) + B_spy * x[n+6]

    # State space with matlab's matrices
    y_ss_mlb = np.append(y_ss_mlb, C_mlb.dot(z_ss_mlb) + D_mlb * x[n+6])
    z_ss_mlb = A_mlb.dot(z_ss_mlb) + B_mlb * x[n+6]


# getting rid of the zero padding in the results
y_lfilter = y_lfilter[6:]
y_df1 = y_df1[6:]
y_df2 = y_df2[6:]


# printing the results 
print "{}\t{}\t{}\t{}\t{}".format('lfilter','matlab ss', 'scipy ss', 'Direct Form I', 'Direct Form II')
for n in range(N-6):
    print "{}\t{}\t{}\t{}\t{}".format(y_lfilter[n], y_ss_mlb[n], y_ss_spy[n], y_df1[n], y_df2[n])

还有输出：

lfilter         matlab ss       scipy ss        Direct Form I   Direct Form II
0.0             0.0             0.0             0.0             0.0
0.0313965294015 0.0314090254837 0.0313965294015 0.0313965294015 0.0313965294015
0.225326252712  0.22531468279   0.0313965294015 0.225326252712  0.225326252712
0.684651781013  0.684650012268  0.0313965294015 0.733485689277  0.733485689277
1.10082931381   1.10080090424   0.0313965294015 1.45129994748   1.45129994748
0.891192957678  0.891058879496  0.0313965294015 2.00124367289   2.00124367289
0.140178897557  0.139981099035  0.0313965294015 2.17642377522   2.17642377522
-0.162384434762 -0.162488434882 0.225326252712  2.24911228252   2.24911228252
0.60258601688   0.602631573263  0.225326252712  2.69643931422   2.69643931422
1.72287292534   1.72291129518   0.225326252712  3.67851039998   3.67851039998
2.00953056605   2.00937857026   0.225326252712  4.8441925268    4.8441925268
1.20855478679   1.20823164284   0.225326252712  5.65255635018   5.65255635018
0.172378732435  0.172080718929  0.225326252712  5.88329450124   5.88329450124
-0.128647387408 -0.128763927074 0.684651781013  5.8276996139    5.8276996139
0.47311062085   0.473146568232  0.684651781013  5.97105082682   5.97105082682
1.25980235112   1.25982698592   0.684651781013  6.48492462347   6.48492462347
1.32273336715   1.32261397627   0.684651781013  7.03788646586   7.03788646586
0.428664985784  0.428426965442  0.684651781013  7.11454966484   7.11454966484
-0.724128943343 -0.724322419906 0.684651781013  6.52441390718   6.52441390718
-1.16886662032  -1.16886884238  1.10082931381   5.59188293911   5.59188293911
-0.639469994539 -0.639296371149 1.10082931381   4.83744942709   4.83744942709
0.153883055505  0.154067363252  1.10082931381   4.46863620556   4.46863620556
0.24752293493   0.247568224184  1.10082931381   4.18930262192   4.18930262192
-0.595875437915 -0.595952759718 1.10082931381   3.51735265599   3.51735265599
-1.64776590859  -1.64780228552  1.10082931381   2.29229811755   2.29229811755
-1.94352867959  -1.94338167159  0.891192957678  0.86412577159   0.86412577159

【问题讨论】：

标签： python ios numpy filtering signal-processing

【解决方案1】：

看看FIR wiki，还有这个公式：

所以首先你自己生成一个汉明窗口（我仍在使用 python，但你可以将它翻译成 c/cpp）：

def getHammingWin(n):
    ham=[0.54-0.46*np.cos(2*np.pi*i/(n-1)) for i in range(n)]
    ham=np.asanyarray(ham)
    ham/=ham.sum()
    return ham

然后是你自己的低通滤波器：

def myLpf(data, hamming):

    N=len(hamming)
    res=[]
    for n, v in enumerate(data):
        y=0
        for i in range(N):
            if n<i:
                break
            y+=hamming[i]*data[n-i]
        res.append(y)
    return np.asanyarray(res)
    pass

【讨论】：

实际上我要使用的巴特沃斯滤波器是 IIR，而不是 FIR。不同之处在于 IIR 具有与先前输出相乘的系数 - 先前的 y[n] 值。问题是为什么以我的方式实现 IIR（称为直接形式 I）并没有给出与 scipy.signal.lfilter 相同的结果。
这个汉明窗 FIR 肯定会做一个低通滤波器，对于较大的 N 值具有较低的截止频率。
@JohnSeales 抱歉，我对那个 butterworth 也不是很熟悉，上次我看到它的用法是signal.filtfilt(b, a, data, )，不是signal.lfilter，所以我猜lfilter 只是使用 FIR 计算？
filtfilt 是一个双向滤波器，在正向使用 lfilter 一次，然后反转信号并再次使用 lfilter，最后将信号反转回正确的方向。 lfilter 代表“线性滤波器”，用于诸如 butterworth 或 FIR 之类的 IIR。 b 向量是前向系数，例如您示例中的汉明窗值。 a 向量是反馈系数。如果将 a 设置为标量 1，则它没有反馈并且是 FIR。如果它是一个向量（通常以 1 开头），那么它就是一个 IIR。
它与 scipy 不匹配的原因可能是因为 SciPy 的 lFilter 是直接 II 过滤器吗？ docs.scipy.org/doc/scipy/reference/generated/… : "过滤函数实现为直接II转置结构"

【解决方案2】：

所以，我终于找到了我正在寻找的加速框架部分。

我首先实现了过滤器以进行下采样；您需要在下采样之前进行过滤以避免混叠。 Accelerate 的vDSP_zrdesamp 是我一直想要的功能。

此外，对于单独的过滤，ipodEQ 音频单元通常是正确的选择：（子类型kAudioUnitSubType_AUiPodEQ）

如果您确实需要从头开始实现过滤器，状态空间形式似乎是最好的。

仍然没有答案：为什么我的直接形式 I 和 II 实现不能按预期工作？

【讨论】：

【解决方案3】：

为什么我的直接形式 I 和 II 实现不能按预期工作？

您在 Direct Form I 和 Direct Form II 中的错误可能是由于数值精度问题。下面的代码在 Direct Transpose Form II 中实现了一个过滤器，女巫在数值上更稳定（我在某个地方读过这个，我不记得了）：

d = [0]*4
filtered = [0]

for x in df['dado_sem_spike'].values:

    y    = ((b[0] * x)              + d[3]) / a[0]
    d[3] =  (b[1] * x) - (a[1] * y) + d[2]
    d[2] =  (b[2] * x) - (a[2] * y) + d[1]
    d[1] =  (b[3] * x) - (a[3] * y) + d[0]
    d[0] =  (b[4] * x) - (a[4] * y)
    filtered.append(y)

我实现了这个表单，因为我的直接表单没有给出好的结果。

【讨论】：