src/spectrum/modcovar.py
def modcovar_marple (X,IP):
    """Fast algorithm for the solution of the modified covariance least squares normal equations.

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

    :param X:    - Array of complex data samples X(1) through X(N)
    :param int IP:   - Order of linear prediction model (integer)

    :return:
        * P    - Real linear prediction variance at order IP
        * A    - Array of complex linear prediction coefficients
        * ISTAT - Integer status indicator at time of exit
              0. for normal exit (no numerical ill-conditioning)
              1. if P is not a positive value
              2. if DELTA' and GAMMA' do not lie in the range 0 to 1
              3. if P' is not a positive value
              4. if DELTA and GAMMA do not lie in the range 0 to 1


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

    .. note:: validation. results similar to test example in Marple but
        starts to differ for    ip~8. with ratio of 0.975 for ip=15 probably
        due to precision.


    :References: [Marple]_
    """
    Pv = []
    N = len(X)
    A = np.zeros(N, dtype=complex)
    D = np.zeros(N, dtype=complex)
    C = np.zeros(N, dtype=complex)
    R = np.zeros(N, dtype=complex)
    #   Initialization
    R1=0.
    for K in range(1, N-1):
        R1=R1 + 2.*(X[K].real**2 + X[K].imag**2)
    R2 = X[0].real**2 + X[0].imag**2
    R3 = X[N-1].real**2 + X[N-1].imag**2
    R4 = 1. / (R1 + 2. * (R2 + R3))
    P = R1 + R2 + R3
    DELTA = 1. - R2 * R4
    GAMMA = 1. - R3 * R4
    LAMBDA = (X[0] * X[N-1]).conjugate()*R4
    C[0] = X[N-1] * R4
    D[0] = X[0].conjugate() * R4

    M = 0
    if (IP ==0):
        P = (.5*R1+R2+R3)/float(N)
        return [], P, []

    #C   Main loop

    for M in range(0, IP):
        #print '---------------------', M
        SAVE1 = 0+0j
        for K in range(M+1, N):
            SAVE1 = SAVE1 + X[K]*X[K-M-1].conjugate()
        SAVE1 *= 2.
        R[M] = SAVE1.conjugate()
        THETA = X[N-1]*D[0]
        PSI=X[N-1]*C[0]
        XI = X[0].conjugate() * D[0]
        if M==0:
            pass
        else:
            for K in range(0, M):
                THETA=THETA+X[N-K-2]*D[K+1]    # Eq. (8.D.45)
                PSI = PSI + X[N-K-2]*C[K+1]    # Eq. (8.D.45)
                XI = XI + X[K+1].conjugate() * D[K+1] # Eq. (8.D.45)
                R[K] = R[K]-X[N-M-1] * X[N+1-M+K-1].conjugate() - X[M].conjugate() * X[M-K-1] # Eq. (8.D.37)
                SAVE1=SAVE1+R[K].conjugate()*A[M-K-1]             # Eq. (8.D.24)
                #print 'withini loop', K, THETA, PSI, XI, R[K], SAVE1
                #print 'x   -------', C[K], C[K+1]
                #print 'x   -------', D[K], D[K+1]
                #print 'x   -------', N-K-2, K+1,N-M-1,  M, M-K-1
                #print 'x   -------', X[N-K-2], X[K+1],X[N-M-1],  X[M], X[M-K-1]


        #  Order update of A vector

        C1 = -SAVE1/P
        A[M]=C1                                   # Eq. (8.D.23)
        P=P*(1.-C1.real**2-C1.imag**2)            # Eq. (8.D.25)
        #print 'C1=' , C1, 'A[M]', A[M], 'P=',P
        if M==0:
            pass
        else:
            #for k in range(0, (M-1)/2+1):
            #print 'BEFORE', A[0:4]
            for K in range(0, (M+1)//2):
                MK = M-K-1
                SAVE1=A[K]
                A[K]=SAVE1+C1*A[MK].conjugate() # Eq. (8.D.22)
                #print 'A uopdate---->',M, K, MK, A[K], A[MK], SAVE1, C1
                if (K != MK):
                    #print 'ENTERI IF', K, MK, A[MK], C1, SAVE1
                    A[MK]=A[MK]+C1*(SAVE1.conjugate()) # Eq. (8.D.22)
            #print 'AFTER', A[0:4]
        #print M+1, IP
        if M+1 == IP:
            #Pv.append(P)
            P=.5*P/float(N-M-1)
            Pv.append(P)
            return A, P, Pv
        else:
            Pv.append(.5*P/float(N-M-1))
            #Pv.append(P)

        #   Time update of C,D vectors and GAMMA,DELTA,LAMBDA scalars

        R1=1./(DELTA*GAMMA-(LAMBDA.real)**2-(LAMBDA.imag)**2)
        C1=(THETA*(LAMBDA.conjugate())+PSI*DELTA)*R1
        C2=(PSI*LAMBDA+THETA*GAMMA)*R1
        C3=(XI*(LAMBDA.conjugate())+THETA*DELTA)*R1
        C4=(THETA*LAMBDA+XI*GAMMA)*R1

        #print 'C1C2C3C4', C1, C2, C3, C4
        #print M, M/2+1, (M-1)/2+1
        for K in range(0, (M)//2+1):
            MK=M-K
            SAVE1=C[K].conjugate()
            SAVE2=D[K].conjugate()
            SAVE3=C[MK].conjugate()
            SAVE4=D[MK].conjugate()
            C[K]=C[K]+C1*SAVE3+C2*SAVE4                # Eq. (8.D.43)
            D[K]=D[K]+C3*SAVE3+C4*SAVE4                # Eq. (8.D.44)
            #print '-------->', K,MK, SAVE1, SAVE2, SAVE3, SAVE4
            if K != MK:
                C[MK]=C[MK]+C1*SAVE1+C2*SAVE2              # Eq. (8.D.43)
                D[MK]=D[MK]+C3*SAVE1+C4*SAVE2              # Eq. (8.D.44)
                #print 'witninloop loop ',K,MK,C[MK], D[MK]
        #print 'after loop'
        #print 'C=',C[0:2]
        #print 'D=',D[0:2]
        R2=PSI.real**2+PSI.imag**2
        R3=THETA.real**2+THETA.imag**2
        R4=XI.real**2+XI.imag**2
        R5=GAMMA-(R2*DELTA+R3*GAMMA+2.*np.real(PSI*LAMBDA*THETA.conjugate()))*R1
        R2=DELTA-(R3*DELTA+R4*GAMMA+2.*np.real(THETA*LAMBDA*XI.conjugate()))*R1
        GAMMA=R5                                     # Eq. (8.D.46)
        DELTA=R2                                     # Eq. (8.D.47)
        LAMBDA=LAMBDA+C3*PSI.conjugate()+C4*THETA.conjugate()  # Eq. (8.D.48)
        #print 'before time updaqt',R2, R3, R4, R5, LAMBDA
        if P <= 0.:
            raise ValueError('Found an invalid negative value for P.')
        if (DELTA > 0. and DELTA <= 1. and GAMMA > 0. and GAMMA <=1.):
            pass
        else:
            raise ValueError('Found an invalid DELTA or GAMMA value.')
        #   Time update of A vector; order updates of C,D vectors aand GAMMA,
        #   DELTA,LAMBDA scalars
        #print 'time update--------------'
        R1=1./P
        R2=1./(DELTA*GAMMA-LAMBDA.real**2-LAMBDA.imag**2) # Eq. (8.D.41)
        EF=X[M+1]
        EB=X[N-M-2]
        #print 'efeb', EF, EB
        #print 'A=',A[0:2]
        for K in range(0, M+1):
            #print 'efeb lloop', K,A[K], X[M-K], X[N-M+K-1]
            EF=EF+A[K]*X[M-K]                        # Eq. (8.D.1)
            EB=EB+A[K].conjugate()*X[N-M+K-1]            # Eq. (8.D.2)
        C1=EB*R1                                     # Eq. (8.D.28)
        C2=EF.conjugate()*R1                              # Eq. (8.D.29)
        C3=(EB.conjugate()*DELTA+EF*LAMBDA)*R2
        C4=(EF*GAMMA+(EB*LAMBDA).conjugate())*R2

        #print 'c1c2c3c4', C1 ,C2, C3, C4
        #print 'efeb',EF, EB

        for K in range(M, -1, -1):
            SAVE1=A[K]
            A[K]=SAVE1+C3*C[K]+C4*D[K]                 # Eq. (8.D.38)
            C[K+1]=C[K]+C1*SAVE1                       # Eq. (8.D.26)
            D[K+1]=D[K]+C2*SAVE1                       # Eq. (8.D.27)
            #print 'ACD',K, A[K], C[K+1], D[K+1]
        C[0]=C1
        D[0]=C2
        R3=EB.real**2+EB.imag**2
        R4=EF.real**2+EF.imag**2
        P=P-(R3*DELTA+R4*GAMMA+2.*np.real(EF*EB*LAMBDA))*R2  # Eq. (8.D.42)
        DELTA=DELTA-R4*R1                            # Eq. (8.D.32)
        GAMMA=GAMMA-R3*R1                            # Eq. (8.D.33)
        LAMBDA=LAMBDA+(EF*EB).conjugate()*R1                # Eq. (8.D.35)

        #print 'final'
        #print R3, R4, P, DELTA, GAMMA, LAMBDA
        #print 'Afinal=',A[0:2]
        #print 'P=',P
        if (P > 0.):
            pass
        else:
            raise ValueError("Found a invalid negative P value ")
        if (DELTA > 0. and DELTA <= 1. and  GAMMA > 0. and GAMMA <= 1.):
            pass
        else:
            raise ValueError("Found an invalid negative GAMMA or DELTA value ")
def modcovar_marple (X,IP):
    """Fast algorithm for the solution of the modified covariance least squares normal equations.

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

    :param X:    - Array of complex data samples X(1) through X(N)
    :param int IP:   - Order of linear prediction model (integer)

    :return:
        * P    - Real linear prediction variance at order IP
        * A    - Array of complex linear prediction coefficients
        * ISTAT - Integer status indicator at time of exit
              0. for normal exit (no numerical ill-conditioning)
              1. if P is not a positive value
              2. if DELTA' and GAMMA' do not lie in the range 0 to 1
              3. if P' is not a positive value
              4. if DELTA and GAMMA do not lie in the range 0 to 1


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

    .. note:: validation. results similar to test example in Marple but
        starts to differ for    ip~8. with ratio of 0.975 for ip=15 probably
        due to precision.


    :References: [Marple]_
    """
    Pv = []
    N = len(X)
    A = np.zeros(N, dtype=complex)
    D = np.zeros(N, dtype=complex)
    C = np.zeros(N, dtype=complex)
    R = np.zeros(N, dtype=complex)
    #   Initialization
    R1=0.
    for K in range(1, N-1):
        R1=R1 + 2.*(X[K].real**2 + X[K].imag**2)
    R2 = X[0].real**2 + X[0].imag**2
    R3 = X[N-1].real**2 + X[N-1].imag**2
    R4 = 1. / (R1 + 2. * (R2 + R3))
    P = R1 + R2 + R3
    DELTA = 1. - R2 * R4
    GAMMA = 1. - R3 * R4
    LAMBDA = (X[0] * X[N-1]).conjugate()*R4
    C[0] = X[N-1] * R4
    D[0] = X[0].conjugate() * R4

    M = 0
    if (IP ==0):
        P = (.5*R1+R2+R3)/float(N)
        return [], P, []

    #C   Main loop

    for M in range(0, IP):
        #print '---------------------', M
        SAVE1 = 0+0j
        for K in range(M+1, N):
            SAVE1 = SAVE1 + X[K]*X[K-M-1].conjugate()
        SAVE1 *= 2.
        R[M] = SAVE1.conjugate()
        THETA = X[N-1]*D[0]
        PSI=X[N-1]*C[0]
        XI = X[0].conjugate() * D[0]
        if M==0:
            pass
        else:
            for K in range(0, M):
                THETA=THETA+X[N-K-2]*D[K+1]    # Eq. (8.D.45)
                PSI = PSI + X[N-K-2]*C[K+1]    # Eq. (8.D.45)
                XI = XI + X[K+1].conjugate() * D[K+1] # Eq. (8.D.45)
                R[K] = R[K]-X[N-M-1] * X[N+1-M+K-1].conjugate() - X[M].conjugate() * X[M-K-1] # Eq. (8.D.37)
                SAVE1=SAVE1+R[K].conjugate()*A[M-K-1]             # Eq. (8.D.24)
                #print 'withini loop', K, THETA, PSI, XI, R[K], SAVE1
                #print 'x   -------', C[K], C[K+1]
                #print 'x   -------', D[K], D[K+1]
                #print 'x   -------', N-K-2, K+1,N-M-1,  M, M-K-1
                #print 'x   -------', X[N-K-2], X[K+1],X[N-M-1],  X[M], X[M-K-1]


        #  Order update of A vector

        C1 = -SAVE1/P
        A[M]=C1                                   # Eq. (8.D.23)
        P=P*(1.-C1.real**2-C1.imag**2)            # Eq. (8.D.25)
        #print 'C1=' , C1, 'A[M]', A[M], 'P=',P
        if M==0:
            pass
        else:
            #for k in range(0, (M-1)/2+1):
            #print 'BEFORE', A[0:4]
            for K in range(0, (M+1)//2):
                MK = M-K-1
                SAVE1=A[K]
                A[K]=SAVE1+C1*A[MK].conjugate() # Eq. (8.D.22)
                #print 'A uopdate---->',M, K, MK, A[K], A[MK], SAVE1, C1
                if (K != MK):
                    #print 'ENTERI IF', K, MK, A[MK], C1, SAVE1
                    A[MK]=A[MK]+C1*(SAVE1.conjugate()) # Eq. (8.D.22)
            #print 'AFTER', A[0:4]
        #print M+1, IP
        if M+1 == IP:
            #Pv.append(P)
            P=.5*P/float(N-M-1)
            Pv.append(P)
            return A, P, Pv
        else:
            Pv.append(.5*P/float(N-M-1))
            #Pv.append(P)

        #   Time update of C,D vectors and GAMMA,DELTA,LAMBDA scalars

        R1=1./(DELTA*GAMMA-(LAMBDA.real)**2-(LAMBDA.imag)**2)
        C1=(THETA*(LAMBDA.conjugate())+PSI*DELTA)*R1
        C2=(PSI*LAMBDA+THETA*GAMMA)*R1
        C3=(XI*(LAMBDA.conjugate())+THETA*DELTA)*R1
        C4=(THETA*LAMBDA+XI*GAMMA)*R1

        #print 'C1C2C3C4', C1, C2, C3, C4
        #print M, M/2+1, (M-1)/2+1
        for K in range(0, (M)//2+1):
            MK=M-K
            SAVE1=C[K].conjugate()
            SAVE2=D[K].conjugate()
            SAVE3=C[MK].conjugate()
            SAVE4=D[MK].conjugate()
            C[K]=C[K]+C1*SAVE3+C2*SAVE4                # Eq. (8.D.43)
            D[K]=D[K]+C3*SAVE3+C4*SAVE4                # Eq. (8.D.44)
            #print '-------->', K,MK, SAVE1, SAVE2, SAVE3, SAVE4
            if K != MK:
                C[MK]=C[MK]+C1*SAVE1+C2*SAVE2              # Eq. (8.D.43)
                D[MK]=D[MK]+C3*SAVE1+C4*SAVE2              # Eq. (8.D.44)
                #print 'witninloop loop ',K,MK,C[MK], D[MK]
        #print 'after loop'
        #print 'C=',C[0:2]
        #print 'D=',D[0:2]
        R2=PSI.real**2+PSI.imag**2
        R3=THETA.real**2+THETA.imag**2
        R4=XI.real**2+XI.imag**2
        R5=GAMMA-(R2*DELTA+R3*GAMMA+2.*np.real(PSI*LAMBDA*THETA.conjugate()))*R1
        R2=DELTA-(R3*DELTA+R4*GAMMA+2.*np.real(THETA*LAMBDA*XI.conjugate()))*R1
        GAMMA=R5                                     # Eq. (8.D.46)
        DELTA=R2                                     # Eq. (8.D.47)
        LAMBDA=LAMBDA+C3*PSI.conjugate()+C4*THETA.conjugate()  # Eq. (8.D.48)
        #print 'before time updaqt',R2, R3, R4, R5, LAMBDA
        if P <= 0.:
            raise ValueError('Found an invalid negative value for P.')
        if (DELTA > 0. and DELTA <= 1. and GAMMA > 0. and GAMMA <=1.):
            pass
        else:
            raise ValueError('Found an invalid DELTA or GAMMA value.')
        #   Time update of A vector; order updates of C,D vectors aand GAMMA,
        #   DELTA,LAMBDA scalars
        #print 'time update--------------'
        R1=1./P
        R2=1./(DELTA*GAMMA-LAMBDA.real**2-LAMBDA.imag**2) # Eq. (8.D.41)
        EF=X[M+1]
        EB=X[N-M-2]
        #print 'efeb', EF, EB
        #print 'A=',A[0:2]
        for K in range(0, M+1):
            #print 'efeb lloop', K,A[K], X[M-K], X[N-M+K-1]
            EF=EF+A[K]*X[M-K]                        # Eq. (8.D.1)
            EB=EB+A[K].conjugate()*X[N-M+K-1]            # Eq. (8.D.2)
        C1=EB*R1                                     # Eq. (8.D.28)
        C2=EF.conjugate()*R1                              # Eq. (8.D.29)
        C3=(EB.conjugate()*DELTA+EF*LAMBDA)*R2
        C4=(EF*GAMMA+(EB*LAMBDA).conjugate())*R2

        #print 'c1c2c3c4', C1 ,C2, C3, C4
        #print 'efeb',EF, EB

        for K in range(M, -1, -1):
            SAVE1=A[K]
            A[K]=SAVE1+C3*C[K]+C4*D[K]                 # Eq. (8.D.38)
            C[K+1]=C[K]+C1*SAVE1                       # Eq. (8.D.26)
            D[K+1]=D[K]+C2*SAVE1                       # Eq. (8.D.27)
            #print 'ACD',K, A[K], C[K+1], D[K+1]
        C[0]=C1
        D[0]=C2
        R3=EB.real**2+EB.imag**2
        R4=EF.real**2+EF.imag**2
        P=P-(R3*DELTA+R4*GAMMA+2.*np.real(EF*EB*LAMBDA))*R2  # Eq. (8.D.42)
        DELTA=DELTA-R4*R1                            # Eq. (8.D.32)
        GAMMA=GAMMA-R3*R1                            # Eq. (8.D.33)
        LAMBDA=LAMBDA+(EF*EB).conjugate()*R1                # Eq. (8.D.35)

        #print 'final'
        #print R3, R4, P, DELTA, GAMMA, LAMBDA
        #print 'Afinal=',A[0:2]
        #print 'P=',P
        if (P > 0.):
            pass
        else:
            raise ValueError("Found a invalid negative P value ")
        if (DELTA > 0. and DELTA <= 1. and  GAMMA > 0. and GAMMA <= 1.):
            pass
        else:
            raise ValueError("Found an invalid negative GAMMA or DELTA value ")
def modcovar_marple ( X , IP ) : Pv = [ ] N = len ( X ) A = np . zeros ( N , dtype = complex ) D = np . zeros ( N , dtype = complex ) C = np . zeros ( N , dtype = complex ) R = np . zeros ( N , dtype = complex ) #   Initialization R1 = 0. for K in range ( 1 , N - 1 ) : R1 = R1 + 2. * ( X [ K ] . real ** 2 + X [ K ] . imag ** 2 ) R2 = X [ 0 ] . real ** 2 + X [ 0 ] . imag ** 2 R3 = X [ N - 1 ] . real ** 2 + X [ N - 1 ] . imag ** 2 R4 = 1. / ( R1 + 2. * ( R2 + R3 ) ) P = R1 + R2 + R3 DELTA = 1. - R2 * R4 GAMMA = 1. - R3 * R4 LAMBDA = ( X [ 0 ] * X [ N - 1 ] ) . conjugate ( ) * R4 C [ 0 ] = X [ N - 1 ] * R4 D [ 0 ] = X [ 0 ] . conjugate ( ) * R4 M = 0 if ( IP == 0 ) : P = ( .5 * R1 + R2 + R3 ) / float ( N ) return [ ] , P , [ ] #C   Main loop for M in range ( 0 , IP ) : #print '---------------------', M SAVE1 = 0 + 0j for K in range ( M + 1 , N ) : SAVE1 = SAVE1 + X [ K ] * X [ K - M - 1 ] . conjugate ( ) SAVE1 *= 2. R [ M ] = SAVE1 . conjugate ( ) THETA = X [ N - 1 ] * D [ 0 ] PSI = X [ N - 1 ] * C [ 0 ] XI = X [ 0 ] . conjugate ( ) * D [ 0 ] if M == 0 : pass else : for K in range ( 0 , M ) : THETA = THETA + X [ N - K - 2 ] * D [ K + 1 ] # Eq. (8.D.45) PSI = PSI + X [ N - K - 2 ] * C [ K + 1 ] # Eq. (8.D.45) XI = XI + X [ K + 1 ] . conjugate ( ) * D [ K + 1 ] # Eq. (8.D.45) R [ K ] = R [ K ] - X [ N - M - 1 ] * X [ N + 1 - M + K - 1 ] . conjugate ( ) - X [ M ] . conjugate ( ) * X [ M - K - 1 ] # Eq. (8.D.37) SAVE1 = SAVE1 + R [ K ] . conjugate ( ) * A [ M - K - 1 ] # Eq. (8.D.24) #print 'withini loop', K, THETA, PSI, XI, R[K], SAVE1 #print 'x   -------', C[K], C[K+1] #print 'x   -------', D[K], D[K+1] #print 'x   -------', N-K-2, K+1,N-M-1,  M, M-K-1 #print 'x   -------', X[N-K-2], X[K+1],X[N-M-1],  X[M], X[M-K-1] #  Order update of A vector C1 = - SAVE1 / P A [ M ] = C1 # Eq. (8.D.23) P = P * ( 1. - C1 . real ** 2 - C1 . imag ** 2 ) # Eq. (8.D.25) #print 'C1=' , C1, 'A[M]', A[M], 'P=',P if M == 0 : pass else : #for k in range(0, (M-1)/2+1): #print 'BEFORE', A[0:4] for K in range ( 0 , ( M + 1 ) // 2 ) : MK = M - K - 1 SAVE1 = A [ K ] A [ K ] = SAVE1 + C1 * A [ MK ] . conjugate ( ) # Eq. (8.D.22) #print 'A uopdate---->',M, K, MK, A[K], A[MK], SAVE1, C1 if ( K != MK ) : #print 'ENTERI IF', K, MK, A[MK], C1, SAVE1 A [ MK ] = A [ MK ] + C1 * ( SAVE1 . conjugate ( ) ) # Eq. (8.D.22) #print 'AFTER', A[0:4] #print M+1, IP if M + 1 == IP : #Pv.append(P) P = .5 * P / float ( N - M - 1 ) Pv . append ( P ) return A , P , Pv else : Pv . append ( .5 * P / float ( N - M - 1 ) ) #Pv.append(P) #   Time update of C,D vectors and GAMMA,DELTA,LAMBDA scalars R1 = 1. / ( DELTA * GAMMA - ( LAMBDA . real ) ** 2 - ( LAMBDA . imag ) ** 2 ) C1 = ( THETA * ( LAMBDA . conjugate ( ) ) + PSI * DELTA ) * R1 C2 = ( PSI * LAMBDA + THETA * GAMMA ) * R1 C3 = ( XI * ( LAMBDA . conjugate ( ) ) + THETA * DELTA ) * R1 C4 = ( THETA * LAMBDA + XI * GAMMA ) * R1 #print 'C1C2C3C4', C1, C2, C3, C4 #print M, M/2+1, (M-1)/2+1 for K in range ( 0 , ( M ) // 2 + 1 ) : MK = M - K SAVE1 = C [ K ] . conjugate ( ) SAVE2 = D [ K ] . conjugate ( ) SAVE3 = C [ MK ] . conjugate ( ) SAVE4 = D [ MK ] . conjugate ( ) C [ K ] = C [ K ] + C1 * SAVE3 + C2 * SAVE4 # Eq. (8.D.43) D [ K ] = D [ K ] + C3 * SAVE3 + C4 * SAVE4 # Eq. (8.D.44) #print '-------->', K,MK, SAVE1, SAVE2, SAVE3, SAVE4 if K != MK : C [ MK ] = C [ MK ] + C1 * SAVE1 + C2 * SAVE2 # Eq. (8.D.43) D [ MK ] = D [ MK ] + C3 * SAVE1 + C4 * SAVE2 # Eq. (8.D.44) #print 'witninloop loop ',K,MK,C[MK], D[MK] #print 'after loop' #print 'C=',C[0:2] #print 'D=',D[0:2] R2 = PSI . real ** 2 + PSI . imag ** 2 R3 = THETA . real ** 2 + THETA . imag ** 2 R4 = XI . real ** 2 + XI . imag ** 2 R5 = GAMMA - ( R2 * DELTA + R3 * GAMMA + 2. * np . real ( PSI * LAMBDA * THETA . conjugate ( ) ) ) * R1 R2 = DELTA - ( R3 * DELTA + R4 * GAMMA + 2. * np . real ( THETA * LAMBDA * XI . conjugate ( ) ) ) * R1 GAMMA = R5 # Eq. (8.D.46) DELTA = R2 # Eq. (8.D.47) LAMBDA = LAMBDA + C3 * PSI . conjugate ( ) + C4 * THETA . conjugate ( ) # Eq. (8.D.48) #print 'before time updaqt',R2, R3, R4, R5, LAMBDA if P <= 0. : raise ValueError ( 'Found an invalid negative value for P.' ) if ( DELTA > 0. and DELTA <= 1. and GAMMA > 0. and GAMMA <= 1. ) : pass else : raise ValueError ( 'Found an invalid DELTA or GAMMA value.' ) #   Time update of A vector; order updates of C,D vectors aand GAMMA, #   DELTA,LAMBDA scalars #print 'time update--------------' R1 = 1. / P R2 = 1. / ( DELTA * GAMMA - LAMBDA . real ** 2 - LAMBDA . imag ** 2 ) # Eq. (8.D.41) EF = X [ M + 1 ] EB = X [ N - M - 2 ] #print 'efeb', EF, EB #print 'A=',A[0:2] for K in range ( 0 , M + 1 ) : #print 'efeb lloop', K,A[K], X[M-K], X[N-M+K-1] EF = EF + A [ K ] * X [ M - K ] # Eq. (8.D.1) EB = EB + A [ K ] . conjugate ( ) * X [ N - M + K - 1 ] # Eq. (8.D.2) C1 = EB * R1 # Eq. (8.D.28) C2 = EF . conjugate ( ) * R1 # Eq. (8.D.29) C3 = ( EB . conjugate ( ) * DELTA + EF * LAMBDA ) * R2 C4 = ( EF * GAMMA + ( EB * LAMBDA ) . conjugate ( ) ) * R2 #print 'c1c2c3c4', C1 ,C2, C3, C4 #print 'efeb',EF, EB for K in range ( M , - 1 , - 1 ) : SAVE1 = A [ K ] A [ K ] = SAVE1 + C3 * C [ K ] + C4 * D [ K ] # Eq. (8.D.38) C [ K + 1 ] = C [ K ] + C1 * SAVE1 # Eq. (8.D.26) D [ K + 1 ] = D [ K ] + C2 * SAVE1 # Eq. (8.D.27) #print 'ACD',K, A[K], C[K+1], D[K+1] C [ 0 ] = C1 D [ 0 ] = C2 R3 = EB . real ** 2 + EB . imag ** 2 R4 = EF . real ** 2 + EF . imag ** 2 P = P - ( R3 * DELTA + R4 * GAMMA + 2. * np . real ( EF * EB * LAMBDA ) ) * R2 # Eq. (8.D.42) DELTA = DELTA - R4 * R1 # Eq. (8.D.32) GAMMA = GAMMA - R3 * R1 # Eq. (8.D.33) LAMBDA = LAMBDA + ( EF * EB ) . conjugate ( ) * R1 # Eq. (8.D.35) #print 'final' #print R3, R4, P, DELTA, GAMMA, LAMBDA #print 'Afinal=',A[0:2] #print 'P=',P if ( P > 0. ) : pass else : raise ValueError ( "Found a invalid negative P value " ) if ( DELTA > 0. and DELTA <= 1. and GAMMA > 0. and GAMMA <= 1. ) : pass else : raise ValueError ( "Found an invalid negative GAMMA or DELTA value " )
Fast algorithm for the solution of the modified covariance least squares normal equations.

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

    :param X:    - Array of complex data samples X(1) through X(N)
    :param int IP:   - Order of linear prediction model (integer)

    :return:
        * P    - Real linear prediction variance at order IP
        * A    - Array of complex linear prediction coefficients
        * ISTAT - Integer status indicator at time of exit
              0. for normal exit (no numerical ill-conditioning)
              1. if P is not a positive value
              2. if DELTA' and GAMMA' do not lie in the range 0 to 1
              3. if P' is not a positive value
              4. if DELTA and GAMMA do not lie in the range 0 to 1


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

    .. note:: validation. results similar to test example in Marple but
        starts to differ for    ip~8. with ratio of 0.975 for ip=15 probably
        due to precision.


    :References: [Marple]_
Fast algorithm for the solution of the modified covariance least squares normal equations .