Commit 60ceb05f authored by thomas dettbarn's avatar thomas dettbarn
Browse files

Code cleanup in biDiagonalReductionLeft_64() and biDiagonalReductionRight_64()... 

Code cleanup in biDiagonalReductionLeft_64() and biDiagonalReductionRight_64() revealed some accuracy issues.
parent 240563be

lib_dec/ivas_svd_dec_fx.c

+18 −25
Original line number Diff line number Diff line
@@ -1005,7 +1005,6 @@ static void biDiagonalReductionLeft_64( 


#define HEADROOM_LEFT_1 1

#define HEADROOM_LEFT_2 ( HEADROOM_LEFT_1 + 1 )



    Word16 iCh, jCh;

    Word32 norm_x;

    Word16 norm_x_e;
@@ -1029,7 +1028,6 @@ static void biDiagonalReductionLeft_64( 
        }

        norm_x_e = W_norm( norm_64 );

        norm_x = W_extract_h( W_shl( norm_64, norm_x_e ) );

        norm_x_e = add( sub( ( HEADROOM_LEFT_1 + HEADROOM_LEFT_1 ), norm_x_e ), 1 );

    }

    IF( norm_x )

    {
@@ -1041,7 +1039,7 @@ static void biDiagonalReductionLeft_64( 
        Word32 r, invVal;

        Word16 r_e, invVal_e;



        ( *g_e ) = norm_x_e;

        ( *g_e ) = add( sub( ( HEADROOM_LEFT_1 + HEADROOM_LEFT_1 ), norm_x_e ), 1 );

        move16();

        ( *g ) = Sqrt32( norm_x, g_e );

        IF( GE_64( singularVectors_Left_64[currChannel][currChannel], 0 ) )
@@ -1049,16 +1047,15 @@ static void biDiagonalReductionLeft_64( 
            ( *g ) = L_negate( *g );

        }

        factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][currChannel], HEADROOM_LEFT_1 ) );

        tmp_e = sub( ( *g_e ), HEADROOM_LEFT_1 );

	tmp_e = shr(sub(norm_x_e,1),1);

        tmpmul = W_mult0_32_32( ( *g ), factor2 );

        tmpmul = W_shl( tmpmul, tmp_e );

        tmpmul = W_shr( tmpmul, tmp_e );

        r_64 = W_sub( tmpmul, norm_64 );

        r_e = W_norm( r_64 );

        r = W_extract_h( W_shl( r_64, r_e ) );

        r_e = sub( add( 1, ( HEADROOM_LEFT_1 + HEADROOM_LEFT_1 ) ), r_e );





        invVal_e = r_e;

        invVal_e = sub( add( 1, ( HEADROOM_LEFT_1 + HEADROOM_LEFT_1 ) ), r_e );

        invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, r, &invVal_e );
@@ -1067,6 +1064,7 @@ static void biDiagonalReductionLeft_64( 


        FOR( iCh = add( currChannel, 1 ); iCh < nChannelsC; iCh++ )

        {

            Word16 magic_shift;

            Word32 factor1;

            Word32 factor2;

            Word32 f; // = norm / r
@@ -1082,11 +1080,10 @@ static void biDiagonalReductionLeft_64( 
            norm_x_e = W_norm( norm_64 );

            norm_x = W_extract_h( W_shl( norm_64, norm_x_e ) );

            f = Mpy_32_32( norm_x, invVal );

	    magic_shift=31+norm_x_e-r_e;

            FOR( jCh = currChannel; jCh < nChannelsL; jCh++ )

            {

                Word16 magic_shift;

                magic_shift = add( add( norm_x_e, 23 ), r_e );

                factor1 = W_extract_l( W_shr( singularVectors_Left_64[jCh][currChannel], HEADROOM_LEFT_2 ) );

                factor1 = W_extract_l( singularVectors_Left_64[jCh][currChannel] );

                singularVectors_Left_64[jCh][iCh] = W_add( singularVectors_Left_64[jCh][iCh], W_shr( W_mult0_32_32( f, factor1 ), magic_shift ) );

            }

        }
@@ -1113,8 +1110,9 @@ static void biDiagonalReductionRight_64( 
    Word64 norm_64;

    Word16 idx;



#define HEADROOM_RIGHT_1 2

#define HEADROOM_RIGHT_1 1

#define HEADROOM_RIGHT_2 ( HEADROOM_RIGHT_1 + 1 )

#define HEADROOM_RIGHT_3 3





    ( *g ) = 0;
@@ -1134,7 +1132,6 @@ static void biDiagonalReductionRight_64( 
        }

        norm_x_e = W_norm( norm_64 );

        norm_x = W_extract_h( W_shl( norm_64, norm_x_e ) );

        norm_x_e = add( sub( ( HEADROOM_RIGHT_1 + HEADROOM_RIGHT_1 ), norm_x_e ), 1 );

        move16();



        IF( norm_x )
@@ -1143,8 +1140,6 @@ static void biDiagonalReductionRight_64( 
            Word32 factor2;

            Word16 tmp_e;

            Word64 tmpmul;

            Word32 tmp_g;

            Word16 tmp_g_e;

            Word16 magic_shift;

            Word64 r_64;

            Word32 r;
@@ -1153,20 +1148,18 @@ static void biDiagonalReductionRight_64( 
            Word32 invVal;

            Word16 invVal_e;



            tmp_g_e = norm_x_e;

            ( *g_e ) = add( sub( ( HEADROOM_RIGHT_1 + HEADROOM_RIGHT_1 ), norm_x_e ), 1 );

            move16();

            tmp_g = Sqrt32( norm_x, &tmp_g_e );

            ( *g ) = Sqrt32( norm_x, g_e );

            IF( GE_64( singularVectors_Left_64[currChannel][idx], 0 ) )

            {

                tmp_g = L_negate( tmp_g );

                ( *g ) = L_negate( *g );

            }

            *g = tmp_g;

            *g_e = tmp_g_e;

            move32();

            move16();

            factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][idx], HEADROOM_RIGHT_1 ) );

            tmp_e = sub( tmp_g_e, HEADROOM_RIGHT_1 );

            tmpmul = W_mult0_32_32( tmp_g, factor2 );

            tmp_e = sub( *g_e, HEADROOM_RIGHT_1 );

            tmpmul = W_mult0_32_32( *g, factor2 );

            tmpmul = W_shl( tmpmul, tmp_e );

            r_64 = W_sub( tmpmul, norm_64 );

            r_e = W_norm( r_64 );
@@ -1176,8 +1169,8 @@ static void biDiagonalReductionRight_64( 
            move16();

            invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx( r ), &invVal_e );



            magic_shift = 32 - tmp_g_e;

            singularVectors_Left_64[currChannel][idx] = W_sub( singularVectors_Left_64[currChannel][idx], W_shr( W_deposit32_h( tmp_g ), magic_shift ) ); // here, the exponent goes up

            magic_shift = 32 - *g_e;

            singularVectors_Left_64[currChannel][idx] = W_sub( singularVectors_Left_64[currChannel][idx], W_shr( W_deposit32_h( *g ), magic_shift ) ); // here, the exponent goes up



            FOR( iCh = idx; iCh < nChannelsL; iCh++ )

            {
@@ -1194,11 +1187,11 @@ static void biDiagonalReductionRight_64( 
                norm_x_e = W_norm( norm_64 );

                norm_x = W_extract_h( W_shl( norm_64, norm_x_e ) );

                f = Mpy_32_32( norm_x, invVal );

                magic_shift = 25 + norm_x_e - r_e; // FIXME: Why does this work?

                magic_shift = 25 + norm_x_e - r_e; // headroom 3 FIXME: Why does this work?



                FOR( jCh = idx; jCh < nChannelsC; jCh++ )

                {

                    factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][jCh], HEADROOM_RIGHT_2 ) );

                    factor2 = W_extract_l( W_shr( singularVectors_Left_64[currChannel][jCh], HEADROOM_RIGHT_3 ) );

                    singularVectors_Left_64[iCh][jCh] = W_add( singularVectors_Left_64[iCh][jCh], W_shr( W_mult0_32_32( f, factor2 ), magic_shift ) );

                }

            }