MythTV: short_term.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 #ifdef K6OPT
00017 #include "k6opt.h"
00018 
00019 #define Short_term_analysis_filtering Short_term_analysis_filteringx
00020 
00021 #endif
00022 /*
00023  *  SHORT TERM ANALYSIS FILTERING SECTION
00024  */
00025 
00026 /* 4.2.8 */
00027 
00028 static void Decoding_of_the_coded_Log_Area_Ratios P2((LARc,LARpp),
00029         word    * LARc,         /* coded log area ratio [0..7]  IN      */
00030         word    * LARpp)        /* out: decoded ..                      */
00031 {
00032         register word   temp1 /* , temp2 */;
00033         register long   ltmp;   /* for GSM_ADD */
00034 
00035         /*  This procedure requires for efficient implementation
00036          *  two tables.
00037          *
00038          *  INVA[1..8] = integer( (32768 * 8) / real_A[1..8])
00039          *  MIC[1..8]  = minimum value of the LARc[1..8]
00040          */
00041 
00042         /*  Compute the LARpp[1..8]
00043          */
00044 
00045         /*      for (i = 1; i <= 8; i++, B++, MIC++, INVA++, LARc++, LARpp++) {
00046          *
00047          *              temp1  = GSM_ADD( *LARc, *MIC ) << 10;
00048          *              temp2  = *B << 1;
00049          *              temp1  = GSM_SUB( temp1, temp2 );
00050          *
00051          *              assert(*INVA != MIN_WORD);
00052          *
00053          *              temp1  = GSM_MULT_R( *INVA, temp1 );
00054          *              *LARpp = GSM_ADD( temp1, temp1 );
00055          *      }
00056          */
00057 
00058 #undef  STEP
00059 #define STEP( B, MIC, INVA )    \
00060                 temp1    = GSM_ADD( *LARc++, MIC ) << 10;       \
00061                 temp1    = GSM_SUB( temp1, B << 1 );            \
00062                 temp1    = GSM_MULT_R( INVA, temp1 );           \
00063                 *LARpp++ = GSM_ADD( temp1, temp1 );
00064 
00065         STEP(      0,  -32,  13107 );
00066         STEP(      0,  -32,  13107 );
00067         STEP(   2048,  -16,  13107 );
00068         STEP(  -2560,  -16,  13107 );
00069 
00070         STEP(     94,   -8,  19223 );
00071         STEP(  -1792,   -8,  17476 );
00072         STEP(   -341,   -4,  31454 );
00073         STEP(  -1144,   -4,  29708 );
00074 
00075         /* NOTE: the addition of *MIC is used to restore
00076          *       the sign of *LARc.
00077          */
00078 }
00079 
00080 /* 4.2.9 */
00081 /* Computation of the quantized reflection coefficients 
00082  */
00083 
00084 /* 4.2.9.1  Interpolation of the LARpp[1..8] to get the LARp[1..8]
00085  */
00086 
00087 /*
00088  *  Within each frame of 160 analyzed speech samples the short term
00089  *  analysis and synthesis filters operate with four different sets of
00090  *  coefficients, derived from the previous set of decoded LARs(LARpp(j-1))
00091  *  and the actual set of decoded LARs (LARpp(j))
00092  *
00093  * (Initial value: LARpp(j-1)[1..8] = 0.)
00094  */
00095 
00096 static void Coefficients_0_12 P3((LARpp_j_1, LARpp_j, LARp),
00097         register word * LARpp_j_1,
00098         register word * LARpp_j,
00099         register word * LARp)
00100 {
00101         register int    i;
00102         register longword ltmp;
00103 
00104         for (i = 1; i <= 8; i++, LARp++, LARpp_j_1++, LARpp_j++) {
00105                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
00106                 *LARp = GSM_ADD( *LARp,  SASR( *LARpp_j_1, 1));
00107         }
00108 }
00109 
00110 static void Coefficients_13_26 P3((LARpp_j_1, LARpp_j, LARp),
00111         register word * LARpp_j_1,
00112         register word * LARpp_j,
00113         register word * LARp)
00114 {
00115         register int i;
00116         register longword ltmp;
00117         for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
00118                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 1), SASR( *LARpp_j, 1 ));
00119         }
00120 }
00121 
00122 static void Coefficients_27_39 P3((LARpp_j_1, LARpp_j, LARp),
00123         register word * LARpp_j_1,
00124         register word * LARpp_j,
00125         register word * LARp)
00126 {
00127         register int i;
00128         register longword ltmp;
00129 
00130         for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
00131                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
00132                 *LARp = GSM_ADD( *LARp, SASR( *LARpp_j, 1 ));
00133         }
00134 }
00135 
00136 
00137 static void Coefficients_40_159 P2((LARpp_j, LARp),
00138         register word * LARpp_j,
00139         register word * LARp)
00140 {
00141         register int i;
00142 
00143         for (i = 1; i <= 8; i++, LARp++, LARpp_j++)
00144                 *LARp = *LARpp_j;
00145 }
00146 
00147 /* 4.2.9.2 */
00148 
00149 static void LARp_to_rp P1((LARp),
00150         register word * LARp)   /* [0..7] IN/OUT  */
00151 /*
00152  *  The input of this procedure is the interpolated LARp[0..7] array.
00153  *  The reflection coefficients, rp[i], are used in the analysis
00154  *  filter and in the synthesis filter.
00155  */
00156 {
00157         register int            i;
00158         register word           temp;
00159         register longword       ltmp;
00160 
00161         for (i = 1; i <= 8; i++, LARp++) {
00162 
00163                 /* temp = GSM_ABS( *LARp );
00164                  *
00165                  * if (temp < 11059) temp <<= 1;
00166                  * else if (temp < 20070) temp += 11059;
00167                  * else temp = GSM_ADD( temp >> 2, 26112 );
00168                  *
00169                  * *LARp = *LARp < 0 ? -temp : temp;
00170                  */
00171 
00172                 if (*LARp < 0) {
00173                         temp = *LARp == MIN_WORD ? MAX_WORD : -(*LARp);
00174                         *LARp = - ((temp < 11059) ? temp << 1
00175                                 : ((temp < 20070) ? temp + 11059
00176                                 :  GSM_ADD( temp >> 2, 26112 )));
00177                 } else {
00178                         temp  = *LARp;
00179                         *LARp =    (temp < 11059) ? temp << 1
00180                                 : ((temp < 20070) ? temp + 11059
00181                                 :  GSM_ADD( temp >> 2, 26112 ));
00182                 }
00183         }
00184 }
00185 
00186 
00187 /* 4.2.10 */
00188 #ifndef Short_term_analysis_filtering
00189 
00190 /* SJB Remark:
00191  * I tried 2 MMX versions of this function, neither is significantly
00192  * faster than the C version which follows.  MMX might be useful if
00193  * one were processing 2 input streams in parallel.
00194  */
00195 static void Short_term_analysis_filtering P4((u0,rp0,k_n,s),
00196         register word * u0,
00197         register word   * rp0,  /* [0..7]       IN      */
00198         register int    k_n,    /*   k_end - k_start    */
00199         register word   * s     /* [0..n-1]     IN/OUT  */
00200 )
00201 /*
00202  *  This procedure computes the short term residual signal d[..] to be fed
00203  *  to the RPE-LTP loop from the s[..] signal and from the local rp[..]
00204  *  array (quantized reflection coefficients).  As the call of this
00205  *  procedure can be done in many ways (see the interpolation of the LAR
00206  *  coefficient), it is assumed that the computation begins with index
00207  *  k_start (for arrays d[..] and s[..]) and stops with index k_end
00208  *  (k_start and k_end are defined in 4.2.9.1).  This procedure also
00209  *  needs to keep the array u0[0..7] in memory for each call.
00210  */
00211 {
00212         register word           * u_top = u0 + 8;
00213         register word           * s_top = s + k_n;
00214 
00215         while (s < s_top) {
00216                 register word           *u, *rp ;
00217                 register longword               di, u_out;
00218                 di = u_out = *s;
00219                 for (rp=rp0, u=u0; u<u_top;) {
00220                         register longword       ui, rpi;
00221                         ui    = *u;
00222                         *u++  = u_out;
00223                         rpi   = *rp++;
00224                         u_out = ui + (((rpi*di)+0x4000)>>15);
00225                         di    = di + (((rpi*ui)+0x4000)>>15);
00226                         /* make the common case fastest: */
00227                         if ((u_out == (word)u_out) && (di == (word)di)) continue;
00228                         /* otherwise do slower fixup (saturation) */
00229                         if (u_out>MAX_WORD) u_out=MAX_WORD;
00230                         else if (u_out<MIN_WORD) u_out=MIN_WORD;
00231                         if (di>MAX_WORD) di=MAX_WORD;
00232                         else if (di<MIN_WORD) di=MIN_WORD;
00233                 }
00234                 *s++ = di;
00235         }
00236 }
00237 #endif
00238 
00239 #if defined(USE_FLOAT_MUL) && defined(FAST)
00240 
00241 static void Fast_Short_term_analysis_filtering P4((u,rp,k_n,s),
00242         register word * u;
00243         register word   * rp,   /* [0..7]       IN      */
00244         register int    k_n,    /*   k_end - k_start    */
00245         register word   * s     /* [0..n-1]     IN/OUT  */
00246 )
00247 {
00248         register int            i;
00249 
00250         float     uf[8],
00251                  rpf[8];
00252 
00253         register float scalef = 3.0517578125e-5;
00254         register float          sav, di, temp;
00255 
00256         for (i = 0; i < 8; ++i) {
00257                 uf[i]  = u[i];
00258                 rpf[i] = rp[i] * scalef;
00259         }
00260         for (; k_n--; s++) {
00261                 sav = di = *s;
00262                 for (i = 0; i < 8; ++i) {
00263                         register float rpfi = rpf[i];
00264                         register float ufi  = uf[i];
00265 
00266                         uf[i] = sav;
00267                         temp  = rpfi * di + ufi;
00268                         di   += rpfi * ufi;
00269                         sav   = temp;
00270                 }
00271                 *s = di;
00272         }
00273         for (i = 0; i < 8; ++i) u[i] = uf[i];
00274 }
00275 #endif /* ! (defined (USE_FLOAT_MUL) && defined (FAST)) */
00276 
00277 /*
00278  * SJB Remark: modified Short_term_synthesis_filtering() below
00279  *  for significant (abt 35%) speedup of decompression.
00280  *    (gcc-2.95, k6 cpu)
00281  *  Please don't change this without benchmarking decompression
00282  *  to see that you haven't harmed speed.
00283  *  This function burns most of CPU time for untoasting.
00284  *  Unfortunately, didn't see any good way to benefit from mmx.
00285  */
00286 static void Short_term_synthesis_filtering P5((S,rrp,k,wt,sr),
00287         struct gsm_state * S,
00288         register word   * rrp,  /* [0..7]       IN      */
00289         register int    k,      /* k_end - k_start      */
00290         register word   * wt,   /* [0..k-1]     IN      */
00291         register word   * sr    /* [0..k-1]     OUT     */
00292 )
00293 {
00294         register word           * v = S->v;
00295         register int            i;
00296         register longword               sri;
00297 
00298         while (k--) {
00299                 sri = *wt++;
00300                 for (i = 8; i--;) {
00301                         register longword               tmp1, tmp2;
00302 
00303                         /* sri = GSM_SUB( sri, gsm_mult_r( rrp[i], v[i] ) );
00304                          */
00305                         tmp1 = rrp[i];
00306                         tmp2 = v[i];
00307 
00308                         tmp2 = (( tmp1 * tmp2 + 16384) >> 15) ;
00309                         /* saturation done below */
00310                         sri  -= tmp2;
00311                         if (sri != (word)sri) {
00312                                 sri = (sri<0)? MIN_WORD:MAX_WORD;
00313                         }
00314                         /* v[i+1] = GSM_ADD( v[i], gsm_mult_r( rrp[i], sri ) );
00315                          */
00316 
00317                         tmp1 = (( tmp1 * sri + 16384) >> 15) ;
00318                         /* saturation done below */
00319                         tmp1 += v[i];
00320                         if (tmp1 != (word)tmp1) {
00321                                 tmp1 = (tmp1<0)? MIN_WORD:MAX_WORD;
00322                         }
00323                         v[i+1] = tmp1;
00324                 }
00325                 *sr++ = v[0] = sri;
00326         }
00327 }
00328 
00329 
00330 #if defined(FAST) && defined(USE_FLOAT_MUL)
00331 
00332 static void Fast_Short_term_synthesis_filtering P5((S,rrp,k,wt,sr),
00333         struct gsm_state * S,
00334         register word   * rrp,  /* [0..7]       IN      */
00335         register int    k,      /* k_end - k_start      */
00336         register word   * wt,   /* [0..k-1]     IN      */
00337         register word   * sr    /* [0..k-1]     OUT     */
00338 )
00339 {
00340         register word           * v = S->v;
00341         register int            i;
00342 
00343         float va[9], rrpa[8];
00344         register float scalef = 3.0517578125e-5, temp;
00345 
00346         for (i = 0; i < 8; ++i) {
00347                 va[i]   = v[i];
00348                 rrpa[i] = (float)rrp[i] * scalef;
00349         }
00350         while (k--) {
00351                 register float sri = *wt++;
00352                 for (i = 8; i--;) {
00353                         sri -= rrpa[i] * va[i];
00354                         if     (sri < -32768.) sri = -32768.;
00355                         else if (sri > 32767.) sri =  32767.;
00356 
00357                         temp = va[i] + rrpa[i] * sri;
00358                         if     (temp < -32768.) temp = -32768.;
00359                         else if (temp > 32767.) temp =  32767.;
00360                         va[i+1] = temp;
00361                 }
00362                 *sr++ = va[0] = sri;
00363         }
00364         for (i = 0; i < 9; ++i) v[i] = va[i];
00365 }
00366 
00367 #endif /* defined(FAST) && defined(USE_FLOAT_MUL) */
00368 
00369 void Gsm_Short_Term_Analysis_Filter P3((S,LARc,s),
00370 
00371         struct gsm_state * S,
00372 
00373         word    * LARc,         /* coded log area ratio [0..7]  IN      */
00374         word    * s             /* signal [0..159]              IN/OUT  */
00375 )
00376 {
00377         word            * LARpp_j       = S->LARpp[ S->j      ];
00378         word            * LARpp_j_1     = S->LARpp[ S->j ^= 1 ];
00379 
00380         word            LARp[8];
00381 int i;
00382 #undef  FILTER
00383 #if     defined(FAST) && defined(USE_FLOAT_MUL)
00384 #       define  FILTER  (* (S->fast                     \
00385                            ? Fast_Short_term_analysis_filtering \
00386                            : Short_term_analysis_filtering      ))
00387 
00388 #else
00389 #       define  FILTER  Short_term_analysis_filtering
00390 #endif
00391 
00392         Decoding_of_the_coded_Log_Area_Ratios( LARc, LARpp_j );
00393 
00394         Coefficients_0_12(  LARpp_j_1, LARpp_j, LARp );
00395         LARp_to_rp( LARp );
00396         FILTER( S->u, LARp, 13, s);
00397 
00398         Coefficients_13_26( LARpp_j_1, LARpp_j, LARp);
00399         LARp_to_rp( LARp );
00400         FILTER( S->u, LARp, 14, s + 13);
00401 
00402         Coefficients_27_39( LARpp_j_1, LARpp_j, LARp);
00403         LARp_to_rp( LARp );
00404         FILTER( S->u, LARp, 13, s + 27);
00405 
00406         Coefficients_40_159( LARpp_j, LARp);
00407         LARp_to_rp( LARp );
00408         FILTER( S->u, LARp, 120, s + 40);
00409         
00410 }
00411 
00412 void Gsm_Short_Term_Synthesis_Filter P4((S, LARcr, wt, s),
00413         struct gsm_state * S,
00414 
00415         word    * LARcr,        /* received log area ratios [0..7] IN  */
00416         word    * wt,           /* received d [0..159]             IN  */
00417 
00418         word    * s             /* signal   s [0..159]            OUT  */
00419 )
00420 {
00421         word            * LARpp_j       = S->LARpp[ S->j     ];
00422         word            * LARpp_j_1     = S->LARpp[ S->j ^=1 ];
00423 
00424         word            LARp[8];
00425 
00426 #undef  FILTER
00427 #if     defined(FAST) && defined(USE_FLOAT_MUL)
00428 
00429 #       define  FILTER  (* (S->fast                     \
00430                            ? Fast_Short_term_synthesis_filtering        \
00431                            : Short_term_synthesis_filtering     ))
00432 #else
00433 #       define  FILTER  Short_term_synthesis_filtering
00434 #endif
00435 
00436         Decoding_of_the_coded_Log_Area_Ratios( LARcr, LARpp_j );
00437 
00438         Coefficients_0_12( LARpp_j_1, LARpp_j, LARp );
00439         LARp_to_rp( LARp );
00440         FILTER( S, LARp, 13, wt, s );
00441 
00442         Coefficients_13_26( LARpp_j_1, LARpp_j, LARp);
00443         LARp_to_rp( LARp );
00444         FILTER( S, LARp, 14, wt + 13, s + 13 );
00445 
00446         Coefficients_27_39( LARpp_j_1, LARpp_j, LARp);
00447         LARp_to_rp( LARp );
00448         FILTER( S, LARp, 13, wt + 27, s + 27 );
00449 
00450         Coefficients_40_159( LARpp_j, LARp );
00451         LARp_to_rp( LARp );
00452         FILTER(S, LARp, 120, wt + 40, s + 40);
00453 }