[Scipy-svn] r5195 - trunk/scipy/signal

scipy-svn@scip... scipy-svn@scip...
Wed Nov 26 12:24:04 CST 2008


Author: oliphant
Date: 2008-11-26 12:24:03 -0600 (Wed, 26 Nov 2008)
New Revision: 5195

Modified:
   trunk/scipy/signal/signaltools.py
Log:
Add lfilter_zi so that filtfilt will work.

Modified: trunk/scipy/signal/signaltools.py
===================================================================
--- trunk/scipy/signal/signaltools.py	2008-11-26 00:35:21 UTC (rev 5194)
+++ trunk/scipy/signal/signaltools.py	2008-11-26 18:24:03 UTC (rev 5195)
@@ -12,7 +12,7 @@
      ravel, size, less_equal, sum, r_, iscomplexobj, take, \
      argsort, allclose, expand_dims, unique, prod, sort, reshape, \
      transpose, dot, any, mean, cosh, arccosh, \
-     arccos, concatenate
+     arccos, concatenate, flipud
 import numpy as np
 from scipy.misc import factorial
 
@@ -1509,6 +1509,33 @@
         ret = transpose(ret,tuple(olddims))
         return ret
 
+def lfilter_zi(b,a):
+    #compute the zi state from the filter parameters. see [Gust96].
+
+    #Based on:
+    # [Gust96] Fredrik Gustafsson, Determining the initial states in
+    #          forward-backward filtering, IEEE Transactions on
+    #          Signal Processing, pp. 988--992, April 1996, 
+    #          Volume 44, Issue 4
+
+    n=max(len(a),len(b))
+    
+    zin = (np.eye(n-1) - np.hstack((-a[1:n,newaxis],
+                                    np.vstack((np.eye(n-2),zeros(n-2))))))
+
+    zid=  b[1:n] - a[1:n]*b[0]
+
+    zi_matrix=linalg.inv(zin)*(np.matrix(zid).transpose())
+    zi_return=[]
+
+    #convert the result into a regular array (not a matrix)
+    for i in range(len(zi_matrix)):
+      zi_return.append(float(zi_matrix[i][0]))
+
+    return array(zi_return)
+   
+
+
 def filtfilt(b,a,x):
     # FIXME:  For now only accepting 1d arrays
     ntaps=max(len(a),len(b))



More information about the Scipy-svn mailing list