python - 如何实现像 scipy.signal.lfilter 这样的过滤器

Question

我用 python 制作了一个原型(prototype)，然后将其转换为 iOS 应用程序。不幸的是，scipy 和 numpy 的所有优秀特性在 objective-C 中都不可用。所以，显然我需要从头开始在 objective-c 中实现一个过滤器。作为第一步，我尝试在 python 中从头开始实现 IIR。如果我能理解如何在 python 中执行此操作，我将能够在 C 中对其进行编码。

作为旁注，对于在 iOS 中进行过滤的资源的任何建议，我将不胜感激。作为习惯使用 matlab 和 python 的 objective-c 的新手，我很震惊，音频工具箱和加速框架和惊人的音频引擎之类的东西没有与 scipy.signal.filtfilt 等效的东西，也没有像 scipy 这样的过滤器设计功能。 signal.butter 等

所以，在下面的代码中，我以五种方式实现了过滤器。 1) scipy.signal.lfilter(用于比较)，2) 使用由 Matlab 的 butter 函数计算的 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

Answer 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