[SciPy-User] scipy.signal.butter function different from Matlab

[gwenael....@aliceadsl.fr]

Hi,

I'm trying to rewrite a Matlab code in Python and I encounter some difficulties
using the butter function.

In Matlab, my code is:

% Parameters:
N=5;
F_inf=88.38834764831843;
F_sup=111.3623397675424;
Fs=32768.00013421773;

% Filter coefficients computing:
[z,p,k]=butter(N,[F_inf F_sup]/(Fs/2));

% Result:
z=[1;1;1;1;1;-1;-1;-1;-1;-1;]
p=[0.999020109086358+0.021203989980732i;...
0.999020109086358-0.021203989980732i;...
0.99789313059316+0.020239448648803i;...
0.99789313059316-0.020239448648803i;...
0.99762168426512+0.018853164086747i;...
0.99762168426512-0.018853164086747i;...
0.998184731408311+0.017659802911797i;...
0.998184731408311-0.017659802911797i;...
0.999249126282689+0.017020854458606i;...
0.999249126282689-0.017020854458606i;]
k=5.147424357763888e-014

In Python, the code is quite similar:
import numpy as np
import scipy.signal as sig

# Parameters:
N=5
F_inf=88.38834764831843
F_sup=111.3623397675424
Fs=32768.00013421773

# Filter coefficients computing:
z,p,k=sig.butter(N,np.array([freqmin,freqmax])/(Fs/2),output='zpk')

# Result:
Error message:
C:\Python26\lib\site-packages\scipy\signal\filter_design.py:221:
BadCoefficients: Badly conditionned filter coefficients (numerator): the results
may be meaningless
"results may be meaningless", BadCoefficients)

z=array([], dtype=float64)
p=array([
0.99464737+0.01603399j,0.99464737-0.01603399j,0.98633302+0.00982676j,0.98633302-0.00982676j,0.98319371+0.j])
k=4.2522321460489923e-11

Does anyone could help me to understand why the results are not the same?

I have also another question. In Matlab, after having computed the filter
coefficients, I apply the filtering to my signal X using the following code:

% To avoid round-off errors, do not use the transfer function. Instead
% get the zpk representation and convert it to second-order sections.
[sos_var,g]=zp2sos(z, p, k);
Hd=dfilt.df2sos(sos_var, g);

% Signal filtering
y = filter(Hd,X);

Is there any Python functions that could do the same?

Thanks

Gwenaël

_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[braingateway]

gwenael....@aliceadsl.fr :

No, there is only naive implementation of filter(). As I know, there is
no SOS convertion exists. The Scipy.signal is only for very very simple
signal processing. In this state, you either implement your own SOS
system convertion code and then use concatenated Scipy.lfilter() call to
achieve your go, or you stick to matlab.

LittleBigBrain

[braingateway]

braingateway :

To make it more clear, if you have 3 SO block, then you need to call 3
times of lfilter() *gain

[Ryan May]

On Thu, Nov 18, 2010 at 4:55 AM, <gwenael....@aliceadsl.fr> wrote:
> % Parameters:
> N=5;
> F_inf=88.38834764831843;
> F_sup=111.3623397675424;
> Fs=32768.00013421773;
>
> % Filter coefficients computing:
> [z,p,k]=butter(N,[F_inf F_sup]/(Fs/2));
>
> % Result:
> z=[1;1;1;1;1;-1;-1;-1;-1;-1;]
> p=[0.999020109086358+0.021203989980732i;...
>   0.999020109086358-0.021203989980732i;...
>   0.99789313059316+0.020239448648803i;...
>   0.99789313059316-0.020239448648803i;...
>   0.99762168426512+0.018853164086747i;...
>   0.99762168426512-0.018853164086747i;...
>   0.998184731408311+0.017659802911797i;...
>   0.998184731408311-0.017659802911797i;...
>   0.999249126282689+0.017020854458606i;...
>   0.999249126282689-0.017020854458606i;]
> k=5.147424357763888e-014

Just because MATLAB is giving you an answer doesn't mean it's not
having the same problems. You're asking it to make a bandpass filter
with a normalized width of 0.0014 using only 5 coefficients--that's
too tight with too few coefficients and will *always* result in
problems. That k value there is supposed to be the gain, which is
almost 0 (140dB of loss!). If you look at a plot of the transfer
function you're generating in MATLAB, you'll see it's complete
garbage. Using sig.buttord() and moving the stopband 0.0001 off of the
passbands edges, 5db attenuation in the passband and 30db in the stop,
I get an order of 26, which is >>5.

The code for that is:

Wn = np.array([F_inf F_sup])/(Fs/2)
sig.buttord(Wn, np.array([Wn[0]-0.0001, Wn[1]+0.0001]), 5, 30)

If we make the passband a bit less crazy:

import matplotlib.pyplot as plt
z,p,k = sig.butter(5, [0.3, 0.4], btype='bandpass', output='zpk')
b,a = sig.zpk2tf(z, p, k)
w,h = sig.freqz(b, a)
plt.plot(w, 10.0 * np.log10(np.abs(h)))
plt.grid()
plt.ylim(-40, 3)

You get a nice plot of the filter response, which is as expected for
only order 5.

As for the other things, it appears that there is no equivalent for SOS.

Ryan

--
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma

[gwenael....@aliceadsl.fr]

Thanks Ryan ,

Appreciate your reply.

Best regards

Gwenaël

Selon Ryan May <rma...@gmail.com>: