I made a prototype in python that I'm converting to an iOS app. Unfortunately, all the nice features of scipy and numpy are not available in objective-C. So, apparently I need to implement a filter in objective C from scratch. As a first step, I'm trying to implement an IIR from scratch in python. If i can understand how to do it in python, I'll be able to code it in C.
我在python里做了一个原型，我正在把它转换成一个iOS应用程序。不幸的是，在objective-C中，所有的scipy和numpy的好特性都是不可用的。所以，显然我需要在objective C中从头开始实现一个过滤器。作为第一步，我尝试在python中从头开始实现IIR。如果我能理解如何在python中使用它，我将能够在C中编码它。
As a side note, I'd appreciate any suggestions for resources on doing filtering in iOS. As a newbie to objective C who is used to matlab and python, I'm shocked that things like Audio Toolboxes and Accelerate Frameworks and Amazing Audio Engines don't have an equivalent to scipy.signal.filtfilt, nor filter design functions like scipy.signal.butter etc.
作为一个附带说明，我非常感谢在iOS中对资源进行过滤的任何建议。作为一个对objective - C的新手，在matlab和python中使用，我很震惊，像Audio toolbox和加速框架和惊人的音频引擎没有一个等同于scipy.signal。filtfilt，或过滤器设计功能，如scipy.signal。黄油等。
So, in the following code I implement the filter in five ways. 1) scipy.signal.lfilter (for comparison), 2) a state space form using A, B, C, D matrices as calculated by Matlab's butter function. 3) a state space form using the A, B, C, D matrices as calculated by scipy.signal.tf2ss. 4) the Direct Form I, 5) the Direct Form II.
因此，在下面的代码中，我以五种方式实现了过滤器。1)scipy.signal。lfilter(用于比较)，2)用Matlab的butter函数计算的a、B、C、D矩阵的状态空间形式。3)用scipy.signal.tf2ss计算的a、B、C、D矩阵的状态空间形式。直接形式I, 5)直接形式II。
As you can see, the state space form using Matlab matrices works well enough for me to use it in my app. However, I'm still seeking to understand why the others don't work so well.
正如你所看到的，使用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])
And the output:
和输出:
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
2 个解决方案
#1
1
So, I finally found the part of the accelerate framework I was looking for.
所以，我终于找到了我正在寻找的加速框架的一部分。
I was implementing the filter in the first place for downsampling; you need to filter before downsampling to avoid aliasing. Accelerate's vDSP_zrdesamp is the function I wanted all along.
我在第一个地方实现了向下采样的过滤器;在向下采样之前，需要过滤，以避免混叠。加速的vDSP_zrdesamp是我一直想要的功能。
Furthermore, for filtering alone, the ipodEQ audio unit is usually the right choice: (subtype kAudioUnitSubType_AUiPodEQ)
此外，对于单独过滤，ipodEQ音频单元通常是正确的选择:(subtype kAudioUnitSubType_AUiPodEQ)
If you actually need to implement an filter by scratch, the state-space form seems the best.
如果您确实需要从头实现一个过滤器，状态空间形式似乎是最好的。
Still unanswered: why don't my direct form I and II implementations work as intended?
仍然没有回答:为什么我的直接表单I和II实现按预期工作?
#2
0
look at FIR wiki, and this formula:
看看FIR的维基，这个公式:
so first you generate a hamming window yourself (i'm still using python but you can translate it to c/cpp):
首先，你自己生成一个汉明窗口(我还在使用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
then your own low pass filter:
然后你自己的低通滤波器:
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
#1
1
So, I finally found the part of the accelerate framework I was looking for.
所以，我终于找到了我正在寻找的加速框架的一部分。
I was implementing the filter in the first place for downsampling; you need to filter before downsampling to avoid aliasing. Accelerate's vDSP_zrdesamp is the function I wanted all along.
我在第一个地方实现了向下采样的过滤器;在向下采样之前，需要过滤，以避免混叠。加速的vDSP_zrdesamp是我一直想要的功能。
Furthermore, for filtering alone, the ipodEQ audio unit is usually the right choice: (subtype kAudioUnitSubType_AUiPodEQ)
此外，对于单独过滤，ipodEQ音频单元通常是正确的选择:(subtype kAudioUnitSubType_AUiPodEQ)
If you actually need to implement an filter by scratch, the state-space form seems the best.
如果您确实需要从头实现一个过滤器，状态空间形式似乎是最好的。
Still unanswered: why don't my direct form I and II implementations work as intended?
仍然没有回答:为什么我的直接表单I和II实现按预期工作?
#2
0
look at FIR wiki, and this formula:
看看FIR的维基，这个公式:
so first you generate a hamming window yourself (i'm still using python but you can translate it to c/cpp):
首先，你自己生成一个汉明窗口(我还在使用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
then your own low pass filter:
然后你自己的低通滤波器:
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