This version of ntf_fir_from_q0 passes the unit tests.

Fixed ntf_fir_from_q0

def ntf_fir_from_q0(q0, H_inf=1.5, normalize="auto", **options):
    """
    Synthesize FIR NTF from quadratic form expressing noise weighting.

    Parameters
    ----------
    q0 : ndarray
        first row of the Toeplitz symmetric matrix defining the quadratic form
    H_inf : real, optional
        Max peak NTF gain, defaults to 1.5, used to enforce the Lee criterion
    normalize : string or real, optional
        Normalization to apply to the quadratic form used in the NTF
        selection. Defaults to 'auto' which means setting the top left entry
        in the matrix Q defining the quadratic form to 1.

    Returns
    -------
    ntf : ndarray
        FIR NTF in zpk form

    Other parameters
    ----------------
    show_progress : bool, optional
        provide extended output, default is True
    fix_pos : bool, optional
        fix quadratic form for positive definiteness. Numerical noise
        may make it not positive definite leading to errors. Default is True
    cvxopt_xxx : various type, optional
        Parameters prefixed by ``cvxopt_`` are passed to the ``cvxopt``
        optimizer. Allowed options are:

        ``cvxopt_maxiters``
            Maximum number of iterations (defaults to 100)
        ``cvxopt_abstol``
            Absolute accuracy (defaults to 1e-7)
        ``cvxopt_reltol``
            Relative accuracy (defaults to 1e-6)
        ``cvxopt_feastol``
            Tolerance for feasibility conditions (dtimeefaults to 1e-6)

        Do not use other options since they could break ``cvxpy`` in
        unexpected ways. Defaults can be set by changing the function
        ``default_options`` attribute.

    See Also
    --------

    Check the documentation of ``cvxopt`` for further information.
    """
    # Manage optional parameters
    opts = ntf_fir_from_q0.default_options.copy()
    opts.update(options)
    o = split_options(opts, ['cvxopt_'], ['show_progress', 'fix_pos'])
    # Do the computation
    order = q0.shape[0]-1
    if normalize == 'auto':
        q0 = q0/q0[0]
    elif normalize is not None:
        q0 = q0*normalize
    Q = la.toeplitz(q0)
    d, v = np.linalg.eigh(Q)
    if o.get('fix_pos', True):
        d = d/np.max(d)
        d[d < 0] = 0.
    qs = np.matrix(mdot(v, np.diag(np.sqrt(d)), np.linalg.inv(v)))
    br = cvx.Variable(order, 1, name='br')
    b = cvx.vstack(1, br)
    target = cvx.Minimize(cvx.norm2(qs*b))
    X = cvx.Semidef(order, name='X')
    A = np.matrix(np.eye(order, order, 1))
    B = np.vstack((np.zeros((order-1, 1)), 1.))
    C = br[::-1].T
    D = np.matrix(1.)
    M1 = A.T*X
    M2 = M1*B
    dim = X.size[1] + M2.size[1] + C.T.size[1]
    M = cvx.vstack(
        cvx.hstack(M1*A-X, M2, C.T),
        cvx.hstack(M2.T, B.T*X*B-H_inf**2, D),
        cvx.hstack(C, D, np.matrix(-1.))
        )
    Mcopy = -cvx.Semidef(order+2)
    constraints = [
        cvx.diag(Mcopy) == cvx.diag(M),
        cvx.upper_tri(Mcopy) == cvx.upper_tri(M),
    ]
    p = cvx.Problem(target, constraints)
    p.solve(solver=cvx.CVXOPT, verbose=o.get('show_progress', True),
            kktsolver="chol",
            **strip_options(o, 'cvxopt_'))
    ntf_ir = np.hstack((1, np.asarray(br.value.T)[0]))
    return (np.roots(ntf_ir), np.zeros(order), 1.)

ntf_fir_from_q0.default_options = {'cvxopt_maxiters': 100,
                                   'cvxopt_abstol': 1e-7,
                                   'cvxopt_reltol': 1e-6,
                                   'cvxopt_feastol': 1e-6,
                                   'show_progress': True,
                                   'fix_pos': True}