Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Matlab filter not compatible with Python lfilter

Good day,

I have been fiddling around porting matlab code to python and I ran into this weird issue. I googled around a bit but found no information that indicates I am doing something wrong.

The core of the issue is Matlab's filter(b, a, data) (which is built-in into matlab) is generating a different output when comparing to Python's scipy.signal.lfilter

This is the issue, performed on an arbitrary pink noise signal

I have filter coefficients given to me by a third party and they are as follows:

a0 = 1
a1 = -1.69065929318241
a2 = 0.73248077421585
b0 = 1.53512485958697
b1 = -2.69169618940638
b2 = 1.19839281085285

In matlab I initialize the numerator/denominator as follows:

a = [a0 a1 a2];
b = [b0 b1 b2];

In python I do it like this:

a = np.array([a0, a1, a2])
b = np.array([b0, b1, b2])

After reading in the signal in both matlab/python I print out the first 15 samples to make sure that you guys know the input is the same

Matlab:

   0.061920166015625
  -0.050170898437500
  -0.117370605468750
  -0.065979003906250
  -0.013854980468750
  -0.042663574218750
   0.107452392578125
  -0.044006347656250
   0.115112304687500
  -0.043457031250000
  -0.028778076171875
  -0.128234863281250
   0.045227050781250
  -0.091796875000000
   0.315063476562500

Python:

[[ 0.06192017]
 [-0.0501709 ]
 [-0.11737061]
 [-0.065979  ]
 [-0.01385498]
 [-0.04266357]
 [ 0.10745239]
 [-0.04400635]
 [ 0.1151123 ]
 [-0.04345703]
 [-0.02877808]
 [-0.12823486]
 [ 0.04522705]
 [-0.09179688]]

then I call the filter functions

Matlab:

out = filter(b,a,data);
out(1:15)

ans =

   0.095055186160338
  -0.082982934483728
  -0.180851002009017
  -0.090458464750796
  -0.004794343459254
  -0.049115794227541
   0.183660200687651
  -0.061428954478571
   0.185550654888710
  -0.070597744360580
  -0.044524076275862
  -0.195036835228527
   0.082983215098531
  -0.133175807494538
   0.499012320158226

Python:

out = lfilter(b,a,data)
print out[0:14]

[[ 0.09505519]
 [-0.07701859]
 [-0.18017853]
 [-0.10128601]
 [-0.02126912]
 [-0.06549391]
 [ 0.16495284]
 [-0.06755524]
 [ 0.17671176]
 [-0.06671197]
 [-0.04417794]
 [-0.19685653]
 [ 0.06942917]
 [-0.14091966]]

Extra info:

Matlab R2012a

2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] -> python

1.6.2 -> numpy

My question is this: Am I doing something wrong or did I just find a bug in an essential and basic function in the scipy package?

King regards,

K

EDIT: below in comments it was suggest to feed it with an impulse (I kept the coeffs) Matlab:

   1.535124676585826
  -0.096323067867721
  -0.088906133468550
  -0.079755185442926
  -0.069716811972987
  -0.059448236072219
  -0.049440488368964
  -0.040042331136521
  -0.031483732058538
  -0.023898026476545
  -0.017342192117849
  -0.011814893332425
  -0.007272136901341
  -0.003640523618135
  -0.000828184619352

Python:

[[ 1.53512468]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]]

There is definitely something going wrong here, this doesn't look like a filter at all...

like image 899
kbroos Avatar asked Jun 05 '13 09:06

kbroos


1 Answers

This is not a bug. Matlab's filter operates on the first dimension of the array, while scipy.signal.lfilter by default operates on the the last dimension.

From your question I see that your data array has a second dimension (perhaps empty?). When you use lfilter it defaults to axis=-1, which will give the answer you got for python. If you want the same behaviour of matlab you need to specify the first axis or squeeze the array (if the second dimension has a size of 1):

out = lfilter(b, a, data, axis=0)
out = lfilter(b, a, np.squeeze(data))

Both of these return the following:

[ 0.09505519
 -0.08298293
 -0.180851
 -0.09045846
 -0.00479434
 -0.04911579
  0.1836602
 -0.06142895
  0.18555065
 -0.07059774
 -0.04452408
 -0.19503684
  0.08298322
 -0.13317581
  0.49901232]
like image 111
tiago Avatar answered Sep 20 '22 23:09

tiago