Added fminbound and some filter_design.

Lib/optimize/optimize.py

@@ -5,17 +5,20 @@
 # guarantee implied provided you keep this notice in all copies.
 # *****END NOTICE************
 
-# A collection of optimization algorithms.  Version 0.3.1
+# A collection of optimization algorithms.  Version 0.4.1
+# CHANGES
+#  Added fminbound (July 2001)
 
 # Minimization routines
 """optimize.py
 
 A collection of general-purpose optimization routines using Numeric
 
-fmin        ---      Nelder-Mead Simplex algorithm (uses only function calls)
-fminBFGS    ---      Quasi-Newton method (uses function and gradient)
+fmin        ---      Nelder-Mead Simplex algorithm (uses only function calls).
+fminBFGS    ---      Quasi-Newton method (uses function and gradient).
 fminNCG     ---      Line-search Newton Conjugate Gradient (uses function, 
-                     gradient and hessian (if it's provided))
+                     gradient and hessian (if it's provided)).
+fminbound   ---      Bounded minimization for scalar functions.
 
 """
 
@@ -25,7 +28,8 @@
 max = MLab.max
 min = MLab.min
 abs = Num.absolute
-__version__="0.3.1"
+sqrt = Num.sqrt
+__version__="0.4.1"
 
 def rosen(x):  # The Rosenbrock function
     return MLab.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
@@ -602,8 +606,166 @@ def fminNCG(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
     else:        
         return xk
 
+
+
+def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500, 
+              full_output=0, disp=1):
+    """Bounded minimization for scalar functions.
+
+    Description:
+
+      Finds a local minimizer of the scalar function func in the interval
+      x1 < xopt < x2.
+
+    Inputs:
+
+      func -- the function to be minimized (must accept scalar input and return
+              scalar output).
+      x1, x2 -- the optimization bounds.
+      args -- extra arguments to pass to function.
+      xtol -- the convergence tolerance.
+      maxfun -- maximum function evaluations.
+      full_output -- Non-zero to return optional outputs. 
+      disp -- Non-zero to print messages.
+              0 : no message printing.
+              1 : non-convergence notification messages only.
+              2 : print a message on convergence too.
+              3 : print iteration results. 
+
+
+    Outputs: (xopt, {fval, ierr, numfunc})
+
+      xopt -- The minimizer of the function over the interval.
+      fval -- The function value at the minimum point.
+      ierr -- An error flag (1 if converged, 0 if maximum number of
+              function calls reached).
+      numfunc -- The number of function calls.
+    """
+
+    if x1 > x2:
+        raise ValueError, "The lower bound exceeds the upper bound."
+
+    flag = 1
+    header = ' Func-count     x          f(x)          Procedure'
+    step='       initial'
+
+    seps = sqrt(2.2e-16)
+    c = 0.5*(3.0-sqrt(5.0))
+    a, b = x1, x2
+    v = a + c*(b-a)
+    w, xf = v, v
+    d = e = 0.0
+    x = xf
+    fx = func(x,*args)
+    num = 1
+    fmin_data = (1, xf, fx)
+    if disp > 2:
+        print (" ")
+        print (header)
+        print "%5.0f   %12.6g %12.6g %s" % (fmin_data + (step,))
+
+    fv = fw = fx
+    xm = 0.5*(a+b)
+    tol1 = seps*abs(xf) + xtol / 3.0
+    tol2 = 2.0*tol1
+
+    maxcount = maxfun
+    # Main loop
+
+    while ( abs(xf-xm) > (tol2 - 0.5*(b-a)) ):
+        gs = 1
+        # Check for parabolic fit
+        if abs(e) > tol1:
+            gs = 0
+            r = (xf-w)*(fx-fv)
+            q = (xf-v)*(fx-fw)
+            p = (xf-v)*q - (xf-w)*r
+            q = 2.0*(q-r)
+            if q > 0.0: p = -p
+            q = abs(q)
+            r = e
+            e = d
+
+            # Check for acceptability of parabola
+            if ( (abs(p) < abs(0.5*q*r)) and (p > q*(a-xf)) and (p < q*(b-xf))):
+                d = (p+0.0) / q;
+                x = xf + d
+                step = '       parabolic'
+
+                if ((x-a) < tol2) or ((b-x) < tol2):
+                    si = Numeric.sign(xm-xf) + ((xm-xf)==0)
+                    d = tol1*si
+            else:      # do a golden section step
+                gs = 1
+
+        if gs:  # Do a golden-section step
+            if xf >= xm:
+                e=a-xf
+            else:
+                e=b-xf
+            d = c*e
+            step = '       golden'
+
+        si = Numeric.sign(d) + (d == 0)
+        x = xf + si*max([abs(d), tol1])
+        fu = func(x,*args)
+        num += 1
+        fmin_data = (num, x, fu)
+        if disp > 2:
+            print "%5.0f   %12.6g %12.6g %s" % (fmin_data + (step,))
+
+        if fu <= fx:
+            if x >= xf:
+                a = xf
+            else:
+                b = xf
+            v, fv = w, fw
+            w, fw = xf, fx
+            xf, fx = x, fu
+        else:
+            if x < xf:
+                a = x
+            else:
+                b = x
+            if (fu <= fw) or (w == xf):
+                v, fv = w, fw
+                w, fw = x, fu
+            elif (fu <= fv) or (v == xf) or (v == w):
+                v, fv = x, fu
+
+        xm = 0.5*(a+b)
+        tol1 = seps*abs(xf) + xtol/3.0
+        tol2 = 2.0*tol1
+
+        if num >= maxcount:
+            flag = 0
+            fval = fx
+            if disp > 0:
+                _endprint(x, flag, fval, maxfun, tol, disp)
+            if full_output:
+                return xf, fval, flag, num
+            else:
+                return xf
+
+    fval = fx
+    if disp > 0:
+        _endprint(x, flag, fval, maxfun, xtol, disp)
 
+    if full_output:
+        return xf, fval, flag, num
+    else:
+        return xf
 
+
+def _endprint(x, flag, fval, maxfun, xtol, disp):
+    if flag == 1:
+        if disp > 1:
+            print "\nOptimization terminated successfully;\n the returned value" + \
+                  " satisfies the termination criteria\n (using xtol = ", xtol, ")"
+    if flag == 0:
+        print "\nMaximum number of function evaluations exceeded --- increase maxfun argument.\n"
+    return
+
 if __name__ == "__main__":
     import string
     import time

Lib/signal/filter_design.py

@@ -0,0 +1,180 @@
+import Numeric
+Num = Numeric
+pi = Numeric.pi
+import scipy.special
+
+def butter(N, Wn, bandtype='band', analog=0, output='ba'):
+    """Butterworth digital and analog filter design.
+
+    Description:
+
+      Design an Nth order lowpass digital Butterworth filter and return the
+      filter coeeficients in (B,A), zero-pole, or state-space form.
+
+    Inputs:
+
+      
+    """
+    pass
+
+def cheby1():
+    pass
+
+def cheby2():
+    pass
+
+def ellip():
+    pass
+
+def maxflat():
+    pass
+
+def yulewalk():
+    pass
+
+
+def band_stop_cost(wp, ind, WP, WS, rs, rp, type):
+    """Band Stop Objective Function for order minimization
+
+    Description:
+
+      Returns the non-integer order for an analog band stop filter.
+
+    Inputs:
+
+      wp -- passband edge
+      ind -- index specifying which passband edge to vary (0 or 1).
+      WP -- two element vector of fixed passband edges.
+      WS -- two element vector of fixed stopband edges.
+      rs -- amount in dB of attenuation in stopband.
+      rp -- amount in dB of ripple in the passband.
+      type -- 'butter', 'cheby', or 'ellip':
+
+    Outputs: (n,)
+
+      n -- filter order (non-integer)
+    """
+
+    WP[ind] = wp
+    WSflip = array([WS(0), -WS(1)])
+    WA = WS*(WP(0)-WP(1)) / (WS**2 - WP(0)*WP(1))
+    WA = min(abs(WA))
+
+    if type == 'butter':
+        n = (Num.log10 ( (10**(0.1*abs(rs)) - 1.0) / \
+                         (10**(0.1*abs(rp))-1)) / (2*Num.log10(WA)))
+    elif type == 'cheby':
+        n = Num.acosh(Num.sqrt((10**(0.1*abs(rs))-1) / \
+                               (10**(0.1*abs(rp))-1))) / acosh(WA)
+    elif type == 'ellip':
+        epsilon = sqrt(10**(0.1*rp)-1)
+        k1 = epsilon / sqrt(10**(0.1*rs)-1)
+        k = 1.0 / WA
+        capk = scipy.special.ellpk([k**2, 1-k**2])
+        capk1 = scipy.special.ellpk([k1**2, 1-k1**2])
+        n = (capk[0]*capk1[1] / (capk[1]*capk1[0]))
+    else:
+        raise ValueError, "Incorrect type: ", type
+
+def buttord(wp, ws, rp, rs, analog=0):
+    """Butterworth filter order selection.
+
+    """
+
+    filter_type = 2*(len(wp)-1)
+    if wp(0) < ws(0):
+        filter_type += 1
+    else:
+        filter_type += 2
+
+    # Pre-warp frequencies
+    if not analog:
+        WP = tan(pi*wp/2)
+        WS = tan(pi*ws/2)
+    else:
+        WP = wp
+        WS = ws
+
+    if ftype == 1:           # low
+        WA = WS / WP
+    elif ftype == 2:          # high
+        WA = WP / WS
+    elif ftype == 3:     # stop
+        wp0 = scipy.optimize.fminbound(band_stop_cost, WP[0], WS[1]-1e-12,
+                                       args=(0,WP,WS,rs,rp,'butter'), disp=0)
+        WP[0] = wp0
+        wp1 = scipy.optimize.fminbound(band_stop_cost, WS[1]+1e-12, WP[1],
+                                       args=(1,WP,WS,rs,rp,'butter'), disp=0)
+        WP[1] = wp1
+        WA = (WS * (WP[0]-WP[1])) / (WS**2 - WP[0]*WP[1])
+    elif ftype == 4:  # pass
+        WA = (WS**2 - WP[0]*WP[1]) / (WS* (WP[0]-WP[1]))
+
+    WA = min(abs(WA))
+
+    order = Num.ceil( Num.log10( (10**(0.1*abs(rs)) - 1) / \
+                                 (10**(0.1*abs(rp)) - 1) ) / \
+                      (2*Num.log10(WA)))
+
+    # Find the butterworth natural frequency Wo (or the "3dB" frequency")
+    # to give exactly rs at WA. W0 will be between 1 and WA
+    W0 = WA / ( ( 10**(0.1*abs(rs))-1)**(1.0/(2*order)))
+
+    # now convert this frequency back from lowpass prototype
+    # to the original analog filter
+
+    if ftype == 1:  # low
+        WN = WO*WP
+    elif ftype == 2: # high
+        WN = WP / W0
+    elif ftype == 3:  # stop
+        WN = zeros(2,Numeric.Float)
+        WN[0] = ((WP[1] - WP[0]) + sqrt((WP[1] - WP[0])**2 + \
+                                        4*W0**2 * WP[0] * WP[1])) / (2*W0)
+        WN[1] = ((WP[1] - WP[0]) - sqrt((WP[1] - WP[0])**2 + \
+                                        4*W0**2 * WP[0] * WP[1])) / (2*W0)
+        WN = Num.sort(Num.abs(WN))
+    elif ftype == 4: # pass
+        W0 = array([-W0, W0],Numeric.Float)
+        WN = -W0 * (WP[1]-WP[0]) / 2.0 + sqrt(W0**2 / 4.0 * (WP[1]-WP[0])**2 + \
+                                              WP[0]*WP[1])
+        WN = Num.sort(Num.abs(WN))
+    else:
+        raise ValueError, "Bad type."
+
+    if not analog:
+        wn = (2.0/pi)*Num.atan(WN)
+    else:
+        wn = WN
+
+    return order, wn
+
+def cheb1ord():
+    pass
+
+def cheb2ord():
+    pass
+
+def ellipord():
+    pass
+
+def besselap():
+    pass
+
+def buttap():
+    pass
+
+def cheb1ap():
+    pass
+
+def cheb2ap():
+    pass
+
+def ellipap():
+    pass
+
+def besself():
+    pass
+
+
+