lpc.c

Go to the documentation of this file.

/*
 * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
 * Universitaet Berlin. See the accompanying file "COPYRIGHT" for
 * details. THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
 */

/* $Header$ */

#include <stdio.h>
#include <assert.h>

#include "private.h"

#include "gsm.h"
#include "proto.h"

#ifdef K6OPT
#include "k6opt.h"
#endif

#undef P

/*
 * 4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
 */

/* 4.2.4 */


static void Autocorrelation P2((s, L_ACF),
 word * s, /* [0..159] IN/OUT */
 longword * L_ACF) /* [0..8] OUT */
/*
 * The goal is to compute the array L_ACF[k]. The signal s[i] must
 * be scaled in order to avoid an overflow situation.
 */
{
#ifndef K6OPT
 register int k, i;
 word temp;
#endif

 word smax, scalauto;

#ifdef USE_FLOAT_MUL
 float float_s[160];
#endif

 /* Dynamic scaling of the array s[0..159]
 */

 /* Search for the maximum.
 */
#ifndef K6OPT
 smax = 0;
 for (k = 0; k <= 159; k++) {
 temp = GSM_ABS( s[k] );
 if (temp > smax) smax = temp;
 }
#else
 {
 longword lmax;
 lmax = k6maxmin(s,160,NULL);
 smax = (lmax > MAX_WORD) ? MAX_WORD : lmax;
 }
#endif
 /* Computation of the scaling factor.
 */
 if (smax == 0) scalauto = 0;
 else {
 assert(smax > 0);
 scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
 }

 /* Scaling of the array s[0...159]
 */

 if (scalauto > 0) {
# ifndef K6OPT

# ifdef USE_FLOAT_MUL
# define SCALE(n) \
 case n: for (k = 0; k <= 159; k++) \
 float_s[k] = (float) \
 (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
 break;
# else
# define SCALE(n) \
 case n: for (k = 0; k <= 159; k++) \
 s[k] = (word)GSM_MULT_R( s[k], 16384 >> (n-1) );\
 break;
# endif /* USE_FLOAT_MUL */

 switch (scalauto) {
 SCALE(1)
 SCALE(2)
 SCALE(3)
 SCALE(4)
 }
# undef SCALE

# else /* K6OPT */
 k6vsraw(s,160,scalauto);
# endif
 }
# ifdef USE_FLOAT_MUL
 else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
# endif

 /* Compute the L_ACF[..].
 */
#ifndef K6OPT
 {
# ifdef USE_FLOAT_MUL
 register float * sp = float_s;
 register float sl = *sp;

# define STEP(k) L_ACF[k] += (longword)(sl * sp[ -(k) ]);
# else
 word * sp = s;
 word sl = *sp;

# define STEP(k) L_ACF[k] += ((longword)sl * sp[ -(k) ]);
# endif

# define NEXTI sl = *++sp


 for (k = 9; k--; L_ACF[k] = 0) ;

 STEP (0);
 NEXTI;
 STEP(0); STEP(1);
 NEXTI;
 STEP(0); STEP(1); STEP(2);
 NEXTI;
 STEP(0); STEP(1); STEP(2); STEP(3);
 NEXTI;
 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
 NEXTI;
 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
 NEXTI;
 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
 NEXTI;
 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);

 for (i = 8; i <= 159; i++) {

 NEXTI;

 STEP(0);
 STEP(1); STEP(2); STEP(3); STEP(4);
 STEP(5); STEP(6); STEP(7); STEP(8);
 }

 for (k = 9; k--; L_ACF[k] <<= 1) ;

 }

#else
 {
 int k;
 for (k=0; k<9; k++) {
 L_ACF[k] = 2*k6iprod(s,s+k,160-k);
 }
 }
#endif
 /* Rescaling of the array s[0..159]
 */
 if (scalauto > 0) {
 assert(scalauto <= 4);
#ifndef K6OPT
 for (k = 160; k--; *s++ <<= scalauto) ;
# else /* K6OPT */
 k6vsllw(s,160,scalauto);
# endif
 }
}

#if defined(USE_FLOAT_MUL) && defined(FAST)

static void Fast_Autocorrelation P2((s, L_ACF),
 word * s, /* [0..159] IN/OUT */
 longword * L_ACF) /* [0..8] OUT */
{
 register int k, i;
 float f_L_ACF[9];
 float scale;

 float s_f[160];
 register float *sf = s_f;

 for (i = 0; i < 160; ++i) sf[i] = s[i];
 for (k = 0; k <= 8; k++) {
 register float L_temp2 = 0;
 register float *sfl = sf - k;
 for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
 f_L_ACF[k] = L_temp2;
 }
 scale = MAX_LONGWORD / f_L_ACF[0];

 for (k = 0; k <= 8; k++) {
 L_ACF[k] = f_L_ACF[k] * scale;
 }
}
#endif /* defined (USE_FLOAT_MUL) && defined (FAST) */

/* 4.2.5 */

static void Reflection_coefficients P2( (L_ACF, r),
 longword * L_ACF, /* 0...8 IN */
 register word * r /* 0...7 OUT */
)
{
 register int i, m, n;
 register word temp;
 word ACF[9]; /* 0..8 */
 word P[ 9]; /* 0..8 */
 word K[ 9]; /* 2..8 */

 /* Schur recursion with 16 bits arithmetic.
 */

 if (L_ACF[0] == 0) {
 for (i = 8; i--; *r++ = 0) ;
 return;
 }

 assert( L_ACF[0] != 0 );
 temp = gsm_norm( L_ACF[0] );

 assert(temp >= 0 && temp < 32);

 /* ? overflow ? */
 for (i = 0; i <= 8; i++) ACF[i] = (word)SASR( L_ACF[i] << temp, 16 );

 /* Initialize array P[..] and K[..] for the recursion.
 */

 for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
 for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];

 /* Compute reflection coefficients
 */
 for (n = 1; n <= 8; n++, r++) {

 temp = P[1];
 temp = GSM_ABS(temp);
 if (P[0] < temp) {
 for (i = n; i <= 8; i++) *r++ = 0;
 return;
 }

 *r = gsm_div( temp, P[0] );

 assert(*r >= 0);
 if (P[1] > 0) *r = -*r; /* r[n] = sub(0, r[n]) */
 assert (*r != MIN_WORD);
 if (n == 8) return;

 /* Schur recursion
 */
 temp = (word)GSM_MULT_R( P[1], *r );
 P[0] = GSM_ADD( P[0], temp );

 for (m = 1; m <= 8 - n; m++) {
 temp = (word)GSM_MULT_R( K[ m ], *r );
 P[m] = GSM_ADD( P[ m+1 ], temp );

 temp = (word)GSM_MULT_R( P[ m+1 ], *r );
 K[m] = GSM_ADD( K[ m ], temp );
 }
 }
}

/* 4.2.6 */

static void Transformation_to_Log_Area_Ratios P1((r),
 register word * r /* 0..7 IN/OUT */
)
/*
 * The following scaling for r[..] and LAR[..] has been used:
 *
 * r[..] = integer( real_r[..]*32768. ); -1 <= real_r < 1.
 * LAR[..] = integer( real_LAR[..] * 16384 );
 * with -1.625 <= real_LAR <= 1.625
 */
{
 register word temp;
 register int i;


 /* Computation of the LAR[0..7] from the r[0..7]
 */
 for (i = 1; i <= 8; i++, r++) {

 temp = *r;
 temp = GSM_ABS(temp);
 assert(temp >= 0);

 if (temp < 22118) {
 temp >>= 1;
 } else if (temp < 31130) {
 assert( temp >= 11059 );
 temp -= 11059;
 } else {
 assert( temp >= 26112 );
 temp -= 26112;
 temp <<= 2;
 }

 *r = *r < 0 ? -temp : temp;
 assert( *r != MIN_WORD );
 }
}

/* 4.2.7 */

static void Quantization_and_coding P1((LAR),
 register word * LAR /* [0..7] IN/OUT */
)
{
 register word temp;


 /* This procedure needs four tables; the following equations
 * give the optimum scaling for the constants:
 *
 * A[0..7] = integer( real_A[0..7] * 1024 )
 * B[0..7] = integer( real_B[0..7] * 512 )
 * MAC[0..7] = maximum of the LARc[0..7]
 * MIC[0..7] = minimum of the LARc[0..7]
 */

# undef STEP
# define STEP( A, B, MAC, MIC ) \
 temp = (word)GSM_MULT( A, *LAR ); \
 temp = GSM_ADD( temp, B ); \
 temp = GSM_ADD( temp, 256 ); \
 temp = (word)SASR( temp, 9 ); \
 *LAR = temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
 LAR++;

 STEP( 20480, 0, 31, -32 );
 STEP( 20480, 0, 31, -32 );
 STEP( 20480, 2048, 15, -16 );
 STEP( 20480, -2560, 15, -16 );

 STEP( 13964, 94, 7, -8 );
 STEP( 15360, -1792, 7, -8 );
 STEP( 8534, -341, 3, -4 );
 STEP( 9036, -1144, 3, -4 );

# undef STEP
}

void Gsm_LPC_Analysis P3((S, s,LARc),
 struct gsm_state *S,
 word * s, /* 0..159 signals IN/OUT */
 word * LARc) /* 0..7 LARc's OUT */
{
 longword L_ACF[9];

#if defined(USE_FLOAT_MUL) && defined(FAST)
 if (S->fast) Fast_Autocorrelation (s, L_ACF );
 else
#endif
 Autocorrelation (s, L_ACF );
 Reflection_coefficients (L_ACF, LARc );
 Transformation_to_Log_Area_Ratios (LARc);
 Quantization_and_coding (LARc);
}