Evaluate filter matlab function using python-CodePudding

Based on the description of both functions, filter(b, a, x) matlab function and lfilter(b, a, x, axis=-1, z_i=None) scipy.signal should give the same result.

I take this example where the results are completely different:

import numpy as np
from scipy.signal import lfilter

bb = np.arange(0, 100, 1)

aa = 1

xx = np.tile(0.9, (100, 1))

yy = lfilter(b=bb, a=aa, x=xx)

print(yy[:10])

array([[0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.]])
print(yy.shape)
(100, 1)

# Same example on Matlab using filter
bb = 0:1:99
aa  =  1
xx =repmat(0.9, 100, 1)
dd = filter(bb, aa, xx)
dd(1:10)
ans =

         0
    0.9000
    2.7000
    5.4000
    9.0000
   13.5000
   18.9000
   25.2000
   32.4000
   40.5000
print(size(dd)) # (100 ,1)

CodePudding user response：

MATLAB and NumPy handle array shapes differently. NumPy has general n-dimensional arrays. scipy.signal.lfilter accepts an n-dimensional array, and applies the filter along each 1-d slice of the input array. Which slice is determined by the axis parameter. By default, lfilter operates on the last axis (axis=-1). You gave lfilter an array with shape (100, 1). Applying lfilter with axis=-1 to that input applies the filter 100 times, once for each row of length 1--certainly not what you want! Instead, you want to apply the filter along axis=0 (which in this 2-d case, means apply lfilter along the columns). If you change the call of lfilter to

yy = lfilter(b=bb, a=aa, x=xx, axis=0)

the values returned will match those of the MATLAB code.

In the long term, I recommend not limiting yourself to 2-d arrays. In this case, it makes more sense to create xx as a 1-d array (e.g. xx = np.full(100, fill_value=0.9)). Then lfilter will operate on the given 1-d array the way you expect, without having to specify the axis.

CodePudding user response：

The answer of @warren is good. Here are more details about the use of functions.

filter(b, a, x, zi, dim) return zf value even if zi is not given as input.

In the example above:

[y, zf] = lfilter(b=bb, a=aa, x=xx, axis=0)

zf(1:10)

ans =

   1.0e 03 *

    4.4550
    4.4541
    4.4523
    4.4496
    4.4460
    4.4415
    4.4361
    4.4298
    4.4226
    4.4145

However in lfilter(a, b, x, zi=None, axis=-1), to get the zf, you have to provide zi. For that, you can use, lfilter_zi(b, x) that construct initial conditions for lfilter for step response steady-state.

z_i = lfilter_zi(b=bb, x=xx)
y, zf = filter(b=bb, x=xx, a=aa, zi=zi, axis=0)

zf[1:10]

[[4455. ]
 [4454.1]
 [4452.3]
 [4449.6]
 [4446. ]
 [4441.5]
 [4436.1]
 [4429.8]
 [4422.6]
 [4414.5]]