addition of minreal for TransferFunction

repagh · repagh · commit bba73aa3b75c · 2012-11-02T16:23:54.000Z
diff --git a/scratch/minreal.py b/scratch/minreal.py
@@ -4,31 +4,34 @@
 
 s = TransferFunction([1, 0], [1])
 
-h = (s+1)*(s+2)/(s+2)/(s**2+s+1)
+h = (s+1)*(s+2)/s/(s+2)/(s**2+s+1)
 tol = None
 
 sqrt_eps = np.sqrt(float_info.epsilon)
 
-for i in range(h.inputs):
-    for j in range(h.outputs):
-        newzeros = []
-        zeros = np.roots(h.num[i][j])
-        poles = np.roots(h.den[i][j])
-        # check all zeros
-        for z in zeros:
-            t = tol or \
-                1000 * max(float_info.epsilon, abs(z) * sqrt_eps)
-            idx = np.where(abs(z - poles) < t)[0]
-            print idx
-            if len(idx):
-                print "found match %s" % (abs(z - poles) < t)
-                poles = np.delete(poles, idx[0])
-            else:
-                newzeros.append(z)
-        print newzeros
-        print poles
+i = 0
+j = 0
+newzeros = []
+zeros = np.roots(h.num[i][j])
+poles = np.roots(h.den[i][j])
+gain = h.num[i][j][0] / h.den[i][j][0]
+# check all zeros
+for z in zeros:
+    t = tol or \
+        1000 * max(float_info.epsilon, abs(z) * sqrt_eps)
+    idx = np.where(abs(z - poles) < t)[0]
+    print idx
+    if len(idx):
+        print "found match %s" % (abs(z - poles) < t)
+        poles = np.delete(poles, idx[0])
+    else:
+        newzeros.append(z)
+
+print newzeros
+print poles
             
-                
-        
+h2 = TransferFunction(gain*np.real(np.poly(newzeros)), np.real(np.poly(poles)))
+
+print h2        
 
         
diff --git a/scratch/svninfo.txt b/scratch/svninfo.txt
@@ -0,0 +1,9 @@
+revisions	work			merged?
+
+218 - 212       addition of frd
+
+213   		scratch code frd	DO NOT MERGE
+
+214		scratch code minreal	DO NOT MERGE
+
+215		minreal for xferfcn
diff --git a/src/xferfcn.py b/src/xferfcn.py
@@ -26,6 +26,7 @@
 TransferFunction.pole
 TransferFunction.zero
 TransferFunction.feedback
+TransferFunction.minreal
 TransferFunction.returnScipySignalLti
 TransferFunction._common_den
 _tfpolyToString
@@ -76,7 +77,8 @@
 
 # External function declarations
 from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
-    polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze
+    polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, where, \
+    delete, real, poly
 from scipy.signal import lti
 from copy import deepcopy
 from lti import Lti
@@ -505,6 +507,46 @@ def feedback(self, other, sign=-1):
         # But this does not work correctly because the state size will be too
         # large.
 
+    def minreal(self, tol=None):
+        """Remove cancelling pole/zero pairs from a transfer function"""
+        # based on octave minreal
+
+        # default accuracy
+        from sys import float_info
+        sqrt_eps = sqrt(float_info.epsilon)
+
+        # pre-allocate arrays
+        num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
+        den = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
+
+        for i in range(self.outputs):
+            for j in range(self.inputs):
+
+                # split up in zeros, poles and gain
+                newzeros = []
+                zeros = roots(self.num[i][j])
+                poles = roots(self.den[i][j])
+                gain = self.num[i][j][0] / self.den[i][j][0]
+                
+                # check all zeros
+                for z in zeros:
+                    t = tol or \
+                        1000 * max(float_info.epsilon, abs(z) * sqrt_eps)
+                    idx = where(abs(z - poles) < t)[0]
+                    if len(idx):
+                        # cancel this zero against one of the poles
+                        poles = delete(poles, idx[0])
+                    else:
+                        # keep this zero
+                        newzeros.append(z)
+                        
+                # keep result
+                num[i][j] = gain * real(poly(newzeros))
+                den[i][j] = real(poly(poles))
+
+        # end result
+        return TransferFunction(num, den)
+
     def returnScipySignalLti(self):
         """Return a list of a list of scipy.signal.lti objects.
         
diff --git a/tests/xferfcn_test.py b/tests/xferfcn_test.py
@@ -433,6 +433,15 @@ def testConvertToTransferFunction(self):
                 np.testing.assert_array_almost_equal(tfsys.num[i][j], num[i][j])
                 np.testing.assert_array_almost_equal(tfsys.den[i][j], den[i][j])
 
+    def testMinreal(self):
+        """Try the minreal function"""
+        s = TransferFunction([1, 0], [1])
+        h = (s+1)*(s+2.00000000001)/(s+2)/(s**2+s+1)
+        hm = h.minreal()
+        hr = (s+1)/(s**2+s+1)
+        np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0])
+        np.testing.assert_array_almost_equal(hm.den[0][0], hr.den[0][0])
+
 def suite():
     return unittest.TestLoader().loadTestsFromTestCase(TestXferFcn)
 
