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 "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     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   vorbis_info        *vi=vd->vi;
 347   codec_setup_info   *ci=(codec_setup_info *)vi->codec_setup;
 348   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 349   vorbis_look_floor0 *look=(vorbis_look_floor0 *)_ogg_calloc(1,sizeof(*look));
 350   look->m=info->order;
 351   look->n=ci->blocksizes[mi->blockflag]/2;
 352   look->ln=info->barkmap;
 353   look->vi=info;
 354
 355   /* the mapping from a linear scale to a smaller bark scale is
 356      straightforward.  We do *not* make sure that the linear mapping
 357      does not skip bark-scale bins; the decoder simply skips them and
 358      the encoder may do what it wishes in filling them.  They're
 359      necessary in some mapping combinations to keep the scale spacing
 360      accurate */
 361   look->linearmap=(int *)_ogg_malloc((look->n+1)*sizeof(*look->linearmap));
 362   for(j=0;j<look->n;j++){
 363
 364     int val=(look->ln*
 365              ((toBARK(info->rate/2*j/look->n)<<11)/toBARK(info->rate/2)))>>11;
 366
 367     if(val>=look->ln)val=look->ln-1; /* guard against the approximation */
 368     look->linearmap[j]=val;
 369   }
 370   look->linearmap[j]=-1;
 371
 372   look->lsp_look=(ogg_int32_t *)_ogg_malloc(look->ln*sizeof(*look->lsp_look));
 373   for(j=0;j<look->ln;j++)
 374     look->lsp_look[j]=vorbis_coslook2_i(0x10000*j/look->ln);
 375
 376   return look;
 377 }
 378
 379 static void *floor0_inverse1(vorbis_block *vb,vorbis_look_floor *i){
 380   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 381   vorbis_info_floor0 *info=look->vi;
 382   int j,k;
 383
 384   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 385   if(ampraw>0){ /* also handles the -1 out of data case */
 386     long maxval=(1<<info->ampbits)-1;
 387     int amp=((ampraw*info->ampdB)<<4)/maxval;
 388     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 389
 390     if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
 391       codec_setup_info  *ci=(codec_setup_info *)vb->vd->vi->codec_setup;
 392       codebook *b=ci->fullbooks+info->books[booknum];
 393       ogg_int32_t last=0;
 394       ogg_int32_t *lsp=(ogg_int32_t *)_vorbis_block_alloc(vb,sizeof(*lsp)*(look->m+1));
 395
 396       for(j=0;j<look->m;j+=b->dim)
 397         if(vorbis_book_decodev_set(b,lsp+j,&vb->opb,b->dim,-24)==-1)goto eop;
 398       for(j=0;j<look->m;){
 399         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 400         last=lsp[j-1];
 401       }
 402
 403       lsp[look->m]=amp;
 404       return(lsp);
 405     }
 406   }
 407  eop:
 408   return(NULL);
 409 }
 410
 411 static int floor0_inverse2(vorbis_block *vb,vorbis_look_floor *i,
 412                            void *memo,ogg_int32_t *out){
 413   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 414   vorbis_info_floor0 *info=look->vi;
 415
 416   if(memo){
 417     ogg_int32_t *lsp=(ogg_int32_t *)memo;
 418     ogg_int32_t amp=lsp[look->m];
 419
 420     /* take the coefficients back to a spectral envelope curve */
 421     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 422                         lsp,look->m,amp,info->ampdB,look->lsp_look);
 423     return(1);
 424   }
 425   memset(out,0,sizeof(*out)*look->n);
 426   return(0);
 427 }
 428
 429 /* export hooks */
 430 vorbis_func_floor floor0_exportbundle={
 431   &floor0_unpack,&floor0_look,&floor0_free_info,
 432   &floor0_free_look,&floor0_inverse1,&floor0_inverse2
 433 };
 434
 435