Zobrazení zdroje stránky Python - lineární filter

The filter function is implemented as a direct II transposed structure. This means that the filter implements

  a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[nb]*x[n-nb]
                        - a[1]*y[n-1] - ... - a[na]*y[n-na]

using the following difference equations:

 y[m] = b[0]*x[m] + z[0,m-1]
 z[0,m] = b[1]*x[m] + z[1,m-1] - a[1]*y[m]
 ...
 z[n-3,m] = b[n-2]*x[m] + z[n-2,m-1] - a[n-2]*y[m]
 z[n-2,m] = b[n-1]*x[m] - a[n-1]*y[m]

where m is the output sample number and n=max(len(a),len(b)) is the model order.

The rational transfer function describing this filter in the z-transform domain is:

                    -1               -nb
        b[0] + b[1]z  + ... + b[nb] z
 Y(z) = ---------------------------------- X(z)
                    -1               -na
        a[0] + a[1]z  + ... + a[na] z

více na docs.scipy.org [http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html]


 #!/bin/python
 # coding:utf-8
 
 from scipy.signal import *
 
 def filter(b, a, x, zi=None):
 	y = zeros(len(x))
 	zf = zeros(len(zi))
 	
 	n = max(len(a),len(b))-1
 	_z = zi	
 	
 	for m in range(len(x)):
 		y[m] = b[0]/a[0] * x[m] + _z[0]
 
 		for i in range(len(zi)):
 			zf[i] = 0.
 			if i == len(zi)-1:
 				zf[i] = zf[i] - a[i+1]/a[0]*y[m]
 			else:
 				zf[i] = zf[i] + _z[i+1] - a[i+1]/a[0]*y[m]
 		_z = zf
  		
 	return [y,zf] 

Ještě nějaké testovací data

 #testování
 b = [4.]
 a = [1.,2.,1.,2.]
 x = [1.,2.,2.,0.,0.,1.,0.,1.]
 zi = [1,2.,0.5]
 
 print filter(b,a,x,zi)
 print lfilter(b,a,x,zi=zi)