home *** CD-ROM | disk | FTP | other *** search
/ Dream 52 / Amiga_Dream_52.iso / Amiga / Workbench / Archivers / gsmPPC.lha / gsmPPC / src / lpc.c < prev    next >
C/C++ Source or Header  |  1998-04-27  |  7KB  |  342 lines
  1. /*
  2.  * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
  3.  * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
  4.  * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
  5.  */
  6.  
  7. /* $Header: /tmp_amd/presto/export/kbs/jutta/src/gsm/RCS/lpc.c,v 1.5 1994/12/30 23:14:54 jutta Exp $ */
  8.  
  9. #include <stdio.h>
  10. #include <assert.h>
  11.  
  12. #include "private.h"
  13.  
  14. #include "gsm.h"
  15. #include "proto.h"
  16.  
  17. #undef    P
  18.  
  19. /*
  20.  *  4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
  21.  */
  22.  
  23. /* 4.2.4 */
  24.  
  25.  
  26. static void Autocorrelation P2((s, L_ACF),
  27.     word     * s,        /* [0..159]    IN/OUT  */
  28.      longword * L_ACF)    /* [0..8]    OUT     */
  29. /*
  30.  *  The goal is to compute the array L_ACF[k].  The signal s[i] must
  31.  *  be scaled in order to avoid an overflow situation.
  32.  */
  33. {
  34.     register int    k, i;
  35.  
  36.     word        temp, smax, scalauto;
  37.  
  38. #ifdef    USE_FLOAT_MUL
  39.     float        float_s[160];
  40. #endif
  41.  
  42.     /*  Dynamic scaling of the array  s[0..159]
  43.      */
  44.  
  45.     /*  Search for the maximum.
  46.      */
  47.     smax = 0;
  48.     for (k = 0; k <= 159; k++) {
  49.         temp = GSM_ABS( s[k] );
  50.         if (temp > smax) smax = temp;
  51.     }
  52.  
  53.     /*  Computation of the scaling factor.
  54.      */
  55.     if (smax == 0) scalauto = 0;
  56.     else {
  57.         assert(smax > 0);
  58.         scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
  59.     }
  60.  
  61.     /*  Scaling of the array s[0...159]
  62.      */
  63.  
  64.     if (scalauto > 0) {
  65.  
  66. # ifdef USE_FLOAT_MUL
  67. #   define SCALE(n)    \
  68.     case n: for (k = 0; k <= 159; k++) \
  69.             float_s[k] = (float)    \
  70.                 (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
  71.         break;
  72. # else 
  73. #   define SCALE(n)    \
  74.     case n: for (k = 0; k <= 159; k++) \
  75.             s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
  76.         break;
  77. # endif /* USE_FLOAT_MUL */
  78.  
  79.         switch (scalauto) {
  80.         SCALE(1)
  81.         SCALE(2)
  82.         SCALE(3)
  83.         SCALE(4)
  84.         }
  85. # undef    SCALE
  86.     }
  87. # ifdef    USE_FLOAT_MUL
  88.     else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
  89. # endif
  90.  
  91.     /*  Compute the L_ACF[..].
  92.      */
  93.     {
  94. # ifdef    USE_FLOAT_MUL
  95.         register float * sp = float_s;
  96.         register float   sl = *sp;
  97.  
  98. #        define STEP(k)     L_ACF[k] += (longword)(sl * sp[ -(k) ]);
  99. # else
  100.         word  * sp = s;
  101.         word    sl = *sp;
  102.  
  103. #        define STEP(k)     L_ACF[k] += ((longword)sl * sp[ -(k) ]);
  104. # endif
  105.  
  106. #    define NEXTI     sl = *++sp
  107.  
  108.  
  109.     for (k = 9; k--; L_ACF[k] = 0) ;
  110.  
  111.     STEP (0);
  112.     NEXTI;
  113.     STEP(0); STEP(1);
  114.     NEXTI;
  115.     STEP(0); STEP(1); STEP(2);
  116.     NEXTI;
  117.     STEP(0); STEP(1); STEP(2); STEP(3);
  118.     NEXTI;
  119.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
  120.     NEXTI;
  121.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
  122.     NEXTI;
  123.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
  124.     NEXTI;
  125.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
  126.  
  127.     for (i = 8; i <= 159; i++) {
  128.  
  129.         NEXTI;
  130.  
  131.         STEP(0);
  132.         STEP(1); STEP(2); STEP(3); STEP(4);
  133.         STEP(5); STEP(6); STEP(7); STEP(8);
  134.     }
  135.  
  136.     for (k = 9; k--; L_ACF[k] <<= 1) ; 
  137.  
  138.     }
  139.     /*   Rescaling of the array s[0..159]
  140.      */
  141.     if (scalauto > 0) {
  142.         assert(scalauto <= 4); 
  143.         for (k = 160; k--; *s++ <<= scalauto) ;
  144.     }
  145. }
  146.  
  147. #if defined(USE_FLOAT_MUL) && defined(FAST)
  148.  
  149. static void Fast_Autocorrelation P2((s, L_ACF),
  150.     word * s,        /* [0..159]    IN/OUT  */
  151.      longword * L_ACF)    /* [0..8]    OUT     */
  152. {
  153.     register int    k, i;
  154.     float f_L_ACF[9];
  155.     float scale;
  156.  
  157.     float          s_f[160];
  158.     register float *sf = s_f;
  159.  
  160.     for (i = 0; i < 160; ++i) sf[i] = s[i];
  161.     for (k = 0; k <= 8; k++) {
  162.         register float L_temp2 = 0;
  163.         register float *sfl = sf - k;
  164.         for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
  165.         f_L_ACF[k] = L_temp2;
  166.     }
  167.     scale = MAX_LONGWORD / f_L_ACF[0];
  168.  
  169.     for (k = 0; k <= 8; k++) {
  170.         L_ACF[k] = f_L_ACF[k] * scale;
  171.     }
  172. }
  173. #endif    /* defined (USE_FLOAT_MUL) && defined (FAST) */
  174.  
  175. /* 4.2.5 */
  176.  
  177. static void Reflection_coefficients P2( (L_ACF, r),
  178.     longword    * L_ACF,        /* 0...8    IN    */
  179.     register word    * r            /* 0...7    OUT     */
  180. )
  181. {
  182.     register int    i, m, n;
  183.     register word    temp;
  184.     register longword ltmp;
  185.     word        ACF[9];    /* 0..8 */
  186.     word        P[  9];    /* 0..8 */
  187.     word        K[  9]; /* 2..8 */
  188.  
  189.     /*  Schur recursion with 16 bits arithmetic.
  190.      */
  191.  
  192.     if (L_ACF[0] == 0) {
  193.         for (i = 8; i--; *r++ = 0) ;
  194.         return;
  195.     }
  196.  
  197.     assert( L_ACF[0] != 0 );
  198.     temp = gsm_norm( L_ACF[0] );
  199.  
  200.     assert(temp >= 0 && temp < 32);
  201.  
  202.     /* ? overflow ? */
  203.     for (i = 0; i <= 8; i++) ACF[i] = SASR( L_ACF[i] << temp, 16 );
  204.  
  205.     /*   Initialize array P[..] and K[..] for the recursion.
  206.      */
  207.  
  208.     for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
  209.     for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
  210.  
  211.     /*   Compute reflection coefficients
  212.      */
  213.     for (n = 1; n <= 8; n++, r++) {
  214.  
  215.         temp = P[1];
  216.         temp = GSM_ABS(temp);
  217.         if (P[0] < temp) {
  218.             for (i = n; i <= 8; i++) *r++ = 0;
  219.             return;
  220.         }
  221.  
  222.         *r = gsm_div( temp, P[0] );
  223.  
  224.         assert(*r >= 0);
  225.         if (P[1] > 0) *r = -*r;        /* r[n] = sub(0, r[n]) */
  226.         assert (*r != MIN_WORD);
  227.         if (n == 8) return; 
  228.  
  229.         /*  Schur recursion
  230.          */
  231.         temp = GSM_MULT_R( P[1], *r );
  232.         P[0] = GSM_ADD( P[0], temp );
  233.  
  234.         for (m = 1; m <= 8 - n; m++) {
  235.             temp     = GSM_MULT_R( K[ m   ],    *r );
  236.             P[m]     = GSM_ADD(    P[ m+1 ],  temp );
  237.  
  238.             temp     = GSM_MULT_R( P[ m+1 ],    *r );
  239.             K[m]     = GSM_ADD(    K[ m   ],  temp );
  240.         }
  241.     }
  242. }
  243.  
  244. /* 4.2.6 */
  245.  
  246. static void Transformation_to_Log_Area_Ratios P1((r),
  247.     register word    * r             /* 0..7       IN/OUT */
  248. )
  249. /*
  250.  *  The following scaling for r[..] and LAR[..] has been used:
  251.  *
  252.  *  r[..]   = integer( real_r[..]*32768. ); -1 <= real_r < 1.
  253.  *  LAR[..] = integer( real_LAR[..] * 16384 );
  254.  *  with -1.625 <= real_LAR <= 1.625
  255.  */
  256. {
  257.     register word    temp;
  258.     register int    i;
  259.  
  260.  
  261.     /* Computation of the LAR[0..7] from the r[0..7]
  262.      */
  263.     for (i = 1; i <= 8; i++, r++) {
  264.  
  265.         temp = *r;
  266.         temp = GSM_ABS(temp);
  267.         assert(temp >= 0);
  268.  
  269.         if (temp < 22118) {
  270.             temp >>= 1;
  271.         } else if (temp < 31130) {
  272.             assert( temp >= 11059 );
  273.             temp -= 11059;
  274.         } else {
  275.             assert( temp >= 26112 );
  276.             temp -= 26112;
  277.             temp <<= 2;
  278.         }
  279.  
  280.         *r = *r < 0 ? -temp : temp;
  281.         assert( *r != MIN_WORD );
  282.     }
  283. }
  284.  
  285. /* 4.2.7 */
  286.  
  287. static void Quantization_and_coding P1((LAR),
  288.     register word * LAR        /* [0..7]    IN/OUT    */
  289. )
  290. {
  291.     register word    temp;
  292.     longword    ltmp;
  293.  
  294.  
  295.     /*  This procedure needs four tables; the following equations
  296.      *  give the optimum scaling for the constants:
  297.      *  
  298.      *  A[0..7] = integer( real_A[0..7] * 1024 )
  299.      *  B[0..7] = integer( real_B[0..7] *  512 )
  300.      *  MAC[0..7] = maximum of the LARc[0..7]
  301.      *  MIC[0..7] = minimum of the LARc[0..7]
  302.      */
  303.  
  304. #    undef STEP
  305. #    define    STEP( A, B, MAC, MIC )        \
  306.         temp = GSM_MULT( A,   *LAR );    \
  307.         temp = GSM_ADD(  temp,   B );    \
  308.         temp = GSM_ADD(  temp, 256 );    \
  309.         temp = SASR(     temp,   9 );    \
  310.         *LAR  =  temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
  311.         LAR++;
  312.  
  313.     STEP(  20480,     0,  31, -32 );
  314.     STEP(  20480,     0,  31, -32 );
  315.     STEP(  20480,  2048,  15, -16 );
  316.     STEP(  20480, -2560,  15, -16 );
  317.  
  318.     STEP(  13964,    94,   7,  -8 );
  319.     STEP(  15360, -1792,   7,  -8 );
  320.     STEP(   8534,  -341,   3,  -4 );
  321.     STEP(   9036, -1144,   3,  -4 );
  322.  
  323. #    undef    STEP
  324. }
  325.  
  326. void Gsm_LPC_Analysis P3((S, s,LARc),
  327.     struct gsm_state *S,
  328.     word          * s,        /* 0..159 signals    IN/OUT    */
  329.         word          * LARc)    /* 0..7   LARc's    OUT    */
  330. {
  331.     longword    L_ACF[9];
  332.  
  333. #if defined(USE_FLOAT_MUL) && defined(FAST)
  334.     if (S->fast) Fast_Autocorrelation (s,      L_ACF );
  335.     else
  336. #endif
  337.     Autocorrelation              (s,      L_ACF    );
  338.     Reflection_coefficients          (L_ACF, LARc    );
  339.     Transformation_to_Log_Area_Ratios (LARc);
  340.     Quantization_and_coding          (LARc);
  341. }
  342.