def compute_filter(self):
		theta = 2 * np.pi * self.frequency / self.nyquist * 2
		zero = np.exp(np.array([1j, -1j]) * theta)
		pole = self.module_pole * zero

		self._filter_b = np.poly(zero)
		self._filter_a = np.poly(pole)

def find_coeff(poles, mode="f"):
    """
    Find LP coefficients from a set of LP roots (poles).

    Parameters
    ----------
    poles : ndarray
        Array of LP roots (poles)
    mode : {'f', 'b'}
        Mode in which LP coefficients should be returned.  'f' for coefficients
        ordered m, m - 1,..., 1. 'b' for coefficients ordered 1, 2, ...., m.

    Returns
    -------
    c : ndarray
        LP coefficients ordered according to `mode`.

    """
    # STABILITY
    # the algorithm used here is numpy poly function which convolves
    # the roots to find the coefficients, the accuracy of this method
    # depends on the dtype of the poles parameter.
    if mode not in ['f', 'b']:
        raise ValueError("mode must be 'f'or 'b'")

    if mode == 'f':  # reverse resulting coefficients
        return np.squeeze(-np.poly(poles)[:0:-1])
    else:   # keep coefficients as is
        return np.squeeze(-np.poly(poles)[1:])

    def apply_to(self, array, samplerate):
        # Nyquist frequency
        nyquist = samplerate / 2.

        # ratio of notch freq. to Nyquist freq.
        freq_ratio = self.frequency / nyquist

        # width of the notch
        notchWidth = 0.01

        # Compute zeros
        zeros = np.array([np.exp(1j * np.pi * freq_ratio),
                         np.exp(-1j * np.pi * freq_ratio)])

        # Compute poles
        poles = (1 - notchWidth) * zeros

        b = np.poly(zeros)    # Get moving average filter coefficients
        a = np.poly(poles)    # Get autoregressive filter coefficients

        try:
            if array.ndim == 1:
                filtered = scipy.signal.filtfilt(b, a, array, padtype="even")
            else:
                filtered = np.empty(array.shape)
                for i in range(array.shape[1]):
                    filtered[:, i] = scipy.signal.filtfilt(
                        b, a, array[:, i], padtype="even")

            return filtered * self.gain + array * (1 - self.gain)
        except ValueError:
            print "Could not apply filter"
            return array

 def myzpk2tf(self, z, p, k):
         z = np.atleast_1d(z)
         k = np.atleast_1d(k)
         if len(z.shape) > 1:
                 temp = np.poly(z[0])
                 b = np.zeros((z.shape[0], z.shape[1] + 1), temp.dtype.char)
                 if len(k) == 1:
                         k = [k[0]] * z.shape[0]
                 for i in range(z.shape[0]):
                         b[i] = k[i] * poly(z[i])
         else:
                 b = k * np.poly(z)
         a = np.atleast_1d(np.poly(p))
         # Use real output if possible. Copied from numpy.poly, since
         # we can't depend on a specific version of numpy.
         if issubclass(b.dtype.type, np.complexfloating):
                 # if complex roots are all complex conjugates, the roots are real.
                 roots = np.asarray(z, complex)
                 pos_roots = np.compress(roots.imag > 0, roots)
                 neg_roots = np.conjugate(np.compress(roots.imag < 0, roots))
                 if len(pos_roots) == len(neg_roots):
                         if np.all(np.sort_complex(neg_roots) == np.sort_complex(pos_roots)):
                                 b = b.real.copy()
         if issubclass(a.dtype.type, np.complexfloating):
                 # if complex roots are all complex conjugates, the roots are real.
                 roots = np.asarray(p, complex)
                 pos_roots = np.compress(roots.imag > 0, roots)
                 neg_roots = np.conjugate(np.compress(roots.imag < 0, roots))
                 if len(pos_roots) == len(neg_roots):
                         if np.all(np.sort_complex(neg_roots) == np.sort_complex(pos_roots)):
                                 a = a.real.copy()
         return b, a

def do_plot(ax, z):
    
    tt=np.linspace(0,1,100)
    yy = z[0]*(1-tt)+z[1]*(tt)
    z1=z[:5]
    z2=z[5:]
    
    P = np.poly(z1) 
    Q = np.poly(z2)
    Z0 = np.roots((P-Q)[1:])
    
    def Z(t):
        F=P*(1-t)+Q*t
        return np.roots(F)
        
    N = 100
    
    ZZ = np.zeros([5, N], np.complex)
    for i in range(N):
        ZZ[:,i]=Z(i/50.0-0.5).transpose()

    rval = []
    rval.extend( ax.plot(Z0.real, Z0.imag, "gx") )
    rval.extend( ax.plot(ZZ.real.transpose(), ZZ.imag.transpose(), "r.") )
    return rval

def lsf2poly(lsf):
    """Convert line spectral frequencies to prediction filter coefficients

    returns a vector a containing the prediction filter coefficients from a vector lsf of line spectral frequencies.

    .. doctest::

        >>> from spectrum import lsf2poly
        >>> lsf = [0.7842 ,   1.5605  ,  1.8776 ,   1.8984,    2.3593]
        >>> a = lsf2poly(lsf)

    # array([  1.00000000e+00,   6.14837835e-01,   9.89884967e-01,
    # 9.31594056e-05,   3.13713832e-03,  -8.12002261e-03 ])

    .. seealso:: poly2lsf, rc2poly, ac2poly, rc2is
    """
    #   Reference: A.M. Kondoz, "Digital Speech: Coding for Low Bit Rate Communications
    #   Systems" John Wiley & Sons 1994 ,Chapter 4

    # Line spectral frequencies must be real.

    lsf = numpy.array(lsf)

    if max(lsf) > numpy.pi or min(lsf) < 0:
        raise ValueError('Line spectral frequencies must be between 0 and pi.')

    p = len(lsf) # model order

    # Form zeros using the LSFs and unit amplitudes
    z  = numpy.exp(1.j * lsf)

    # Separate the zeros to those belonging to P and Q
    rQ = z[0::2]
    rP = z[1::2]

    # Include the conjugates as well
    rQ = numpy.concatenate((rQ, rQ.conjugate()))
    rP = numpy.concatenate((rP, rP.conjugate()))

    # Form the polynomials P and Q, note that these should be real
    Q  = numpy.poly(rQ);
    P  = numpy.poly(rP);

    # Form the sum and difference filters by including known roots at z = 1 and
    # z = -1

    if p%2:
        # Odd order: z = +1 and z = -1 are roots of the difference filter, P1(z)
        P1 = numpy.convolve(P, [1, 0, -1])
        Q1 = Q
    else:
        # Even order: z = -1 is a root of the sum filter, Q1(z) and z = 1 is a
        # root of the difference filter, P1(z)
        P1 = numpy.convolve(P, [1, -1])
        Q1 = numpy.convolve(Q, [1,  1])

    # Prediction polynomial is formed by averaging P1 and Q1

    a = .5 * (P1+Q1)
    return a[0:-1:1] # do not return last element

def zpk2tf(z, p, k):
    """Return polynomial transfer function representation from zeros
    and poles

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

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

    """
    z = atleast_1d(z)
    k = atleast_1d(k)
    if len(z.shape) > 1:
        temp = poly(z[0])
        b = zeros((z.shape[0], z.shape[1] + 1), temp.dtype.char)
        if len(k) == 1:
            k = [k[0]] * z.shape[0]
        for i in range(z.shape[0]):
            b[i] = k[i] * poly(z[i])
    else:
        b = k * poly(z)
    a = atleast_1d(poly(p))
    return b, a

def lsf_to_lpc(all_lsf):
    if len(all_lsf.shape) < 2:
        all_lsf = all_lsf[None]
    order = all_lsf.shape[1]
    all_lpc = np.zeros((len(all_lsf), order + 1))
    for i in range(len(all_lsf)):
        lsf = all_lsf[i]
        zeros = np.exp(1j * lsf)
        sum_zeros = zeros[::2]
        diff_zeros = zeros[1::2]
        sum_zeros = np.hstack((sum_zeros, np.conj(sum_zeros)))
        diff_zeros = np.hstack((diff_zeros, np.conj(diff_zeros)))
        sum_filt = np.poly(sum_zeros)
        diff_filt = np.poly(diff_zeros)

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

        lpc = .5 * (deconv_sum + deconv_diff)
        # Last coefficient is 0 and not returned
        all_lpc[i] = lpc[:-1]
    return np.squeeze(all_lpc)

def test_infnorm():
	"""Test function for infnorm()"""
	# FIXME M/P test needed
	num, den = np.poly([3, 0.3, 1]), np.poly([2, 0.5, .25])
	H = (num, den)
	Hinf, fmax = infnorm(H)
	assert np.allclose((Hinf, fmax), (np.array([ 1.84888889]), np.array([ 0.50000001])), 
	                   atol=1e-8, rtol=1e-5)

    def impz(self, zeros, poles, L):
        from scipy.signal import lfilter

        a = np.poly(poles)
        b = np.poly(zeros)
        d = np.zeros(L)
        d[0] = 1
        h = lfilter(b, a, d)
        return h

 def test_impL1_num_den(self):
     """Test function for impL1() 2/4"""
     sys1 = (np.array([-.4]), np.array([0, 0]), 1) 
     num = np.poly(sys1[0])
     den = np.poly(sys1[1])
     num[0] *= sys1[2]
     tf = (num, den)
     r3 = ds.impL1(tf, n=10)
     self.assertTrue(np.allclose(self.r2, r3, atol=1e-8, rtol=1e-4))

def getTheta(ar_roots, ma_roots, sigma):
	'''Return theta for Libcarma'''
	ma_poly = sigma*np.poly(ma_roots) #TODO Change * to /
	ar_poly = np.poly(ar_roots)
	theta = np.concatenate((ar_poly[1:], ma_poly))
	p = len(ar_roots)
	q = len(ma_roots)
		
	return theta, p, q

def dual_band(F, P, E, eps, eps_R=1, x1=0.5, map_type=1):
    r"""
    Function to give modified F, P, and E polynomials after lowpass 
    prototype dual band transformation.
    
    :param F:       Polynomial F, i.e., numerator of S11 (in s-domain)
    :param P:       Polynomial P, i.e., numerator of S21 (in s-domain)
    :param E:       polynomial E is the denominator of S11 and S21 (in s-domain)
    :param eps:     Constant term associated with S21
    :param eps_R:   Constant term associated with S11
    :param x1:
    :param map_type:
                  
    :rtype:    
    """

    if (map_type == 1) or (map_type == 2):

        s = sp.Symbol("s")
        if map_type == 1:
            a = -2j / (1 - x1 ** 2)
            b = -1j * (1 + x1 ** 2) / (1 - x1 ** 2)
        elif map_type == 2:
            a = 2j / (1 - x1 ** 2)
            b = 1j * (1 + x1 ** 2) / (1 - x1 ** 2)
        s1 = a * s ** 2 + b

        F = sp.Poly(F.ravel().tolist(), s)
        F1 = sp.simplify(F.subs(s, s1))
        F1 = sp.Poly(F1, s).all_coeffs()
        F1 = I_to_i(F1)

        E = sp.Poly(E.ravel().tolist(), s)
        E1 = sp.simplify(E.subs(s, s1))
        E1 = sp.Poly(E1, s).all_coeffs()
        E1 = I_to_i(E1)

        P = sp.Poly(P.ravel().tolist(), s)
        P1 = sp.simplify(P.subs(s, s1))
        P1 = sp.Poly(P1, s).all_coeffs()
        P1 = I_to_i(P1)

    elif map_type == 3:

        F_roots = np.roots(F.ravel().tolist())
        P_roots = np.roots(P.ravel().tolist())
        F1_roots = Lee_roots_map(F_roots, x1)
        P1_roots = Lee_roots_map(P_roots, x1)
        P1_roots = np.concatenate((P1_roots, np.array([0, 0])))
        F1 = np.poly(F1_roots)
        P1 = np.poly(P1_roots)
        F1 = np.reshape(F1, (len(F1), -1))
        P1 = np.reshape(P1, (len(P1), -1))
        E1 = poly_E(eps, eps_R, F1, P1)[0]

    return F1, P1, E1

def ss2tf(A, B, C, D, input=0):
    """State-space to transfer function.

    Parameters
    ----------
    A, B, C, D : ndarray
        State-space representation of linear system.
    input : int, optional
        For multiple-input systems, the input to use.

    Returns
    -------
    num : 2-D ndarray
        Numerator(s) of the resulting transfer function(s).  `num` has one row
        for each of the system's outputs. Each row is a sequence representation
        of the numerator polynomial.
    den : 1-D ndarray
        Denominator of the resulting transfer function(s).  `den` is a sequence
        representation of the denominator polynomial.

    """
    # transfer function is C (sI - A)**(-1) B + D
    A, B, C, D = map(asarray, (A, B, C, D))
    # Check consistency and make them all rank-2 arrays
    A, B, C, D = abcd_normalize(A, B, C, D)

    nout, nin = D.shape
    if input >= nin:
        raise ValueError("System does not have the input specified.")

    # make MOSI from possibly MOMI system.
    if B.shape[-1] != 0:
        B = B[:, input]
    B.shape = (B.shape[0], 1)
    if D.shape[-1] != 0:
        D = D[:, input]

    try:
        den = poly(A)
    except ValueError:
        den = 1

    if (product(B.shape, axis=0) == 0) and (product(C.shape, axis=0) == 0):
        num = numpy.ravel(D)
        if (product(D.shape, axis=0) == 0) and (product(A.shape, axis=0) == 0):
            den = []
        return num, den

    num_states = A.shape[0]
    type_test = A[:, 0] + B[:, 0] + C[0, :] + D
    num = numpy.zeros((nout, num_states + 1), type_test.dtype)
    for k in range(nout):
        Ck = atleast_2d(C[k, :])
        num[k] = poly(A - dot(B, Ck)) + (D[k] - 1) * den

    return num, den

 def solve(self):
     if len(self.pole) == 0:
         self.a = np.ones((1))
     else:
         self.a = np.poly(self.pole)
     if len(self.zero) == 0:
         self.b = np.ones((1))
     else:
         self.b = np.poly(self.zero)
     self.solved = True

    def init(self, freq=50.0, mod=0.9):
        theta = 2 * np.pi * 50 / self.input.sample_rate
        zero = np.exp(np.array([1j, -1j]) * theta)
        pole = mod * zero

        a, b = np.poly(pole), np.poly(zero)
        #notch_ab = numpy.poly(zero), numpy.poly(pole)
        #notch_ab = scipy.signal.iirfilter(32, [30.0 / 125], btype='low')

        self.gr_block = filter.iir_filter_ffd(b, a, oldstyle=False)

def z_coeff(Poles,Zeros,fs,g,fg,fo = 'none'):
    if fg == np.inf:
        fg = fs/2
    if fo == 'none':
        beta = 1.0
    else:
        beta = f_warp(fo,fs)/fo
    a = np.poly(z_from_f(beta*np.array(Poles),fs))
    b = np.poly(z_from_f(beta*np.array(Zeros),fs))
    gain = 10.**(g/20.)/abs(Fz_at_f(beta*np.array(Poles),beta*np.array(Zeros),fg,fs))
    
    return (a,b*gain)

def invresz(r, p, k, tol=1e-3, rtype='avg'):
    """Compute b(z) and a(z) from partial fraction expansion: r,p,k

    If M = len(b) and N = len(a)

            b(z)     b[0] + b[1] z**(-1) + ... + b[M-1] z**(-M+1)
    H(z) = ------ = ----------------------------------------------
            a(z)     a[0] + a[1] z**(-1) + ... + a[N-1] z**(-N+1)

                 r[0]                   r[-1]
         = --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ...
           (1-p[0]z**(-1))         (1-p[-1]z**(-1))

    If there are any repeated roots (closer than tol), then the partial
    fraction expansion has terms like

               r[i]              r[i+1]                    r[i+n-1]
          -------------- + ------------------ + ... + ------------------
          (1-p[i]z**(-1))  (1-p[i]z**(-1))**2         (1-p[i]z**(-1))**n

    See also
    --------
    residuez, poly, polyval, unique_roots

    """
    extra = asarray(k)
    p, indx = cmplx_sort(p)
    r = take(r, indx, 0)
    pout, mult = unique_roots(p, tol=tol, rtype=rtype)
    p = []
    for k in range(len(pout)):
        p.extend([pout[k]] * mult[k])
    a = atleast_1d(poly(p))
    if len(extra) > 0:
        b = polymul(extra, a)
    else:
        b = [0]
    indx = 0
    brev = asarray(b)[::-1]
    for k in range(len(pout)):
        temp = []
        # Construct polynomial which does not include any of this root
        for l in range(len(pout)):
            if l != k:
                temp.extend([pout[l]] * mult[l])
        for m in range(mult[k]):
            t2 = temp[:]
            t2.extend([pout[k]] * (mult[k] - m - 1))
            brev = polyadd(brev, (r[indx] * poly(t2))[::-1])
            indx += 1
    b = real_if_close(brev[::-1])
    return b, a

def invres(r,p,k,tol=1e-3,rtype='avg'):
    """Compute b(s) and a(s) from partial fraction expansion: r,p,k

    If M = len(b) and N = len(a)

            b(s)     b[0] x**(M-1) + b[1] x**(M-2) + ... + b[M-1]
    H(s) = ------ = ----------------------------------------------
            a(s)     a[0] x**(N-1) + a[1] x**(N-2) + ... + a[N-1]

             r[0]       r[1]             r[-1]
         = -------- + -------- + ... + --------- + k(s)
           (s-p[0])   (s-p[1])         (s-p[-1])

    If there are any repeated roots (closer than tol), then the partial
    fraction expansion has terms like

            r[i]      r[i+1]              r[i+n-1]
          -------- + ----------- + ... + -----------
          (s-p[i])  (s-p[i])**2          (s-p[i])**n

    See Also
    --------
    residue, poly, polyval, unique_roots

    """
    extra = k
    p, indx = cmplx_sort(p)
    r = take(r,indx,0)
    pout, mult = unique_roots(p,tol=tol,rtype=rtype)
    p = []
    for k in range(len(pout)):
        p.extend([pout[k]]*mult[k])
    a = atleast_1d(poly(p))
    if len(extra) > 0:
        b = polymul(extra,a)
    else:
        b = [0]
    indx = 0
    for k in range(len(pout)):
        temp = []
        for l in range(len(pout)):
            if l != k:
                temp.extend([pout[l]]*mult[l])
        for m in range(mult[k]):
            t2 = temp[:]
            t2.extend([pout[k]]*(mult[k]-m-1))
            b = polyadd(b,r[indx]*poly(t2))
            indx += 1
    b = real_if_close(b)
    while allclose(b[0], 0, rtol=1e-14) and (b.shape[-1] > 1):
        b = b[1:]
    return b, a

def zp_expr(s,Z,P):
    pZ = np.poly(Z)
    pP = np.poly(P)

    a = np.polyval(pZ, s)
    b = np.polyval(pP, s)

    print "zp_expr"

    print "  a",a
    print "  b",b

    return a / b