John Seales John Seales - 1 year ago 83
Python Question

How to implement a filter like scipy.signal.lfilter

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.

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.

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.

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.

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

Answer Source

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.

Furthermore, for filtering alone, the ipodEQ audio unit is usually the right choice: (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?

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download