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