Tremor/floor0.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis 'TREMOR' CODEC SOURCE CODE.   *
   4  *                                                                  *
   5  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   6  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   7  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   8  *                                                                  *
   9  * THE OggVorbis 'TREMOR' SOURCE CODE IS (C) COPYRIGHT 1994-2002    *
  10  * BY THE Xiph.Org FOUNDATION http://www.xiph.org/                  *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: floor backend 0 implementation
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <string.h>
  20 #include <math.h>
  21 #include <ogg/ogg.h>
  22 #include "ivorbiscodec.h"
  23 #include "codec_internal.h"
  24 #include "registry.h"
  25 #include "codebook.h"
  26 #include "misc.h"
  27 #include "block.h"
  28
  29 #define LSP_FRACBITS 14
  30
  31 typedef struct {
  32   long n;
  33   int ln;
  34   int  m;
  35   int *linearmap;
  36
  37   vorbis_info_floor0 *vi;
  38   ogg_int32_t *lsp_look;
  39
  40 } vorbis_look_floor0;
  41
  42 /*************** LSP decode ********************/
  43
  44 #include "lsp_lookup.h"
  45
  46 /* interpolated 1./sqrt(p) where .5 <= a < 1. (.100000... to .111111...) in
  47    16.16 format
  48    returns in m.8 format */
  49
  50 static long ADJUST_SQRT2[2]={8192,5792};
  51 STIN ogg_int32_t vorbis_invsqlook_i(long a,long e){
  52   long i=(a&0x7fff)>>(INVSQ_LOOKUP_I_SHIFT-1);
  53   long d=a&INVSQ_LOOKUP_I_MASK;                              /*  0.10 */
  54   long val=INVSQ_LOOKUP_I[i]-                                /*  1.16 */
  55     ((INVSQ_LOOKUP_IDel[i]*d)>>INVSQ_LOOKUP_I_SHIFT);        /* result 1.16 */
  56   val*=ADJUST_SQRT2[e&1];
  57   e=(e>>1)+21;
  58   return(val>>e);
  59 }
  60
  61 /* interpolated lookup based fromdB function, domain -140dB to 0dB only */
  62 /* a is in n.12 format */
  63 STIN ogg_int32_t vorbis_fromdBlook_i(long a){
  64   int i=(-a)>>(12-FROMdB2_SHIFT);
  65   if(i<0) return 0x7fffffff;
  66   if(i>=(FROMdB_LOOKUP_SZ<<FROMdB_SHIFT))return 0;
  67
  68   return FROMdB_LOOKUP[i>>FROMdB_SHIFT] * FROMdB2_LOOKUP[i&FROMdB2_MASK];
  69 }
  70
  71 /* interpolated lookup based cos function, domain 0 to PI only */
  72 /* a is in 0.16 format, where 0==0, 2^^16-1==PI, return 0.14 */
  73 STIN ogg_int32_t vorbis_coslook_i(long a){
  74   int i=a>>COS_LOOKUP_I_SHIFT;
  75   int d=a&COS_LOOKUP_I_MASK;
  76   return COS_LOOKUP_I[i]- ((d*(COS_LOOKUP_I[i]-COS_LOOKUP_I[i+1]))>>
  77                            COS_LOOKUP_I_SHIFT);
  78 }
  79
  80 /* interpolated lookup based cos function */
  81 /* a is in 0.16 format, where 0==0, 2^^16==PI, return .LSP_FRACBITS */
  82 STIN ogg_int32_t vorbis_coslook2_i(long a){
  83   a=a&0x1ffff;
  84
  85   if(a>0x10000)a=0x20000-a;
  86   {
  87     int i=a>>COS_LOOKUP_I_SHIFT;
  88     int d=a&COS_LOOKUP_I_MASK;
  89     a=((COS_LOOKUP_I[i]<<COS_LOOKUP_I_SHIFT)-
  90        d*(COS_LOOKUP_I[i]-COS_LOOKUP_I[i+1]))>>
  91       (COS_LOOKUP_I_SHIFT-LSP_FRACBITS+14);
  92   }
  93
  94   return(a);
  95 }
  96
  97 static const int barklook[28]={
  98   0,100,200,301,          405,516,635,766,
  99   912,1077,1263,1476,     1720,2003,2333,2721,
 100   3184,3742,4428,5285,    6376,7791,9662,12181,
 101   15624,20397,27087,36554
 102 };
 103
 104 /* used in init only; interpolate the long way */
 105 STIN ogg_int32_t toBARK(int n){
 106   int i;
 107   for(i=0;i<27;i++)
 108     if(n>=barklook[i] && n<barklook[i+1])break;
 109
 110   if(i==27){
 111     return 27<<15;
 112   }else{
 113     int gap=barklook[i+1]-barklook[i];
 114     int del=n-barklook[i];
 115
 116     return((i<<15)+((del<<15)/gap));
 117   }
 118 }
 119
 120 static const unsigned char MLOOP_1[64]={
 121    0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
 122   14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
 123   15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
 124   15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
 125 };
 126
 127 static const unsigned char MLOOP_2[64]={
 128   0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
 129   8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
 130   9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
 131   9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
 132 };
 133
 134 static const unsigned char MLOOP_3[8]={0,1,2,2,3,3,3,3};
 135
 136 void vorbis_lsp_to_curve(ogg_int32_t *curve,int *map,int n,int ln,
 137                          ogg_int32_t *lsp,int m,
 138                          ogg_int32_t amp,
 139                          ogg_int32_t ampoffset,
 140                          ogg_int32_t *icos){
 141
 142   /* 0 <= m < 256 */
 143
 144   /* set up for using all int later */
 145   int i;
 146   int ampoffseti=ampoffset*4096;
 147   int ampi=amp;
 148   ogg_int32_t *ilsp=(ogg_int32_t *)alloca(m*sizeof(*ilsp));
 149   /* lsp is in 8.24, range 0 to PI; coslook wants it in .16 0 to 1*/
 150   for(i=0;i<m;i++){
 151 #ifndef _LOW_ACCURACY_
 152     ogg_int32_t val=MULT32(lsp[i],0x517cc2);
 153 #else
 154     ogg_int32_t val=((lsp[i]>>10)*0x517d)>>14;
 155 #endif
 156
 157     /* safeguard against a malicious stream */
 158     if(val<0 || (val>>COS_LOOKUP_I_SHIFT)>=COS_LOOKUP_I_SZ){
 159       memset(curve,0,sizeof(*curve)*n);
 160       return;
 161     }
 162
 163     ilsp[i]=vorbis_coslook_i(val);
 164   }
 165
 166   i=0;
 167   while(i<n){
 168     int j,k=map[i];
 169     ogg_uint32_t pi=46341; /* 2**-.5 in 0.16 */
 170     ogg_uint32_t qi=46341;
 171     ogg_int32_t qexp=0,shift;
 172     ogg_int32_t wi=icos[k];
 173
 174 #ifdef _V_LSP_MATH_ASM
 175     lsp_loop_asm(&qi,&pi,&qexp,ilsp,wi,m);
 176
 177     pi=((pi*pi)>>16);
 178     qi=((qi*qi)>>16);
 179
 180     if(m&1){
 181       qexp= qexp*2-28*((m+1)>>1)+m;
 182       pi*=(1<<14)-((wi*wi)>>14);
 183       qi+=pi>>14;
 184     }else{
 185       qexp= qexp*2-13*m;
 186
 187       pi*=(1<<14)-wi;
 188       qi*=(1<<14)+wi;
 189
 190       qi=(qi+pi)>>14;
 191     }
 192
 193     if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
 194       qi>>=1; qexp++;
 195     }else
 196       lsp_norm_asm(&qi,&qexp);
 197
 198 #else
 199
 200     j=1;
 201     if(m>1){
 202       qi*=labs(ilsp[0]-wi);
 203       pi*=labs(ilsp[1]-wi);
 204
 205       for(j+=2;j<m;j+=2){
 206         if(!(shift=MLOOP_1[(pi|qi)>>25]))
 207           if(!(shift=MLOOP_2[(pi|qi)>>19]))
 208             shift=MLOOP_3[(pi|qi)>>16];
 209         qi=(qi>>shift)*labs(ilsp[j-1]-wi);
 210         pi=(pi>>shift)*labs(ilsp[j]-wi);
 211         qexp+=shift;
 212       }
 213     }
 214     if(!(shift=MLOOP_1[(pi|qi)>>25]))
 215       if(!(shift=MLOOP_2[(pi|qi)>>19]))
 216         shift=MLOOP_3[(pi|qi)>>16];
 217
 218     /* pi,qi normalized collectively, both tracked using qexp */
 219
 220     if(m&1){
 221       /* odd order filter; slightly assymetric */
 222       /* the last coefficient */
 223       qi=(qi>>shift)*labs(ilsp[j-1]-wi);
 224       pi=(pi>>shift)<<14;
 225       qexp+=shift;
 226
 227       if(!(shift=MLOOP_1[(pi|qi)>>25]))
 228         if(!(shift=MLOOP_2[(pi|qi)>>19]))
 229           shift=MLOOP_3[(pi|qi)>>16];
 230
 231       pi>>=shift;
 232       qi>>=shift;
 233       qexp+=shift-14*((m+1)>>1);
 234
 235       pi=((pi*pi)>>16);
 236       qi=((qi*qi)>>16);
 237       qexp=qexp*2+m;
 238
 239       pi*=(1<<14)-((wi*wi)>>14);
 240       qi+=pi>>14;
 241
 242     }else{
 243       /* even order filter; still symmetric */
 244
 245       /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
 246          worth tracking step by step */
 247
 248       pi>>=shift;
 249       qi>>=shift;
 250       qexp+=shift-7*m;
 251
 252       pi=((pi*pi)>>16);
 253       qi=((qi*qi)>>16);
 254       qexp=qexp*2+m;
 255
 256       pi*=(1<<14)-wi;
 257       qi*=(1<<14)+wi;
 258       qi=(qi+pi)>>14;
 259
 260     }
 261
 262
 263     /* we've let the normalization drift because it wasn't important;
 264        however, for the lookup, things must be normalized again.  We
 265        need at most one right shift or a number of left shifts */
 266
 267     if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
 268       qi>>=1; qexp++;
 269     }else
 270       while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
 271         qi<<=1; qexp--;
 272       }
 273
 274 #endif
 275
 276     amp=vorbis_fromdBlook_i(ampi*                     /*  n.4         */
 277                             vorbis_invsqlook_i(qi,qexp)-
 278                                                       /*  m.8, m+n<=8 */
 279                             ampoffseti);              /*  8.12[0]     */
 280
 281 #ifdef _LOW_ACCURACY_
 282     amp>>=9;
 283 #endif
 284     curve[i]= MULT31_SHIFT15(curve[i],amp);
 285     while(map[++i]==k) curve[i]= MULT31_SHIFT15(curve[i],amp);
 286   }
 287 }
 288
 289 /*************** vorbis decode glue ************/
 290
 291 static void floor0_free_info(vorbis_info_floor *i){
 292   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 293   if(info){
 294     memset(info,0,sizeof(*info));
 295     _ogg_free(info);
 296   }
 297 }
 298
 299 static void floor0_free_look(vorbis_look_floor *i){
 300   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 301   if(look){
 302
 303     if(look->linearmap)_ogg_free(look->linearmap);
 304     if(look->lsp_look)_ogg_free(look->lsp_look);
 305     memset(look,0,sizeof(*look));
 306     _ogg_free(look);
 307   }
 308 }
 309
 310 static vorbis_info_floor *floor0_unpack (vorbis_info *vi,oggpack_buffer *opb){
 311   codec_setup_info     *ci=(codec_setup_info *)vi->codec_setup;
 312   int j;
 313
 314   vorbis_info_floor0 *info=(vorbis_info_floor0 *)_ogg_malloc(sizeof(*info));
 315   info->order=oggpack_read(opb,8);
 316   info->rate=oggpack_read(opb,16);
 317   info->barkmap=oggpack_read(opb,16);
 318   info->ampbits=oggpack_read(opb,6);
 319   info->ampdB=oggpack_read(opb,8);
 320   info->numbooks=oggpack_read(opb,4)+1;
 321
 322   if(info->order<1)goto err_out;
 323   if(info->rate<1)goto err_out;
 324   if(info->barkmap<1)goto err_out;
 325   if(info->numbooks<1)goto err_out;
 326
 327   for(j=0;j<info->numbooks;j++){
 328     info->books[j]=oggpack_read(opb,8);
 329     if(info->books[j]<0 || info->books[j]>=ci->books)goto err_out;
 330     if(ci->book_param[info->books[j]]->maptype==0)goto err_out;
 331     if(ci->book_param[info->books[j]]->dim<1)goto err_out;
 332   }
 333   return(info);
 334
 335  err_out:
 336   floor0_free_info(info);
 337   return(NULL);
 338 }
 339
 340 /* initialize Bark scale and normalization lookups.  We could do this
 341    with static tables, but Vorbis allows a number of possible
 342    combinations, so it's best to do it computationally.
 343
 344    The below is authoritative in terms of defining scale mapping.
 345    Note that the scale depends on the sampling rate as well as the
 346    linear block and mapping sizes */
 347
 348 static vorbis_look_floor *floor0_look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 349                               vorbis_info_floor *i){
 350   int j;
 351   vorbis_info        *vi=vd->vi;
 352   codec_setup_info   *ci=(codec_setup_info *)vi->codec_setup;
 353   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 354   vorbis_look_floor0 *look=(vorbis_look_floor0 *)_ogg_calloc(1,sizeof(*look));
 355   look->m=info->order;
 356   look->n=ci->blocksizes[mi->blockflag]/2;
 357   look->ln=info->barkmap;
 358   look->vi=info;
 359
 360   /* the mapping from a linear scale to a smaller bark scale is
 361      straightforward.  We do *not* make sure that the linear mapping
 362      does not skip bark-scale bins; the decoder simply skips them and
 363      the encoder may do what it wishes in filling them.  They're
 364      necessary in some mapping combinations to keep the scale spacing
 365      accurate */
 366   look->linearmap=(int *)_ogg_malloc((look->n+1)*sizeof(*look->linearmap));
 367   for(j=0;j<look->n;j++){
 368
 369     int val=(look->ln*
 370              ((toBARK(info->rate/2*j/look->n)<<11)/toBARK(info->rate/2)))>>11;
 371
 372     if(val>=look->ln)val=look->ln-1; /* guard against the approximation */
 373     look->linearmap[j]=val;
 374   }
 375   look->linearmap[j]=-1;
 376
 377   look->lsp_look=(ogg_int32_t *)_ogg_malloc(look->ln*sizeof(*look->lsp_look));
 378   for(j=0;j<look->ln;j++)
 379     look->lsp_look[j]=vorbis_coslook2_i(0x10000*j/look->ln);
 380
 381   return look;
 382 }
 383
 384 static void *floor0_inverse1(vorbis_block *vb,vorbis_look_floor *i){
 385   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 386   vorbis_info_floor0 *info=look->vi;
 387   int j,k;
 388
 389   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 390   if(ampraw>0){ /* also handles the -1 out of data case */
 391     long maxval=(1<<info->ampbits)-1;
 392     int amp=((ampraw*info->ampdB)<<4)/maxval;
 393     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 394
 395     if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
 396       codec_setup_info  *ci=(codec_setup_info *)vb->vd->vi->codec_setup;
 397       codebook *b=ci->fullbooks+info->books[booknum];
 398       ogg_int32_t last=0;
 399       ogg_int32_t *lsp=(ogg_int32_t *)_vorbis_block_alloc(vb,sizeof(*lsp)*(look->m+1));
 400
 401       if(vorbis_book_decodev_set(b,lsp,&vb->opb,look->m,-24)==-1)goto eop;
 402       for(j=0;j<look->m;){
 403         for(k=0;j<look->m && k<b->dim;k++,j++)lsp[j]+=last;
 404         last=lsp[j-1];
 405       }
 406
 407       lsp[look->m]=amp;
 408       return(lsp);
 409     }
 410   }
 411  eop:
 412   return(NULL);
 413 }
 414
 415 static int floor0_inverse2(vorbis_block *vb,vorbis_look_floor *i,
 416                            void *memo,ogg_int32_t *out){
 417   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 418   vorbis_info_floor0 *info=look->vi;
 419
 420   if(memo){
 421     ogg_int32_t *lsp=(ogg_int32_t *)memo;
 422     ogg_int32_t amp=lsp[look->m];
 423
 424     /* take the coefficients back to a spectral envelope curve */
 425     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 426                         lsp,look->m,amp,info->ampdB,look->lsp_look);
 427     return(1);
 428   }
 429   memset(out,0,sizeof(*out)*look->n);
 430   return(0);
 431 }
 432
 433 /* export hooks */
 434 vorbis_func_floor floor0_exportbundle={
 435   &floor0_unpack,&floor0_look,&floor0_free_info,
 436   &floor0_free_look,&floor0_inverse1,&floor0_inverse2
 437 };
 438
 439