Commit 055407fd authored by thomas dettbarn's avatar thomas dettbarn
Browse files

initial draft of the fix point code

parent 94bfc2bd
Loading
Loading
Loading
Loading
+233 −2
Original line number Diff line number Diff line
@@ -29,7 +29,7 @@
   the United Nations Convention on Contracts on the International Sales of Goods.

*******************************************************************************************************/

#define	MYCHANGES
#include <stdint.h>
#include "options.h"
#include "prot_fx.h"
@@ -66,6 +66,36 @@ static void HouseholderReduction_fx(
    Word32 *eps_x_fx,        /* exp(eps_x_fx_e) */
    Word16 *eps_x_fx_e );

#ifdef	MYCHANGES

static void biDiagonalReductionLeft_fx(
    Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
    Word32 singularValues[MAX_OUTPUT_CHANNELS],    /* exp(singularValues_e) */
    Word32 secDiag[MAX_OUTPUT_CHANNELS],           /* exp(secDiag_e) */
    Word16 singularVectors2_e[][MAX_OUTPUT_CHANNELS],
    Word16 singularValues_e[MAX_OUTPUT_CHANNELS],
    Word16 *secDiag_e,
    const Word16 nChannelsL,  /* Q0 */
    const Word16 nChannelsC,  /* Q0 */
    const Word16 currChannel, /* Q0 */
    Word32 g, /* Q31 */
    Word16 g_e /* Q31 */
);

static void biDiagonalReductionRight_fx(
    Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
    Word32 secDiag[MAX_OUTPUT_CHANNELS],           /* exp(secDiag_e) */
    Word16 singularVectors2_e[][MAX_OUTPUT_CHANNELS],
    Word16 *secDiag_e,
    const Word16 nChannelsL,  /* Q0 */
    const Word16 nChannelsC,  /* Q0 */
    const Word16 currChannel, /* Q0 */
    Word32 *g, /* Q31 */
    Word16 *g_e /* Q31 */
);            // Q31

#else

static void biDiagonalReductionLeft_fx(
    Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
    Word32 singularValues[MAX_OUTPUT_CHANNELS],    /* exp(singularValues_e) */
@@ -93,7 +123,7 @@ static void biDiagonalReductionRight_fx(
    Word16 *sig_x_e,
    Word32 *g /* Q31 */
);            // Q31

#endif
static void singularVectorsAccumulationLeft_fx(
    Word32 singularVectors_Left[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) as Input, Q31 as output */
    Word32 singularValues[MAX_OUTPUT_CHANNELS],         /* exp(singularValues_e) */
@@ -822,11 +852,18 @@ static void HouseholderReduction_fx(
{
    Word16 nCh;
    // float g = 0.0f, sig_x = 0.0f;// to be removed
#ifdef	MYCHANGES
    Word32 g_fx = 0;
    move32();
    Word16 g_fx_e = 0;
    move16();
#else
    Word32 g_fx = 0, sig_x_fx = 0;
    move32();
    move32();
    Word16 sig_x_fx_e = 0;
    move16();
#endif

    Word16 iCh, jCh;
    Word16 singularVectors_Left_fx_e[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
@@ -842,8 +879,14 @@ static void HouseholderReduction_fx(
    /* Bidiagonal Reduction for every channel */
    FOR( nCh = 0; nCh < nChannelsC; nCh++ ) /* nChannelsC */
    {
#ifdef	MYCHANGES
        biDiagonalReductionLeft_fx( singularVectors_Left_fx, singularValues_fx, secDiag_fx, singularVectors_Left_fx_e, singularValues_fx_e, secDiag_fx_e, nChannelsL, nChannelsC, nCh, g_fx, g_fx_e );
        biDiagonalReductionRight_fx( singularVectors_Left_fx, secDiag_fx, singularVectors_Left_fx_e, secDiag_fx_e, nChannelsL, nChannelsC, nCh, &g_fx, &g_fx_e );

#else
        biDiagonalReductionLeft_fx( singularVectors_Left_fx, singularValues_fx, secDiag_fx, singularVectors_Left_fx_e, singularValues_fx_e, secDiag_fx_e, nChannelsL, nChannelsC, nCh, &sig_x_fx, &sig_x_fx_e, &g_fx );
        biDiagonalReductionRight_fx( singularVectors_Left_fx, secDiag_fx, singularVectors_Left_fx_e, secDiag_fx_e, nChannelsL, nChannelsC, nCh, &sig_x_fx, &sig_x_fx_e, &g_fx );
#endif

        Word16 L_temp_e;
        Word32 L_temp = BASOP_Util_Add_Mant32Exp( L_abs( singularValues_fx[nCh] ), singularValues_fx_e[nCh], L_abs( secDiag_fx[nCh] ), secDiag_fx_e[nCh], &L_temp_e ); /* exp(L_temp_e) */
@@ -870,7 +913,92 @@ static void HouseholderReduction_fx(
 *
 *
 *-------------------------------------------------------------------------*/
#ifdef	MYCHANGES
static void biDiagonalReductionLeft_fx(
		Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
		Word32 singularValues[MAX_OUTPUT_CHANNELS],    /* exp(singularValues_e) */
		Word32 secDiag[MAX_OUTPUT_CHANNELS],           /* exp(secDiag_e) */
		Word16 singularVectors2_e[][MAX_OUTPUT_CHANNELS],
		Word16 singularValues_e[MAX_OUTPUT_CHANNELS],
		Word16 *secDiag_e,
		const Word16 nChannelsL,  /* Q0 */
		const Word16 nChannelsC,  /* Q0 */
		const Word16 currChannel, /* Q0 */
		Word32 g, /* Q31 */
		Word16 g_e /* Q31 */
		)
{
	Word16 iCh, jCh;
	Word32 norm_x, f, r;
	Word16 norm_x_e, f_e, r_e;
	Word32 invR;
	Word16 invR_e;

	secDiag[currChannel] = g;
	secDiag_e[currChannel] = g_e;
	g=0;
	move32();
	g_e=0;
	move16();

	if (currChannel < nChannelsL )
	{
		norm_x=0;
		move32();
		norm_x_e=0;
		move16();
		for ( jCh = currChannel; jCh < nChannelsL; jCh++)
		{
			norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[jCh][currChannel], singularVectors[jCh][currChannel] ), shl( singularVectors2_e[jCh][currChannel], 1 ), &norm_x_e ); /* exp(norm_x_e) */
		}
		if (norm_x)
		{
			printf("norm: %08X  < %2d  ",norm_x,norm_x_e );
			g_e = norm_x_e;
			move32();
			g = Sqrt32(norm_x,&g_e);
			printf("sqrt: %08X < %2d\n",g,g_e);
//			g= L_shl_r(g, g_e);
			if ( singularVectors[currChannel][currChannel] >= 0)
			{
				g = L_negate(g);
			}
			move32();
			r=BASOP_Util_Add_Mant32Exp( Mpy_32_32( g, singularVectors[currChannel][currChannel] ), singularVectors2_e[currChannel][currChannel], -norm_x, norm_x_e, &r_e);
			singularVectors[currChannel][currChannel] = BASOP_Util_Add_Mant32Exp(singularVectors[currChannel][currChannel], singularVectors2_e[currChannel][currChannel], -g, 0, &singularVectors2_e[currChannel][currChannel]);
			move32();

			invR= BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx(r), &invR_e);
			FOR( iCh = currChannel + 1; iCh < nChannelsC; iCh++ ) /* nChannelsC */
			{

				norm_x = 0;
				move32();
				norm_x_e = 0;
				move16();
				FOR (jCh = currChannel; jCh< nChannelsL; jCh++)
				{
					norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[jCh][currChannel], singularVectors[jCh][iCh] ), add( singularVectors2_e[jCh][currChannel], singularVectors2_e[jCh][iCh] ), &norm_x_e ); /* exp(norm_x_e) */
				}
				f= Mpy_32_32( norm_x, invR );
				f_e = add(invR_e, norm_x);

				FOR( jCh = currChannel; jCh < nChannelsL; jCh++ ) /* nChannelsL */
				{
					singularVectors[jCh][iCh] = BASOP_Util_Add_Mant32Exp( singularVectors[jCh][iCh], singularVectors2_e[jCh][iCh], Mpy_32_32( f, singularVectors[jCh][currChannel] ), add( f_e, singularVectors2_e[jCh][currChannel] ), &singularVectors2_e[jCh][iCh] );
					move32();
				}
			}
		}
		singularValues[currChannel] = g;
		move32();
		singularValues_e[currChannel] = g_e;
		move16();
	}
	return;
}

#else
static void biDiagonalReductionLeft_fx(
    Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
    Word32 singularValues[MAX_OUTPUT_CHANNELS],    /* exp(singularValues_e) */
@@ -910,6 +1038,8 @@ static void biDiagonalReductionLeft_fx(

        FOR( jCh = idx; jCh < nChannelsL; jCh++ ) /* nChannelsL */
        {
		printf("1> singularVectors2_e[currChannel][jCh]:%d\n",singularVectors2_e[currChannel][jCh]);
		fflush(stdout);
            ( *sig_x ) = BASOP_Util_Add_Mant32Exp( *sig_x, *sig_x_e, L_abs( singularVectors[jCh][currChannel] ), singularVectors2_e[jCh][currChannel], sig_x_e ); /* exp(sig_x_e) */
        }

@@ -953,6 +1083,8 @@ static void biDiagonalReductionLeft_fx(
            r = BASOP_Util_Add_Mant32Exp( Mpy_32_32( ( *g ), singularVectors[currChannel][idx] ), singularVectors2_e[currChannel][idx], -norm_x, norm_x_e, &r_e );                                      /* exp(r_e) */
            singularVectors[currChannel][idx] = BASOP_Util_Add_Mant32Exp( singularVectors[currChannel][idx], singularVectors2_e[currChannel][idx], -( *g ), 0, &singularVectors2_e[currChannel][idx] ); /* sing_exp */
            move32();
		printf("2> singularVectors2_e[currChannel][jCh]:%d\n",singularVectors2_e[currChannel][idx]);
		fflush(stdout);

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

@@ -974,6 +1106,8 @@ static void biDiagonalReductionLeft_fx(
                {
                    singularVectors[jCh][iCh] = BASOP_Util_Add_Mant32Exp( singularVectors[jCh][iCh], singularVectors2_e[jCh][iCh], Mpy_32_32( f, singularVectors[jCh][currChannel] ), add( f_e, singularVectors2_e[jCh][currChannel] ), &singularVectors2_e[jCh][iCh] );
                    move32();
		printf("3> singularVectors2_e[currChannel][jCh]:%d\n",singularVectors2_e[jCh][idx]);
		fflush(stdout);
                }
            }

@@ -993,16 +1127,110 @@ static void biDiagonalReductionLeft_fx(
        singularValues_e[currChannel] = *sig_x_e;
        move16();
    }
	printf("-----------------\n");

    return;
}
#endif

/*-------------------------------------------------------------------------
 * biDiagonalReductionRight()
 *
 *
 *-------------------------------------------------------------------------*/
#ifdef	MYCHANGES
static void biDiagonalReductionRight_fx(
    Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
    Word32 secDiag[MAX_OUTPUT_CHANNELS],           /* exp(secDiag_exp[]) */
    Word16 singularVectors2_e[][MAX_OUTPUT_CHANNELS],
    Word16 *secDiag_exp,
    const Word16 nChannelsL,  /* Q0 */
    const Word16 nChannelsC,  /* Q0 */
    const Word16 currChannel, /* Q0 */
    Word32 *g, /* Q31 */
    Word16 *g_e /* Q31 */
)
{
	Word16 iCh, jCh, idx;
	Word32 norm_x;
	Word16 norm_x_e;
	Word32 r;
	Word16 r_e;
	Word32 invR;
	Word16 invR_e;

	*g = 0;
	move32();
	*g_e = 0;
	move16();
	IF( LT_16( currChannel, nChannelsL ) && NE_16( currChannel, sub( nChannelsC, 1 ) ) ) /* i <=m && i !=n */
	{
		idx = add(currChannel, 1);
		norm_x = 0;
		move32();
		norm_x_e = 0;
		move16();
		FOR( jCh = idx; jCh < nChannelsC; jCh++ ) /*nChannelsC */
		{
			norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[currChannel][jCh], singularVectors[currChannel][jCh] ), shl( singularVectors2_e[currChannel][jCh], 1 ), &norm_x_e ); /* exp(norm_x_e) */
		}
		if ( norm_x)
		{
			IF( GT_16( norm_x_e, 0 ) )
			{
				norm_x = MAX_32;
				move32();
				norm_x_e = 0;
				move16();
			}
			*g_e = norm_x_e;
			move32();
			*g = Sqrt32(norm_x,g_e);
			*g= L_shl_r((*g), *g_e);
			if ( singularVectors[currChannel][currChannel] >= 0)
			{
				*g = L_negate(*g);
			}
			move32();
			r=BASOP_Util_Add_Mant32Exp( Mpy_32_32( *g, singularVectors[currChannel][currChannel] ), singularVectors2_e[currChannel][currChannel], -norm_x, norm_x_e, &r_e);
			singularVectors[currChannel][currChannel] = BASOP_Util_Add_Mant32Exp(singularVectors[currChannel][currChannel], singularVectors2_e[currChannel][currChannel], -( *g ), 0, &singularVectors2_e[currChannel][currChannel]);
			move32();

			invR= BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx(r), &invR_e);
			FOR (jCh = idx; jCh < nChannelsC; jCh++)
			{
				Word16 temp_e;
				temp_e = norm_l(singularVectors[currChannel][jCh]);
				secDiag[jCh]=Mpy_32_32( L_shl(singularVectors[currChannel][jCh], temp_e), invR);
				move32();
				secDiag_exp[jCh]=add(sub(invR_e, temp_e), sub(singularVectors2_e[currChannel][jCh], r_e ) );
				move16();
			}
			FOR (iCh = currChannel + 1; iCh < nChannelsL; iCh++ )
			{
				norm_x = 0;
				move32();
				norm_x_e = 0;
				move16();
				FOR ( jCh = idx; jCh < nChannelsC ; jCh++ )
				{
					 norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[iCh][jCh], singularVectors[currChannel][jCh] ), add( singularVectors2_e[iCh][jCh], singularVectors2_e[currChannel][jCh] ), &norm_x_e );	
				}
				FOR ( jCh = idx; jCh < nChannelsC ; jCh++ )
				{
					printf("normx:%d secdiag:%d\n", norm_x_e, secDiag_exp[jCh]);
					fflush(stdout);
                    			singularVectors[iCh][jCh] = BASOP_Util_Add_Mant32Exp( singularVectors[iCh][jCh], singularVectors2_e[iCh][jCh], Mpy_32_32( norm_x, secDiag[jCh] ), add( norm_x_e, secDiag_exp[jCh] ), &singularVectors2_e[iCh][jCh] ); /* exp(sing_exp2) */
				}	
			}
		}
	}
	return;
}



#else
static void biDiagonalReductionRight_fx(
    Word32 singularVectors[][MAX_OUTPUT_CHANNELS], /* exp(singularVectors_e) */
    Word32 secDiag[MAX_OUTPUT_CHANNELS],           /* exp(secDiag_exp[]) */
@@ -1050,6 +1278,7 @@ static void biDiagonalReductionRight_fx(
            FOR( jCh = idx; jCh < nChannelsC; jCh++ ) /*nChannelsC */
            {
                temp_e = norm_l( singularVectors[currChannel][jCh] );
		printf("singularVectors2_e[currChannel][jCh]:%d\n",singularVectors2_e[currChannel][jCh]);
                singularVectors[currChannel][jCh] = Mpy_32_32( L_shl( singularVectors[currChannel][jCh], temp_e ), invVal ); /* exp(sing_exp + (singularVectors_e - sig_x_e) */
                move32();
                singularVectors2_e[currChannel][jCh] = add( sub( invVal_e, temp_e ), sub( singularVectors2_e[currChannel][jCh], *sig_x_e ) );
@@ -1123,6 +1352,8 @@ static void biDiagonalReductionRight_fx(

    return;
}
#endif


/*-------------------------------------------------------------------------
 * singularVectorsAccumulationLeft()