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
  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 static inline 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 static inline 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 static inline 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 static inline 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 static inline 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     qi*=labs(ilsp[0]-wi);
 201     pi*=labs(ilsp[1]-wi);
 202
 203     for(j=3;j<m;j+=2){
 204       if(!(shift=MLOOP_1[(pi|qi)>>25]))
 205         if(!(shift=MLOOP_2[(pi|qi)>>19]))
 206           shift=MLOOP_3[(pi|qi)>>16];
 207       qi=(qi>>shift)*labs(ilsp[j-1]-wi);
 208       pi=(pi>>shift)*labs(ilsp[j]-wi);
 209       qexp+=shift;
 210     }
 211     if(!(shift=MLOOP_1[(pi|qi)>>25]))
 212       if(!(shift=MLOOP_2[(pi|qi)>>19]))
 213         shift=MLOOP_3[(pi|qi)>>16];
 214
 215     /* pi,qi normalized collectively, both tracked using qexp */
 216
 217     if(m&1){
 218       /* odd order filter; slightly assymetric */
 219       /* the last coefficient */
 220       qi=(qi>>shift)*labs(ilsp[j-1]-wi);
 221       pi=(pi>>shift)<<14;
 222       qexp+=shift;
 223
 224       if(!(shift=MLOOP_1[(pi|qi)>>25]))
 225         if(!(shift=MLOOP_2[(pi|qi)>>19]))
 226           shift=MLOOP_3[(pi|qi)>>16];
 227
 228       pi>>=shift;
 229       qi>>=shift;
 230       qexp+=shift-14*((m+1)>>1);
 231
 232       pi=((pi*pi)>>16);
 233       qi=((qi*qi)>>16);
 234       qexp=qexp*2+m;
 235
 236       pi*=(1<<14)-((wi*wi)>>14);
 237       qi+=pi>>14;
 238
 239     }else{
 240       /* even order filter; still symmetric */
 241
 242       /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
 243          worth tracking step by step */
 244
 245       pi>>=shift;
 246       qi>>=shift;
 247       qexp+=shift-7*m;
 248
 249       pi=((pi*pi)>>16);
 250       qi=((qi*qi)>>16);
 251       qexp=qexp*2+m;
 252
 253       pi*=(1<<14)-wi;
 254       qi*=(1<<14)+wi;
 255       qi=(qi+pi)>>14;
 256
 257     }
 258
 259
 260     /* we've let the normalization drift because it wasn't important;
 261        however, for the lookup, things must be normalized again.  We
 262        need at most one right shift or a number of left shifts */
 263
 264     if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
 265       qi>>=1; qexp++;
 266     }else
 267       while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
 268         qi<<=1; qexp--;
 269       }
 270
 271 #endif
 272
 273     amp=vorbis_fromdBlook_i(ampi*                     /*  n.4         */
 274                             vorbis_invsqlook_i(qi,qexp)-
 275                                                       /*  m.8, m+n<=8 */
 276                             ampoffseti);              /*  8.12[0]     */
 277
 278 #ifdef _LOW_ACCURACY_
 279     amp>>=9;
 280 #endif
 281     curve[i]= MULT31_SHIFT15(curve[i],amp);
 282     while(map[++i]==k) curve[i]= MULT31_SHIFT15(curve[i],amp);
 283   }
 284 }
 285
 286 /*************** vorbis decode glue ************/
 287
 288 static void floor0_free_info(vorbis_info_floor *i){
 289   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 290   if(info){
 291     memset(info,0,sizeof(*info));
 292     _ogg_free(info);
 293   }
 294 }
 295
 296 static void floor0_free_look(vorbis_look_floor *i){
 297   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 298   if(look){
 299
 300     if(look->linearmap)_ogg_free(look->linearmap);
 301     if(look->lsp_look)_ogg_free(look->lsp_look);
 302     memset(look,0,sizeof(*look));
 303     _ogg_free(look);
 304   }
 305 }
 306
 307 static vorbis_info_floor *floor0_unpack (vorbis_info *vi,oggpack_buffer *opb){
 308   codec_setup_info     *ci=(codec_setup_info *)vi->codec_setup;
 309   int j;
 310
 311   vorbis_info_floor0 *info=(vorbis_info_floor0 *)_ogg_malloc(sizeof(*info));
 312   info->order=oggpack_read(opb,8);
 313   info->rate=oggpack_read(opb,16);
 314   info->barkmap=oggpack_read(opb,16);
 315   info->ampbits=oggpack_read(opb,6);
 316   info->ampdB=oggpack_read(opb,8);
 317   info->numbooks=oggpack_read(opb,4)+1;
 318
 319   if(info->order<1)goto err_out;
 320   if(info->rate<1)goto err_out;
 321   if(info->barkmap<1)goto err_out;
 322   if(info->numbooks<1)goto err_out;
 323
 324   for(j=0;j<info->numbooks;j++){
 325     info->books[j]=oggpack_read(opb,8);
 326     if(info->books[j]<0 || info->books[j]>=ci->books)goto err_out;
 327   }
 328   return(info);
 329
 330  err_out:
 331   floor0_free_info(info);
 332   return(NULL);
 333 }
 334
 335 /* initialize Bark scale and normalization lookups.  We could do this
 336    with static tables, but Vorbis allows a number of possible
 337    combinations, so it's best to do it computationally.
 338
 339    The below is authoritative in terms of defining scale mapping.
 340    Note that the scale depends on the sampling rate as well as the
 341    linear block and mapping sizes */
 342
 343 static vorbis_look_floor *floor0_look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 344                               vorbis_info_floor *i){
 345   int j;
 346   ogg_int32_t scale;
 347   vorbis_info        *vi=vd->vi;
 348   codec_setup_info   *ci=(codec_setup_info *)vi->codec_setup;
 349   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 350   vorbis_look_floor0 *look=(vorbis_look_floor0 *)_ogg_calloc(1,sizeof(*look));
 351   look->m=info->order;
 352   look->n=ci->blocksizes[mi->blockflag]/2;
 353   look->ln=info->barkmap;
 354   look->vi=info;
 355
 356   /* the mapping from a linear scale to a smaller bark scale is
 357      straightforward.  We do *not* make sure that the linear mapping
 358      does not skip bark-scale bins; the decoder simply skips them and
 359      the encoder may do what it wishes in filling them.  They're
 360      necessary in some mapping combinations to keep the scale spacing
 361      accurate */
 362   look->linearmap=(int *)_ogg_malloc((look->n+1)*sizeof(*look->linearmap));
 363   for(j=0;j<look->n;j++){
 364
 365     int val=(look->ln*
 366              ((toBARK(info->rate/2*j/look->n)<<11)/toBARK(info->rate/2)))>>11;
 367
 368     if(val>=look->ln)val=look->ln-1; /* guard against the approximation */
 369     look->linearmap[j]=val;
 370   }
 371   look->linearmap[j]=-1;
 372
 373   look->lsp_look=(ogg_int32_t *)_ogg_malloc(look->ln*sizeof(*look->lsp_look));
 374   for(j=0;j<look->ln;j++)
 375     look->lsp_look[j]=vorbis_coslook2_i(0x10000*j/look->ln);
 376
 377   return look;
 378 }
 379
 380 static void *floor0_inverse1(vorbis_block *vb,vorbis_look_floor *i){
 381   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 382   vorbis_info_floor0 *info=look->vi;
 383   int j,k;
 384
 385   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 386   if(ampraw>0){ /* also handles the -1 out of data case */
 387     long maxval=(1<<info->ampbits)-1;
 388     int amp=((ampraw*info->ampdB)<<4)/maxval;
 389     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 390
 391     if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
 392       codec_setup_info  *ci=(codec_setup_info *)vb->vd->vi->codec_setup;
 393       codebook *b=ci->fullbooks+info->books[booknum];
 394       ogg_int32_t last=0;
 395       ogg_int32_t *lsp=(ogg_int32_t *)_vorbis_block_alloc(vb,sizeof(*lsp)*(look->m+1));
 396
 397       for(j=0;j<look->m;j+=b->dim)
 398         if(vorbis_book_decodev_set(b,lsp+j,&vb->opb,b->dim,-24)==-1)goto eop;
 399       for(j=0;j<look->m;){
 400         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 401         last=lsp[j-1];
 402       }
 403
 404       lsp[look->m]=amp;
 405       return(lsp);
 406     }
 407   }
 408  eop:
 409   return(NULL);
 410 }
 411
 412 static int floor0_inverse2(vorbis_block *vb,vorbis_look_floor *i,
 413                            void *memo,ogg_int32_t *out){
 414   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 415   vorbis_info_floor0 *info=look->vi;
 416
 417   if(memo){
 418     ogg_int32_t *lsp=(ogg_int32_t *)memo;
 419     ogg_int32_t amp=lsp[look->m];
 420
 421     /* take the coefficients back to a spectral envelope curve */
 422     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 423                         lsp,look->m,amp,info->ampdB,look->lsp_look);
 424     return(1);
 425   }
 426   memset(out,0,sizeof(*out)*look->n);
 427   return(0);
 428 }
 429
 430 /* export hooks */
 431 vorbis_func_floor floor0_exportbundle={
 432   &floor0_unpack,&floor0_look,&floor0_free_info,
 433   &floor0_free_look,&floor0_inverse1,&floor0_inverse2
 434 };
 435
 436