Python numpy.roots函数代码示例

def verify_roots_of_generated_polynomial_on_unit_circle(num_of_lamda, lamda_lst):

    # we need to flip the coefficient left and right, because
    # numpy.roots function expect the coefficient in the following format:
    # p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]

    # numpy.fliplr require 2D array, for 1D array, we need the follow trick
    lamda_new = np.fliplr([lamda_lst])[0]

    roots = np.roots(lamda_new)

    # get the absolute of the roots to see if they are on unit circle
    roots_abs = np.absolute(roots)

    # since leeya polynomial requires one symmetry, that is p[0] = p[n], p[1] = p[n-1], etc...
    # above polynomial will not meet this symmetry, let's just create a symmetric version

    lamda_symmetric = np.concatenate([lamda_new, lamda_lst])
    roots_symmetric = np.roots(lamda_symmetric)

    # get the absolute of the roots to see if they are on unit circle
    roots_symmetric_abs = np.absolute(roots_symmetric)

    # print "num_of_lamda", num_of_lamda,  'roots_abs =', roots_abs
    print "num_of_lamda", num_of_lamda,  'roots_symmetric_abs =', roots_symmetric_abs
    # print "num_of_lamda", num_of_lamda,  'roots_symmetric =', roots_symmetric

    return roots, roots_abs, roots_symmetric, roots_symmetric_abs

def test_calculateTF():
	"""Test function for calculateTF()"""
	from ._utils import cplxpair
	ABCD = [[1.000000000000000, 0., 0., 0.044408783846879, -0.044408783846879],
	        [0.999036450096481, 0.997109907515262, -0.005777399147297, 0., 0.499759089304780],
	        [0.499759089304780, 0.999036450096481, 0.997109907515262,  0., -0.260002096136488],
	        [0,                 0,                 1.000000000000000,  0, -0.796730400347216]]
	ABCD = np.array(ABCD)
	ntf, stf = calculateTF(ABCD)
	ntf_zeros, ntf_poles = np.roots(ntf.num), np.roots(ntf.den)
	stf_zeros, stf_poles = np.roots(stf.num), np.roots(stf.den)
	mntf_poles = np.array((1.498975311463384, 1.102565142679772, 0.132677264750882))
	mntf_zeros = np.array((0.997109907515262 + 0.075972576202904j, 
	                       0.997109907515262 - 0.075972576202904j, 
	                       1.000000000000000 + 0.000000000000000j)
	                      )
	mstf_zeros = np.array((-0.999999999999996,))
	mstf_poles = np.array((1.498975311463384, 1.102565142679772, 0.132677264750882))
	# for some reason, sometimes the zeros are in different order.
	ntf_zeros, mntf_zeros = cplxpair(ntf_zeros), cplxpair(mntf_zeros)
	stf_zeros, mstf_zeros = cplxpair(stf_zeros), cplxpair(mstf_zeros)
	ntf_poles, mntf_poles = cplxpair(ntf_poles), cplxpair(mntf_poles)
	stf_poles, mstf_poles = cplxpair(stf_poles), cplxpair(mstf_poles)
	assert np.allclose(ntf_zeros, mntf_zeros, rtol=1e-5, atol=1e-8)
	assert np.allclose(ntf_poles, mntf_poles, rtol=1e-5, atol=1e-8)
	assert np.allclose(stf_zeros, mstf_zeros, rtol=1e-5, atol=1e-8)
	assert np.allclose(stf_poles, mstf_poles, rtol=1e-5, atol=1e-8)

def daub(p):
    """
    The coefficients for the FIR low-pass filter producing Daubechies wavelets.

    p>=1 gives the order of the zero at f=1/2.
    There are 2p filter coefficients.

    Parameters
    ----------
    p : int
        Order of the zero at f=1/2, can have values from 1 to 34.

    """
    sqrt = np.sqrt
    if p < 1:
        raise ValueError("p must be at least 1.")
    if p==1:
        c = 1/sqrt(2)
        return np.array([c,c])
    elif p==2:
        f = sqrt(2)/8
        c = sqrt(3)
        return f*np.array([1+c,3+c,3-c,1-c])
    elif p==3:
        tmp  = 12*sqrt(10)
        z1 = 1.5 + sqrt(15+tmp)/6 - 1j*(sqrt(15)+sqrt(tmp-15))/6
        z1c = np.conj(z1)
        f = sqrt(2)/8
        d0 = np.real((1-z1)*(1-z1c))
        a0 = np.real(z1*z1c)
        a1 = 2*np.real(z1)
        return f/d0*np.array([a0, 3*a0-a1, 3*a0-3*a1+1, a0-3*a1+3, 3-a1, 1])
    elif p<35:
        # construct polynomial and factor it
        if p<35:
            P = [comb(p-1+k,k,exact=1) for k in range(p)][::-1]
            yj = np.roots(P)
        else:  # try different polynomial --- needs work
            P = [comb(p-1+k,k,exact=1)/4.0**k for k in range(p)][::-1]
            yj = np.roots(P) / 4
        # for each root, compute two z roots, select the one with |z|>1
        # Build up final polynomial
        c = np.poly1d([1,1])**p
        q = np.poly1d([1])
        for k in range(p-1):
            yval = yj[k]
            part = 2*sqrt(yval*(yval-1))
            const = 1-2*yval
            z1 = const + part
            if (abs(z1)) < 1:
                z1 = const - part
            q = q * [1,-z1]

        q = c * np.real(q)
        # Normalize result
        q = q / np.sum(q) * sqrt(2)
        return q.c[::-1]
    else:
        raise ValueError("Polynomial factorization does not work "
              "well for p too large.")

    def splitBimodal(self, x, y, largepoly=30):
        p = np.polyfit(x, y, largepoly) # polynomial coefficients for fit

        extrema = np.roots(np.polyder(p))
        extrema = extrema[np.isreal(extrema)]
        extrema = extrema[(extrema - x[1]) * (x[-2] - extrema) > 0] # exclude the endpoints due false maxima during fitting
        try:
            root_vals = [sum([p[::-1][i]*(root**i) for i in range(len(p))]) for root in extrema]
            peaks = extrema[np.argpartition(root_vals, -2)][-2:] # find two peaks of bimodal distribution

            mid, = np.where((x - peaks[0])* (peaks[1] - x) > 0)
             # want data points between the peaks
        except:
            warnings.warn("Peak finding failed!")
            return None

        try:
            p_mid = np.polyfit(x[mid], y[mid], 2) # fit middle section to a parabola
            midpoint = np.roots(np.polyder(p_mid))[0]
        except:
            warnings.warn("Polynomial fit between peaks of distribution poorly conditioned. Falling back on using the minimum! May result in inaccurate split determination.")
            if len(mid) == 0:
                return None

            midx = np.argmin(y[mid])
            midpoint = x[mid][midx]

        return midpoint

    def setUp(self):
        ABCD = [[1.0, 0.0, 0.0, 0.044408783846879, -0.044408783846879],
                [0.999036450096481, 0.997109907515262, -0.005777399147297,
                 0.0, 0.499759089304780],
                [0.499759089304780, 0.999036450096481, 0.997109907515262,
                 0.0, -0.260002096136488],
                [0.0, 0.0, 1.0,  0.0, -0.796730400347216]]
        ABCD = np.array(ABCD)
        (ntf, stf) = ds.calculateTF(ABCD)
        (ntf_zeros, ntf_poles) = (np.roots(ntf.num), np.roots(ntf.den))
        (stf_zeros, stf_poles) = (np.roots(stf.num), np.roots(stf.den))
        mntf_poles = np.array((1.498975311463384, 1.102565142679772,
                               0.132677264750882))
        mntf_zeros = np.array((0.997109907515262 + 0.075972576202904j,
                               0.997109907515262 - 0.075972576202904j,
                               1.000000000000000 + 0.000000000000000j))
        mstf_zeros = np.array((-0.999999999999996,))
        mstf_poles = np.array((1.498975311463384, 1.102565142679772,
                               0.132677264750882))

        # for some reason, sometimes the zeros are in different order.
        (self.ntf_zeros, self.mntf_zeros) = (cplxpair(ntf_zeros),
                                             cplxpair(mntf_zeros))
        (self.stf_zeros, self.mstf_zeros) = (cplxpair(stf_zeros),
                                             cplxpair(mstf_zeros))
        (self.ntf_poles, self.mntf_poles) = (cplxpair(ntf_poles),
                                             cplxpair(mntf_poles))
        (self.stf_poles, self.mstf_poles) = (cplxpair(stf_poles),
                                             cplxpair(mstf_poles))

def real_roots(poly_coeffs, x1=None, x2=None, warn=False):

# Supress warnings (default)
     if (warn==False):
          warnings.simplefilter('ignore', np.RankWarning)

# Evaluate roots, keeping only the real parts or those without an imaginary component
     re_roots = \
         np.roots(poly_coeffs)[np.roots(poly_coeffs).imag == 0.].real

# Go through limit possibilities, returning the appropriate values
# If no limits were given then return all real roots
     if (x1==None and x2==None):
          return re_roots
# The following are cases where either or both limits are given
     elif (x2==None):  # If only lower limit was given
          return re_roots[(re_roots >= x1)]
     elif (x1==None):  # If only upper limit was given
          return re_roots[(re_roots <= x2)]
     else:             # If both limits are given
          # Check that x1 < x2 and fix if necessary
          if (x1 > x2):
               temp_x = x1
               x1 = x2
               x2 = temp_x
          return re_roots[(re_roots >= x1) & (re_roots <= x2)]

def tf2zpk(b, a):
    """Return zero, pole, gain (z,p,k) representation from a numerator,
    denominator representation of a linear filter.

    Parameters
    ----------
    b : ndarray
        Numerator polynomial.
    a : ndarray
        Denominator polynomial.

    Returns
    -------
    z : ndarray
        Zeros of the transfer function.
    p : ndarray
        Poles of the transfer function.
    k : float
        System gain.

    Notes
    -----
    If some values of `b` are too close to 0, they are removed. In that case,
    a BadCoefficients warning is emitted.

    """
    b, a = normalize(b, a)
    b = (b + 0.0) / a[0]
    a = (a + 0.0) / a[0]
    k = b[0]
    b /= b[0]
    z = roots(b)
    p = roots(a)
    return z, p, k

def isstable(b,a,ftype='digital'):
    """Determine whether IIR filter (b,a) is stable
    
    Parameters
    ----------
        b: ndarray
            filter numerator coefficients
        a: ndarray
            filter denominator coefficients
        ftype: string
            type of filter (`digital` or `analog`)
    Returns
    -------
        stable: bool
            whether filter is stable or not
                
    """

    if ftype=='digital':
        v = np.roots(a)
        if np.any(np.abs(v)>1.0):
            return False
        else:
            return True
    elif ftype=='analog':
        v = np.roots(a)
        if np.any(np.real(v)<0):
            return False
        else:
            return True

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf)

def fit_edge_hist(bins, counts, fwhm_guess=10.0):
    if len(bins) == len(counts)+1: bins = bins[:-1]+0.5*(bins[1]-bins[0]) # convert bin edge to bin centers if neccesary
    pfit = np.polyfit(bins, counts, 3)
    edgeGuess = np.roots(np.polyder(pfit,2))
    try:
        preGuessX, postGuessX = np.sort(np.roots(np.polyder(pfit,1)))
    except:
        raise ValueError("failed to generate guesses")
    use = bins>(edgeGuess+2*fwhm_guess)
    if np.sum(use)>4:
        pfit2 = np.polyfit(bins[use], counts[use],1)
        slope_guess = pfit2[0]
    else:
        slope_guess=1
    pGuess = np.array([edgeGuess, np.polyval(pfit,preGuessX), np.polyval(pfit,postGuessX),fwhm_guess,slope_guess],dtype='float64')

    try:
        pOut = curve_fit(edge_model, bins, counts, pGuess)
    except:
        return (0,0,0,0,0,0)
    (edgeCenter, preHeight, postHeight, fwhm, bgSlope) = pOut[0]
    model_counts = edge_model(bins, edgeCenter, preHeight, postHeight, fwhm, bgSlope)
    num_degree_of_freedom = float(len(bins)-1-5) # num points - 1 - number of fitted parameters
    chi2 = np.sum(((counts - model_counts)**2)/model_counts)/num_degree_of_freedom
    return (edgeCenter, preHeight, postHeight, fwhm, bgSlope, chi2)

def lsf(fir_filt):
  """
  Find the Line Spectral Frequencies (LSF) from a given FIR filter.

  Parameters
  ----------
  filt :
    A LTI FIR filter as a LinearFilter object.

  Returns
  -------
  A tuple with all LSFs in rad/sample, alternating from the forward prediction
  and backward prediction filters, starting with the lowest LSF value.

  """
  den = fir_filt.denominator
  if len(den) != 1:
    raise ValueError("Filter has feedback")
  elif den[0] != 1: # So we don't have to worry with the denominator anymore
    fir_filt /= den[0]

  from numpy import roots
  rev_filt = ZFilter(fir_filt.numerator[::-1]) * z ** -1
  P = fir_filt + rev_filt
  Q = fir_filt - rev_filt
  roots_p = roots(P.numerator[::-1])
  roots_q = roots(Q.numerator[::-1])
  lsf_p = sorted(phase(roots_p))
  lsf_q = sorted(phase(roots_q))
  return reduce(operator.concat, xzip(*sorted([lsf_p, lsf_q])), tuple())

def polyxval(poly_coeffs, y, x1=None, x2=None, warn=False):
     if (x1==None or x2==None):
          print "Must assign range [x1, x2] in which to search for x values."
          return None

     if(x1==x2):
          print "x1 must not equal x2."
          return None

     if (x1 > x2):
          temp_x = x1
          x1 = x2
          x2 = temp_x

# Supress warnings (default)
     if (warn==False):
          warnings.simplefilter('ignore', np.RankWarning)

     poly_coeffs_y = poly_coeffs
     # subtract y-value from zeroth order coefficient (i.e. no x's)
     poly_coeffs_y[len(poly_coeffs_y)-1] -= y
     re_roots = \
         np.roots(poly_coeffs_y)[np.roots(poly_coeffs_y).imag == 0.].real

# restrict solution to range [x1, x2]
     x_val = re_roots[(re_roots >= x1) & (re_roots <= x2)]

     return x_val

def plot_range(num, den):

    # The corner frequencies
    zero = sort(abs(roots(num)))    
    pole = sort(abs(roots(den)))
    
    # Calculate the minimum and maximum corner frequencies needed
    if len(pole) == 0:
        corner_min = zero[0]
        corner_max = zero[-1]
    
    elif len(zero) == 0:
        corner_min = pole[0]
        corner_max = pole[-1]
        
    elif len(zero) > 0 and len(pole) > 0:
        corner_min = min(zero[0], pole[0])
        corner_max = max(zero[-1], pole[-1]) 
    
    else:
        corner_min, corner_max = 0.1, 10
    
    # start from 2 decades lower than the lowest corner 
    # and end at 2 decades above the highest corner
    freq_range = [10 ** (floor(log10(corner_min)) - 1), 
                   10 ** (floor(log10(corner_max)) + 2)]
    
    return freq_range

def ARLineSpectra(ar):
    """
    Convert AR coeffs to LSPs

    From wikipedia:
    A palindromic polynomial (i.e., P) of odd degree has -1 as a root.
    An antipalindromic polynomial (i.e., Q) has 1 as a root.
    An antipalindromic polynomial of even degree has -1 and 1 as roots
    """
    order = ar.shape[-1]
    ret = np.zeros(ar.shape)
    for a, o in core.refiter([ar, ret], core.newshape(ar.shape)):
        p = np.ones((order+2))
        q = np.ones((order+2))
        q[-1] = -1.0
        for i in range(order):
            p[i+1] = -a[i] - a[order-i-1]
            q[i+1] = -a[i] + a[order-i-1]
        pr = np.roots(p)
        qr = np.roots(q)

        j = 0
        an = np.ndarray((order+2))
        for i in range(len(pr)):
            if np.imag(pr[i]) >= 0.0:
                an[j] = np.angle(pr[i])
                j += 1
            if np.imag(qr[i]) >= 0.0:
                an[j] = np.angle(qr[i])
                j += 1
        # The angle list (an) will always contain both 0 and pi; they
        # will move to the ends after the sort
        o[...] = np.sort(an)[1:-1]
    return ret;

def asymptote(num, den):

    # create a Python list for the zeros and the poles of the system
    zero = list(sort(abs(roots(num))))    
    pole = list(sort(abs(roots(den))))
    
    #calculate the low frequency gain -- type 0 system
    lf_gain = 20 * log10(abs((num[-1] / den[-1])))
    
    # create an empty matrix to contain the corner frequencies and
    # the corresponding slope indicator (+1 or -1)
    corners = zeros((len(zero) + len(pole) + 2, 2))
    
    starting_freq, end_freq = plot_range(num, den)
    corners[0] = [starting_freq, 0]
    corners[-1] = [end_freq, 0]
    
    # take the first elements from the list of poles and zeros
    # compare them and assign the slope indicator
    # delete the corresponding valuefrom the original list of poles and zeros
    for count in range(len(zero) + len(pole)): 
            
        if len(zero) > 0:
            a = zero[0]
        else:
            a = inf
            
        if len(pole) > 0:
            b = pole[0]
        else:
            b = inf
        
        c = min(a, b)
                
        if c == a:
            corners[count + 1] = [c, 1]
            if len(zero) > 0:
                zero.pop(0)

        if c == b:
            corners[count + 1] = [c, -1]
            if len(pole) > 0:
                pole.pop(0)
         
    # now calculate the gains at the corners using 
    # gain = +/- 20log10(upper_corner / lower_corner) 
    asymptotic_gain = zeros_like(corners)   
    asymptotic_gain[0, 1] = lf_gain
    asymptotic_gain[1, 1] = lf_gain
    
    gain = lf_gain
    multiplier = cumsum(corners[:, 1])
    for k in range(2, len(corners)):    
        gain += multiplier[k-1] * 20 * log10(corners[k, 0] / corners[k-1, 0])
        asymptotic_gain[k, 1] = gain
        
    asymptotic_gain[:, 0] = corners[:, 0]
                   
    return asymptotic_gain