MythTV: lpc.c Source File

00001 /*
00002  * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
00003  * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
00004  * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
00005  */
00006 
00007 /* $Header$ */
00008 
00009 #include <stdio.h>
00010 #include <assert.h>
00011 
00012 #include "private.h"
00013 
00014 #include "gsm.h"
00015 #include "proto.h"
00016 
00017 #ifdef K6OPT
00018 #include "k6opt.h"
00019 #endif
00020 
00021 #undef  P
00022 
00023 /*
00024  *  4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
00025  */
00026 
00027 /* 4.2.4 */
00028 
00029 
00030 static void Autocorrelation P2((s, L_ACF),
00031         word     * s,           /* [0..159]     IN/OUT  */
00032         longword * L_ACF)       /* [0..8]       OUT     */
00033 /*
00034  *  The goal is to compute the array L_ACF[k].  The signal s[i] must
00035  *  be scaled in order to avoid an overflow situation.
00036  */
00037 {
00038         register int    k, i;
00039 
00040         word            temp, smax, scalauto;
00041 
00042 #ifdef  USE_FLOAT_MUL
00043         float           float_s[160];
00044 #endif
00045 
00046         /*  Dynamic scaling of the array  s[0..159]
00047          */
00048 
00049         /*  Search for the maximum.
00050          */
00051 #ifndef K6OPT
00052         smax = 0;
00053         for (k = 0; k <= 159; k++) {
00054                 temp = GSM_ABS( s[k] );
00055                 if (temp > smax) smax = temp;
00056         }
00057 #else
00058         {
00059                 longword lmax;
00060                 lmax = k6maxmin(s,160,NULL);
00061                 smax = (lmax > MAX_WORD) ? MAX_WORD : lmax;
00062         }
00063 #endif
00064         /*  Computation of the scaling factor.
00065          */
00066         if (smax == 0) scalauto = 0;
00067         else {
00068                 assert(smax > 0);
00069                 scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
00070         }
00071 
00072         /*  Scaling of the array s[0...159]
00073          */
00074 
00075         if (scalauto > 0) {
00076 #       ifndef K6OPT
00077 
00078 # ifdef USE_FLOAT_MUL
00079 #   define SCALE(n)     \
00080         case n: for (k = 0; k <= 159; k++) \
00081                         float_s[k] = (float)    \
00082                                 (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
00083                 break;
00084 # else 
00085 #   define SCALE(n)     \
00086         case n: for (k = 0; k <= 159; k++) \
00087                         s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
00088                 break;
00089 # endif /* USE_FLOAT_MUL */
00090 
00091                 switch (scalauto) {
00092                 SCALE(1)
00093                 SCALE(2)
00094                 SCALE(3)
00095                 SCALE(4)
00096                 }
00097 # undef SCALE
00098 
00099 #       else /* K6OPT */
00100                 k6vsraw(s,160,scalauto);
00101 #       endif
00102         }
00103 # ifdef USE_FLOAT_MUL
00104         else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
00105 # endif
00106 
00107         /*  Compute the L_ACF[..].
00108          */
00109 #ifndef K6OPT
00110         {
00111 # ifdef USE_FLOAT_MUL
00112                 register float * sp = float_s;
00113                 register float   sl = *sp;
00114 
00115 #               define STEP(k)   L_ACF[k] += (longword)(sl * sp[ -(k) ]);
00116 # else
00117                 word  * sp = s;
00118                 word    sl = *sp;
00119 
00120 #               define STEP(k)   L_ACF[k] += ((longword)sl * sp[ -(k) ]);
00121 # endif
00122 
00123 #       define NEXTI     sl = *++sp
00124 
00125 
00126         for (k = 9; k--; L_ACF[k] = 0) ;
00127 
00128         STEP (0);
00129         NEXTI;
00130         STEP(0); STEP(1);
00131         NEXTI;
00132         STEP(0); STEP(1); STEP(2);
00133         NEXTI;
00134         STEP(0); STEP(1); STEP(2); STEP(3);
00135         NEXTI;
00136         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
00137         NEXTI;
00138         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
00139         NEXTI;
00140         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
00141         NEXTI;
00142         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
00143 
00144         for (i = 8; i <= 159; i++) {
00145 
00146                 NEXTI;
00147 
00148                 STEP(0);
00149                 STEP(1); STEP(2); STEP(3); STEP(4);
00150                 STEP(5); STEP(6); STEP(7); STEP(8);
00151         }
00152 
00153         for (k = 9; k--; L_ACF[k] <<= 1) ; 
00154 
00155         }
00156 
00157 #else
00158         {
00159                 int k;
00160                 for (k=0; k<9; k++) {
00161                         L_ACF[k] = 2*k6iprod(s,s+k,160-k);
00162                 }
00163         }
00164 #endif
00165         /*   Rescaling of the array s[0..159]
00166          */
00167         if (scalauto > 0) {
00168                 assert(scalauto <= 4); 
00169 #ifndef K6OPT
00170                 for (k = 160; k--; *s++ <<= scalauto) ;
00171 #       else /* K6OPT */
00172                 k6vsllw(s,160,scalauto);
00173 #       endif
00174         }
00175 }
00176 
00177 #if defined(USE_FLOAT_MUL) && defined(FAST)
00178 
00179 static void Fast_Autocorrelation P2((s, L_ACF),
00180         word * s,               /* [0..159]     IN/OUT  */
00181         longword * L_ACF)       /* [0..8]       OUT     */
00182 {
00183         register int    k, i;
00184         float f_L_ACF[9];
00185         float scale;
00186 
00187         float          s_f[160];
00188         register float *sf = s_f;
00189 
00190         for (i = 0; i < 160; ++i) sf[i] = s[i];
00191         for (k = 0; k <= 8; k++) {
00192                 register float L_temp2 = 0;
00193                 register float *sfl = sf - k;
00194                 for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
00195                 f_L_ACF[k] = L_temp2;
00196         }
00197         scale = MAX_LONGWORD / f_L_ACF[0];
00198 
00199         for (k = 0; k <= 8; k++) {
00200                 L_ACF[k] = f_L_ACF[k] * scale;
00201         }
00202 }
00203 #endif  /* defined (USE_FLOAT_MUL) && defined (FAST) */
00204 
00205 /* 4.2.5 */
00206 
00207 static void Reflection_coefficients P2( (L_ACF, r),
00208         longword        * L_ACF,                /* 0...8        IN      */
00209         register word   * r                     /* 0...7        OUT     */
00210 )
00211 {
00212         register int    i, m, n;
00213         register word   temp;
00214         register longword ltmp;
00215         word            ACF[9]; /* 0..8 */
00216         word            P[  9]; /* 0..8 */
00217         word            K[  9]; /* 2..8 */
00218 
00219         /*  Schur recursion with 16 bits arithmetic.
00220          */
00221 
00222         if (L_ACF[0] == 0) {
00223                 for (i = 8; i--; *r++ = 0) ;
00224                 return;
00225         }
00226 
00227         assert( L_ACF[0] != 0 );
00228         temp = gsm_norm( L_ACF[0] );
00229 
00230         assert(temp >= 0 && temp < 32);
00231 
00232         /* ? overflow ? */
00233         for (i = 0; i <= 8; i++) ACF[i] = SASR( L_ACF[i] << temp, 16 );
00234 
00235         /*   Initialize array P[..] and K[..] for the recursion.
00236          */
00237 
00238         for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
00239         for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
00240 
00241         /*   Compute reflection coefficients
00242          */
00243         for (n = 1; n <= 8; n++, r++) {
00244 
00245                 temp = P[1];
00246                 temp = GSM_ABS(temp);
00247                 if (P[0] < temp) {
00248                         for (i = n; i <= 8; i++) *r++ = 0;
00249                         return;
00250                 }
00251 
00252                 *r = gsm_div( temp, P[0] );
00253 
00254                 assert(*r >= 0);
00255                 if (P[1] > 0) *r = -*r;         /* r[n] = sub(0, r[n]) */
00256                 assert (*r != MIN_WORD);
00257                 if (n == 8) return; 
00258 
00259                 /*  Schur recursion
00260                  */
00261                 temp = GSM_MULT_R( P[1], *r );
00262                 P[0] = GSM_ADD( P[0], temp );
00263 
00264                 for (m = 1; m <= 8 - n; m++) {
00265                         temp     = GSM_MULT_R( K[ m   ],    *r );
00266                         P[m]     = GSM_ADD(    P[ m+1 ],  temp );
00267 
00268                         temp     = GSM_MULT_R( P[ m+1 ],    *r );
00269                         K[m]     = GSM_ADD(    K[ m   ],  temp );
00270                 }
00271         }
00272 }
00273 
00274 /* 4.2.6 */
00275 
00276 static void Transformation_to_Log_Area_Ratios P1((r),
00277         register word   * r                     /* 0..7    IN/OUT */
00278 )
00279 /*
00280  *  The following scaling for r[..] and LAR[..] has been used:
00281  *
00282  *  r[..]   = integer( real_r[..]*32768. ); -1 <= real_r < 1.
00283  *  LAR[..] = integer( real_LAR[..] * 16384 );
00284  *  with -1.625 <= real_LAR <= 1.625
00285  */
00286 {
00287         register word   temp;
00288         register int    i;
00289 
00290 
00291         /* Computation of the LAR[0..7] from the r[0..7]
00292          */
00293         for (i = 1; i <= 8; i++, r++) {
00294 
00295                 temp = *r;
00296                 temp = GSM_ABS(temp);
00297                 assert(temp >= 0);
00298 
00299                 if (temp < 22118) {
00300                         temp >>= 1;
00301                 } else if (temp < 31130) {
00302                         assert( temp >= 11059 );
00303                         temp -= 11059;
00304                 } else {
00305                         assert( temp >= 26112 );
00306                         temp -= 26112;
00307                         temp <<= 2;
00308                 }
00309 
00310                 *r = *r < 0 ? -temp : temp;
00311                 assert( *r != MIN_WORD );
00312         }
00313 }
00314 
00315 /* 4.2.7 */
00316 
00317 static void Quantization_and_coding P1((LAR),
00318         register word * LAR     /* [0..7]       IN/OUT  */
00319 )
00320 {
00321         register word   temp;
00322         longword        ltmp;
00323 
00324 
00325         /*  This procedure needs four tables; the following equations
00326          *  give the optimum scaling for the constants:
00327          *  
00328          *  A[0..7] = integer( real_A[0..7] * 1024 )
00329          *  B[0..7] = integer( real_B[0..7] *  512 )
00330          *  MAC[0..7] = maximum of the LARc[0..7]
00331          *  MIC[0..7] = minimum of the LARc[0..7]
00332          */
00333 
00334 #       undef STEP
00335 #       define  STEP( A, B, MAC, MIC )          \
00336                 temp = GSM_MULT( A,   *LAR );   \
00337                 temp = GSM_ADD(  temp,   B );   \
00338                 temp = GSM_ADD(  temp, 256 );   \
00339                 temp = SASR(     temp,   9 );   \
00340                 *LAR  =  temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
00341                 LAR++;
00342 
00343         STEP(  20480,     0,  31, -32 );
00344         STEP(  20480,     0,  31, -32 );
00345         STEP(  20480,  2048,  15, -16 );
00346         STEP(  20480, -2560,  15, -16 );
00347 
00348         STEP(  13964,    94,   7,  -8 );
00349         STEP(  15360, -1792,   7,  -8 );
00350         STEP(   8534,  -341,   3,  -4 );
00351         STEP(   9036, -1144,   3,  -4 );
00352 
00353 #       undef   STEP
00354 }
00355 
00356 void Gsm_LPC_Analysis P3((S, s,LARc),
00357         struct gsm_state *S,
00358         word             * s,           /* 0..159 signals       IN/OUT  */
00359         word             * LARc)        /* 0..7   LARc's        OUT     */
00360 {
00361         longword        L_ACF[9];
00362 
00363 #if defined(USE_FLOAT_MUL) && defined(FAST)
00364         if (S->fast) Fast_Autocorrelation (s,     L_ACF );
00365         else
00366 #endif
00367         Autocorrelation                   (s,     L_ACF );
00368         Reflection_coefficients           (L_ACF, LARc  );
00369         Transformation_to_Log_Area_Ratios (LARc);
00370         Quantization_and_coding           (LARc);
00371 }