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 [1]

#!/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)