def arcovar_marple(x, order):
    r"""Estimate AR model parameters using covariance method

    This implementation is based on [Marple]_. This code is far more
    complicated and slower than :func:`arcovar` function, which is now the official version.
    See :func:`arcovar` for a detailed description of Covariance method.

    This function should be used in place of arcovar only if order<=4, for
    which :func:`arcovar` does not work.

    Fast algorithm for the solution of the covariance least squares normal
    equations from Marple.

    :param array X:  Array of complex data samples
    :param int oder: Order of linear prediction model

    :return:
        * AF   - Array of complex forward linear prediction coefficients
        * PF   - Real forward linear prediction variance at order IP
        * AB   - Array of complex backward linear prediction coefficients
        * PB   - Real backward linear prediction variance at order IP
        * PV   - store linear prediction coefficients

    .. note:: this code and the original code in Marple diverge for ip>10.
        it seems that this is related to single precision used with
        complex type in fortran whereas numpy uses double precision for
        complex type.

    :validation: the AR parameters are the same as those returned by
        a completely different function :func:`arcovar`.

    :References: [Marple]_
    """
    assert len(x) >= order, "X must be dimensioned >=N"

    #   ----------------------------------------------------- Initialization
    x = np.array(x)
    N = len(x)


    # Equations 8.C.42
    r0 = sum(abs(x)**2.)
    r1 = abs(x[0])**2
    rN = abs(x[N-1])**2

    pf = r0 - r1
    pb = r0 - rN

    delta = 1. - r1 / r0
    gamma = 1. - rN / r0

    c = np.zeros(N, dtype=complex)
    d = np.zeros(N, dtype=complex)
    r = np.zeros(N, dtype=complex)
    af = np.zeros(N, dtype=complex)
    ab = np.zeros(N, dtype=complex)

    c[0] = x[N-1].conjugate() / r0
    d[0] = x[0].conjugate() / r0

    # special case
    if order == 0:
        pf = r0 / float(N)
        pb = pf
        return af, pf, ab, pb, 0

    # ----------------------------------------------------------  MAIN LOOP

    #ip +1 because we want to enter in the loop to run the first part of the code.
    pbv = []
    for m in range(0, order+1):
        logging.debug('----------------------------m=', m)
        logging.debug(c[0:2])
        logging.debug(d[0:2])
        r1 = 1./pf
        r2 = 1./pb
        r3 = 1./delta
        r4 = 1./gamma

        #logging.debug('starting r1r2r3r4=', r1, r2, r3, r4, pf, pb, delta, gamma)

        #Order update: AF and AB vectors ; time update: C and D vectors
        temp = 0.+0.j
        for k in range(m+1, N):
            temp = temp + x[k]*x[k-m-1].conjugate()

        r[m] =  temp.conjugate()
        theta = x[0] * c[m]
        #print(('theta', theta))
        #    print(('c=', c[0:2]))
        #    print(('d=', d[0:2]))
        if m == 0:
            pass
        else:
            for k in range(0, m):
                theta = theta + x[m-k] * c[k]                   # Eq. (8.C.39)
                r[k] = r[k] - x[N-m-1] * x[N-m+k].conjugate() # Eq. (8.C.32)
                temp = temp + af[m-k-1] * r[k].conjugate()
                #print 'loop1 k=', k
                #print '         theta=',theta, 'r[k]=',r[k], 'temp=', temp
                #print '         c=',c[k], 'af=',af[m-k-1]
        """if m > 0:
            if debug:
                print((m, N-m))
                print(('Xk=0',x[m-0],x[N-m-1], x[N-m+0]))
            if m > 1:
                if debug:
                    print('Xk=1',x[m-1],x[N-m-1], x[N-m+1])
        """
        c1 = -temp * r2
        c2 = -r1 * temp.conjugate()
        c3 = theta * r3
        c4 = r4 *theta.conjugate()
        #if debug:
        #    print('c1 c2 c3 c4 before af=',c1 ,c2 ,c3 ,c4)
        af[m] = c1                    #             ! Eq. (8.C.19)
        ab[m] = c2                    #             ! Eq. (8.C.22)
        save = c[m]
        c[m] = save + c3*d[m]
        d[m] = d[m] + c4*save
        #if debug:
        #    print('res',m,'af[m]=',af[m], ab[m], save, 'temp=',temp)

        if m == 0:
            pass
        else:
            #if debug:print('af before', af[0:2])
            for k in range(0, m):
                save = af[k]
                af[k] = save + c1 * ab[m-k-1] # Eq. (8.C.18)
                ab[m-k-1] = ab[m-k-1] + c2 * save   # Eq. (8.C.21)

                save = c[k]
                c[k] = save + c3*d[k]       # Eq. (8.C.37)
                d[k] = d[k] + c4*save       # Eq. (8.C.38)
                #if debug:
                #    print('loop2 k=', k)
                #    print('      af[k]=', af[k])
                #    print('      ab[m-k-1]=', ab[m-k-1])
                #    print('      c[k]=', c[k])
                #    print('      d[k]=', d[k])
            #if debug:
            #    print('af after=', af[0:2])
            #    print('ab=', ab[0:2])

        r5 = temp.real**2 + temp.imag**2
        pf = pf - r5*r2         # Eq. (8.C.20)
        pb = pb - r5*r1         # Eq. (8.C.23)
        r5 = theta.real**2 + theta.imag**2
        delta = delta - r5*r4               # Eq. (8.C.39)
        gamma = gamma - r5*r3               # Eq. (8.C.40)
        #if debug:
        #    print('r5r2r1deltagamma', r5, r2, r1 , delta, gamma)
        #    print('pf before norm', pf, pb, N-m-1)
        if m != order-1:
            pass
        else:
            pf = pf / float(N-m-1)
            pb = pb / float(N-m-1)
            #if debug:
            #    print('ENDING', N-m-1)
            break

        #if debug:
        #    print('pf and pb', pf, pb)
        if pf > 0 and pb > 0:
            pass
        else:
            ValueError("Negative PF or PB value")
        if (delta > 0. and delta <=1 and gamma > 0. and gamma <=1):
            pass
        else:
            ValueError("Invalid delta or gamma value")

        #C   Time update:  AF and AB vectors; order update:  C and D vectors

        r1 = 1./pf
        r2 = 1./pb
        r3 = 1./delta
        r4 = 1./gamma
        #if debug:
        #    print('--------time update', r1, r2, r3, r4, m+1, N-m-1, x[m+1], x[N-m-2])
        ef = x[m+1]
        eb = x[(N-1)-m-1]


        for k in range(0,m+1):
            #print 'k=', k, 'ef=', ef, ' eb=',eb,' af=',af[k], ab[k]
            #print x[m-k],x[N-m+k-1]
            ef = ef + af[k] * x[m-k]             # Eq. (8.C.1)
            eb = eb + ab[k] * x[N-m+k-1]                   # Eq. (8.C.2)

        #ef = sum(af)

        #if debug:
        #    print('efweb', ef , eb)
        c1 = ef*r3
        c2 = eb*r4
        c3 = eb.conjugate() * r2
        c4 = ef.conjugate() * r1
        #if debug:
        #    print('c1c2c3c4', c1, c2, c3, c4)
        #    print('af before', af[0:2])
        for k in range(m, -1, -1):
            save = af[k]
            af[k] = save + c1 * d[k]                    #  Eq. (8.C.33)
            d[k+1] = d[k] + c4 * save                    # Eq. (8.C.25)
            save = ab[k]
            ab[k] = save + c2 * c[m-k]                 # Eq. (8.C.35)
            c[m-k] = c[m-k] + c3 * save              # Eq. (8.C.24)
        #if debug:
        #    print('af after', af[0:2])
        #    print('d', d[0:2])
        #    print('ab', ab[0:2])
        #    print('c', c[0:2])
        #if debug:print('Pb before', pf, pb)
        c[m+1] = c3
        d[0] = c4
        #r5 = abs(ef)**2
        r5 = ef.real**2 + ef.imag**2
        pf = pf - r5 * r3                              # Eq. (8.C.34)

        delta = delta-r5 * r1                        # Eq. (8.C.30)
        #r5 = abs(eb)**2
        r5 = eb.real**2 + eb.imag**2
        pb = pb - r5 * r4                              # Eq. (8.C.36)
        #if debug:
        #    print('Pb---------------------', m, pb, r5, r4)
        gamma = gamma-r5*r2                        # Eq. (8.C.31)
        pbv.append(pb)

        if (pf > 0. and pb > 0.):
            pass
        else:
            ValueError("Negative PF or PB value")
        #if debug:
        #    print(delta, gamma)
        if (delta > 0. and delta <= 1.) and (gamma > 0. and gamma <= 1.):
            pass
        else:
            ValueError("Invalid delta or gamma value")


    #af=array of forward coeff
    #ab=array of barward coeff
    #pb=backward variance
    #pf=forward variance
    return af, pf, ab, pb, pbv
def arcovar_marple(x, order):
    r"""Estimate AR model parameters using covariance method

    This implementation is based on [Marple]_. This code is far more
    complicated and slower than :func:`arcovar` function, which is now the official version.
    See :func:`arcovar` for a detailed description of Covariance method.

    This function should be used in place of arcovar only if order<=4, for
    which :func:`arcovar` does not work.

    Fast algorithm for the solution of the covariance least squares normal
    equations from Marple.

    :param array X:  Array of complex data samples
    :param int oder: Order of linear prediction model

    :return:
        * AF   - Array of complex forward linear prediction coefficients
        * PF   - Real forward linear prediction variance at order IP
        * AB   - Array of complex backward linear prediction coefficients
        * PB   - Real backward linear prediction variance at order IP
        * PV   - store linear prediction coefficients

    .. note:: this code and the original code in Marple diverge for ip>10.
        it seems that this is related to single precision used with
        complex type in fortran whereas numpy uses double precision for
        complex type.

    :validation: the AR parameters are the same as those returned by
        a completely different function :func:`arcovar`.

    :References: [Marple]_
    """
    assert len(x) >= order, "X must be dimensioned >=N"

    #   ----------------------------------------------------- Initialization
    x = np.array(x)
    N = len(x)


    # Equations 8.C.42
    r0 = sum(abs(x)**2.)
    r1 = abs(x[0])**2
    rN = abs(x[N-1])**2

    pf = r0 - r1
    pb = r0 - rN

    delta = 1. - r1 / r0
    gamma = 1. - rN / r0

    c = np.zeros(N, dtype=complex)
    d = np.zeros(N, dtype=complex)
    r = np.zeros(N, dtype=complex)
    af = np.zeros(N, dtype=complex)
    ab = np.zeros(N, dtype=complex)

    c[0] = x[N-1].conjugate() / r0
    d[0] = x[0].conjugate() / r0

    # special case
    if order == 0:
        pf = r0 / float(N)
        pb = pf
        return af, pf, ab, pb, 0

    # ----------------------------------------------------------  MAIN LOOP

    #ip +1 because we want to enter in the loop to run the first part of the code.
    pbv = []
    for m in range(0, order+1):
        logging.debug('----------------------------m=', m)
        logging.debug(c[0:2])
        logging.debug(d[0:2])
        r1 = 1./pf
        r2 = 1./pb
        r3 = 1./delta
        r4 = 1./gamma

        #logging.debug('starting r1r2r3r4=', r1, r2, r3, r4, pf, pb, delta, gamma)

        #Order update: AF and AB vectors ; time update: C and D vectors
        temp = 0.+0.j
        for k in range(m+1, N):
            temp = temp + x[k]*x[k-m-1].conjugate()

        r[m] =  temp.conjugate()
        theta = x[0] * c[m]
        #print(('theta', theta))
        #    print(('c=', c[0:2]))
        #    print(('d=', d[0:2]))
        if m == 0:
            pass
        else:
            for k in range(0, m):
                theta = theta + x[m-k] * c[k]                   # Eq. (8.C.39)
                r[k] = r[k] - x[N-m-1] * x[N-m+k].conjugate() # Eq. (8.C.32)
                temp = temp + af[m-k-1] * r[k].conjugate()
                #print 'loop1 k=', k
                #print '         theta=',theta, 'r[k]=',r[k], 'temp=', temp
                #print '         c=',c[k], 'af=',af[m-k-1]
        """if m > 0:
            if debug:
                print((m, N-m))
                print(('Xk=0',x[m-0],x[N-m-1], x[N-m+0]))
            if m > 1:
                if debug:
                    print('Xk=1',x[m-1],x[N-m-1], x[N-m+1])
        """
        c1 = -temp * r2
        c2 = -r1 * temp.conjugate()
        c3 = theta * r3
        c4 = r4 *theta.conjugate()
        #if debug:
        #    print('c1 c2 c3 c4 before af=',c1 ,c2 ,c3 ,c4)
        af[m] = c1                    #             ! Eq. (8.C.19)
        ab[m] = c2                    #             ! Eq. (8.C.22)
        save = c[m]
        c[m] = save + c3*d[m]
        d[m] = d[m] + c4*save
        #if debug:
        #    print('res',m,'af[m]=',af[m], ab[m], save, 'temp=',temp)

        if m == 0:
            pass
        else:
            #if debug:print('af before', af[0:2])
            for k in range(0, m):
                save = af[k]
                af[k] = save + c1 * ab[m-k-1] # Eq. (8.C.18)
                ab[m-k-1] = ab[m-k-1] + c2 * save   # Eq. (8.C.21)

                save = c[k]
                c[k] = save + c3*d[k]       # Eq. (8.C.37)
                d[k] = d[k] + c4*save       # Eq. (8.C.38)
                #if debug:
                #    print('loop2 k=', k)
                #    print('      af[k]=', af[k])
                #    print('      ab[m-k-1]=', ab[m-k-1])
                #    print('      c[k]=', c[k])
                #    print('      d[k]=', d[k])
            #if debug:
            #    print('af after=', af[0:2])
            #    print('ab=', ab[0:2])

        r5 = temp.real**2 + temp.imag**2
        pf = pf - r5*r2         # Eq. (8.C.20)
        pb = pb - r5*r1         # Eq. (8.C.23)
        r5 = theta.real**2 + theta.imag**2
        delta = delta - r5*r4               # Eq. (8.C.39)
        gamma = gamma - r5*r3               # Eq. (8.C.40)
        #if debug:
        #    print('r5r2r1deltagamma', r5, r2, r1 , delta, gamma)
        #    print('pf before norm', pf, pb, N-m-1)
        if m != order-1:
            pass
        else:
            pf = pf / float(N-m-1)
            pb = pb / float(N-m-1)
            #if debug:
            #    print('ENDING', N-m-1)
            break

        #if debug:
        #    print('pf and pb', pf, pb)
        if pf > 0 and pb > 0:
            pass
        else:
            ValueError("Negative PF or PB value")
        if (delta > 0. and delta <=1 and gamma > 0. and gamma <=1):
            pass
        else:
            ValueError("Invalid delta or gamma value")

        #C   Time update:  AF and AB vectors; order update:  C and D vectors

        r1 = 1./pf
        r2 = 1./pb
        r3 = 1./delta
        r4 = 1./gamma
        #if debug:
        #    print('--------time update', r1, r2, r3, r4, m+1, N-m-1, x[m+1], x[N-m-2])
        ef = x[m+1]
        eb = x[(N-1)-m-1]


        for k in range(0,m+1):
            #print 'k=', k, 'ef=', ef, ' eb=',eb,' af=',af[k], ab[k]
            #print x[m-k],x[N-m+k-1]
            ef = ef + af[k] * x[m-k]             # Eq. (8.C.1)
            eb = eb + ab[k] * x[N-m+k-1]                   # Eq. (8.C.2)

        #ef = sum(af)

        #if debug:
        #    print('efweb', ef , eb)
        c1 = ef*r3
        c2 = eb*r4
        c3 = eb.conjugate() * r2
        c4 = ef.conjugate() * r1
        #if debug:
        #    print('c1c2c3c4', c1, c2, c3, c4)
        #    print('af before', af[0:2])
        for k in range(m, -1, -1):
            save = af[k]
            af[k] = save + c1 * d[k]                    #  Eq. (8.C.33)
            d[k+1] = d[k] + c4 * save                    # Eq. (8.C.25)
            save = ab[k]
            ab[k] = save + c2 * c[m-k]                 # Eq. (8.C.35)
            c[m-k] = c[m-k] + c3 * save              # Eq. (8.C.24)
        #if debug:
        #    print('af after', af[0:2])
        #    print('d', d[0:2])
        #    print('ab', ab[0:2])
        #    print('c', c[0:2])
        #if debug:print('Pb before', pf, pb)
        c[m+1] = c3
        d[0] = c4
        #r5 = abs(ef)**2
        r5 = ef.real**2 + ef.imag**2
        pf = pf - r5 * r3                              # Eq. (8.C.34)

        delta = delta-r5 * r1                        # Eq. (8.C.30)
        #r5 = abs(eb)**2
        r5 = eb.real**2 + eb.imag**2
        pb = pb - r5 * r4                              # Eq. (8.C.36)
        #if debug:
        #    print('Pb---------------------', m, pb, r5, r4)
        gamma = gamma-r5*r2                        # Eq. (8.C.31)
        pbv.append(pb)

        if (pf > 0. and pb > 0.):
            pass
        else:
            ValueError("Negative PF or PB value")
        #if debug:
        #    print(delta, gamma)
        if (delta > 0. and delta <= 1.) and (gamma > 0. and gamma <= 1.):
            pass
        else:
            ValueError("Invalid delta or gamma value")


    #af=array of forward coeff
    #ab=array of barward coeff
    #pb=backward variance
    #pf=forward variance
    return af, pf, ab, pb, pbv
def arcovar_marple ( x , order ) : assert len ( x ) >= order , "X must be dimensioned >=N" #   ----------------------------------------------------- Initialization x = np . array ( x ) N = len ( x ) # Equations 8.C.42 r0 = sum ( abs ( x ) ** 2. ) r1 = abs ( x [ 0 ] ) ** 2 rN = abs ( x [ N - 1 ] ) ** 2 pf = r0 - r1 pb = r0 - rN delta = 1. - r1 / r0 gamma = 1. - rN / r0 c = np . zeros ( N , dtype = complex ) d = np . zeros ( N , dtype = complex ) r = np . zeros ( N , dtype = complex ) af = np . zeros ( N , dtype = complex ) ab = np . zeros ( N , dtype = complex ) c [ 0 ] = x [ N - 1 ] . conjugate ( ) / r0 d [ 0 ] = x [ 0 ] . conjugate ( ) / r0 # special case if order == 0 : pf = r0 / float ( N ) pb = pf return af , pf , ab , pb , 0 # ----------------------------------------------------------  MAIN LOOP #ip +1 because we want to enter in the loop to run the first part of the code. pbv = [ ] for m in range ( 0 , order + 1 ) : logging . debug ( '----------------------------m=' , m ) logging . debug ( c [ 0 : 2 ] ) logging . debug ( d [ 0 : 2 ] ) r1 = 1. / pf r2 = 1. / pb r3 = 1. / delta r4 = 1. / gamma #logging.debug('starting r1r2r3r4=', r1, r2, r3, r4, pf, pb, delta, gamma) #Order update: AF and AB vectors ; time update: C and D vectors temp = 0. + 0.j for k in range ( m + 1 , N ) : temp = temp + x [ k ] * x [ k - m - 1 ] . conjugate ( ) r [ m ] = temp . conjugate ( ) theta = x [ 0 ] * c [ m ] #print(('theta', theta)) #    print(('c=', c[0:2])) #    print(('d=', d[0:2])) if m == 0 : pass else : for k in range ( 0 , m ) : theta = theta + x [ m - k ] * c [ k ] # Eq. (8.C.39) r [ k ] = r [ k ] - x [ N - m - 1 ] * x [ N - m + k ] . conjugate ( ) # Eq. (8.C.32) temp = temp + af [ m - k - 1 ] * r [ k ] . conjugate ( ) #print 'loop1 k=', k #print '         theta=',theta, 'r[k]=',r[k], 'temp=', temp #print '         c=',c[k], 'af=',af[m-k-1] """if m > 0:
            if debug:
                print((m, N-m))
                print(('Xk=0',x[m-0],x[N-m-1], x[N-m+0]))
            if m > 1:
                if debug:
                    print('Xk=1',x[m-1],x[N-m-1], x[N-m+1])
        """ c1 = - temp * r2 c2 = - r1 * temp . conjugate ( ) c3 = theta * r3 c4 = r4 * theta . conjugate ( ) #if debug: #    print('c1 c2 c3 c4 before af=',c1 ,c2 ,c3 ,c4) af [ m ] = c1 #             ! Eq. (8.C.19) ab [ m ] = c2 #             ! Eq. (8.C.22) save = c [ m ] c [ m ] = save + c3 * d [ m ] d [ m ] = d [ m ] + c4 * save #if debug: #    print('res',m,'af[m]=',af[m], ab[m], save, 'temp=',temp) if m == 0 : pass else : #if debug:print('af before', af[0:2]) for k in range ( 0 , m ) : save = af [ k ] af [ k ] = save + c1 * ab [ m - k - 1 ] # Eq. (8.C.18) ab [ m - k - 1 ] = ab [ m - k - 1 ] + c2 * save # Eq. (8.C.21) save = c [ k ] c [ k ] = save + c3 * d [ k ] # Eq. (8.C.37) d [ k ] = d [ k ] + c4 * save # Eq. (8.C.38) #if debug: #    print('loop2 k=', k) #    print('      af[k]=', af[k]) #    print('      ab[m-k-1]=', ab[m-k-1]) #    print('      c[k]=', c[k]) #    print('      d[k]=', d[k]) #if debug: #    print('af after=', af[0:2]) #    print('ab=', ab[0:2]) r5 = temp . real ** 2 + temp . imag ** 2 pf = pf - r5 * r2 # Eq. (8.C.20) pb = pb - r5 * r1 # Eq. (8.C.23) r5 = theta . real ** 2 + theta . imag ** 2 delta = delta - r5 * r4 # Eq. (8.C.39) gamma = gamma - r5 * r3 # Eq. (8.C.40) #if debug: #    print('r5r2r1deltagamma', r5, r2, r1 , delta, gamma) #    print('pf before norm', pf, pb, N-m-1) if m != order - 1 : pass else : pf = pf / float ( N - m - 1 ) pb = pb / float ( N - m - 1 ) #if debug: #    print('ENDING', N-m-1) break #if debug: #    print('pf and pb', pf, pb) if pf > 0 and pb > 0 : pass else : ValueError ( "Negative PF or PB value" ) if ( delta > 0. and delta <= 1 and gamma > 0. and gamma <= 1 ) : pass else : ValueError ( "Invalid delta or gamma value" ) #C   Time update:  AF and AB vectors; order update:  C and D vectors r1 = 1. / pf r2 = 1. / pb r3 = 1. / delta r4 = 1. / gamma #if debug: #    print('--------time update', r1, r2, r3, r4, m+1, N-m-1, x[m+1], x[N-m-2]) ef = x [ m + 1 ] eb = x [ ( N - 1 ) - m - 1 ] for k in range ( 0 , m + 1 ) : #print 'k=', k, 'ef=', ef, ' eb=',eb,' af=',af[k], ab[k] #print x[m-k],x[N-m+k-1] ef = ef + af [ k ] * x [ m - k ] # Eq. (8.C.1) eb = eb + ab [ k ] * x [ N - m + k - 1 ] # Eq. (8.C.2) #ef = sum(af) #if debug: #    print('efweb', ef , eb) c1 = ef * r3 c2 = eb * r4 c3 = eb . conjugate ( ) * r2 c4 = ef . conjugate ( ) * r1 #if debug: #    print('c1c2c3c4', c1, c2, c3, c4) #    print('af before', af[0:2]) for k in range ( m , - 1 , - 1 ) : save = af [ k ] af [ k ] = save + c1 * d [ k ] #  Eq. (8.C.33) d [ k + 1 ] = d [ k ] + c4 * save # Eq. (8.C.25) save = ab [ k ] ab [ k ] = save + c2 * c [ m - k ] # Eq. (8.C.35) c [ m - k ] = c [ m - k ] + c3 * save # Eq. (8.C.24) #if debug: #    print('af after', af[0:2]) #    print('d', d[0:2]) #    print('ab', ab[0:2]) #    print('c', c[0:2]) #if debug:print('Pb before', pf, pb) c [ m + 1 ] = c3 d [ 0 ] = c4 #r5 = abs(ef)**2 r5 = ef . real ** 2 + ef . imag ** 2 pf = pf - r5 * r3 # Eq. (8.C.34) delta = delta - r5 * r1 # Eq. (8.C.30) #r5 = abs(eb)**2 r5 = eb . real ** 2 + eb . imag ** 2 pb = pb - r5 * r4 # Eq. (8.C.36) #if debug: #    print('Pb---------------------', m, pb, r5, r4) gamma = gamma - r5 * r2 # Eq. (8.C.31) pbv . append ( pb ) if ( pf > 0. and pb > 0. ) : pass else : ValueError ( "Negative PF or PB value" ) #if debug: #    print(delta, gamma) if ( delta > 0. and delta <= 1. ) and ( gamma > 0. and gamma <= 1. ) : pass else : ValueError ( "Invalid delta or gamma value" ) #af=array of forward coeff #ab=array of barward coeff #pb=backward variance #pf=forward variance return af , pf , ab , pb , pbv
r"""Estimate AR model parameters using covariance method

    This implementation is based on [Marple]_. This code is far more
    complicated and slower than :func:`arcovar` function, which is now the official version.
    See :func:`arcovar` for a detailed description of Covariance method.

    This function should be used in place of arcovar only if order<=4, for
    which :func:`arcovar` does not work.

    Fast algorithm for the solution of the covariance least squares normal
    equations from Marple.

    :param array X:  Array of complex data samples
    :param int oder: Order of linear prediction model

    :return:
        * AF   - Array of complex forward linear prediction coefficients
        * PF   - Real forward linear prediction variance at order IP
        * AB   - Array of complex backward linear prediction coefficients
        * PB   - Real backward linear prediction variance at order IP
        * PV   - store linear prediction coefficients

    .. note:: this code and the original code in Marple diverge for ip>10.
        it seems that this is related to single precision used with
        complex type in fortran whereas numpy uses double precision for
        complex type.

    :validation: the AR parameters are the same as those returned by
        a completely different function :func:`arcovar`.

    :References: [Marple]_
r Estimate AR model parameters using covariance method