lib/floor0.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: floor backend 0 implementation
  14  last mod: $Id: floor0.c,v 1.48 2001/12/20 01:00:26 segher Exp $
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <string.h>
  20 #include <math.h>
  21 #include <ogg/ogg.h>
  22 #include "vorbis/codec.h"
  23 #include "codec_internal.h"
  24 #include "registry.h"
  25 #include "lpc.h"
  26 #include "lsp.h"
  27 #include "codebook.h"
  28 #include "scales.h"
  29 #include "misc.h"
  30 #include "os.h"
  31
  32 #include "misc.h"
  33 #include <stdio.h>
  34
  35 typedef struct {
  36   long n;
  37   int ln;
  38   int  m;
  39   int *linearmap;
  40
  41   vorbis_info_floor0 *vi;
  42   lpc_lookup lpclook;
  43   float *lsp_look;
  44
  45   long bits;
  46   long frames;
  47 } vorbis_look_floor0;
  48
  49 /* infrastructure for finding fit */
  50 static long _f0_fit(codebook *book,
  51                     float *orig,
  52                     float *workfit,
  53                     int cursor){
  54   int dim=book->dim;
  55   float norm,base=0.f;
  56   int i,best=0;
  57   float *lsp=workfit+cursor;
  58
  59   if(cursor)base=workfit[cursor-1];
  60   norm=orig[cursor+dim-1]-base;
  61
  62   for(i=0;i<dim;i++)
  63     lsp[i]=(orig[i+cursor]-base);
  64   best=_best(book,lsp,1);
  65
  66   memcpy(lsp,book->valuelist+best*dim,dim*sizeof(*lsp));
  67   for(i=0;i<dim;i++)
  68     lsp[i]+=base;
  69   return(best);
  70 }
  71
  72 /***********************************************/
  73
  74 static vorbis_info_floor *floor0_copy_info (vorbis_info_floor *i){
  75   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
  76   vorbis_info_floor0 *ret=_ogg_malloc(sizeof(*ret));
  77   memcpy(ret,info,sizeof(*ret));
  78   return(ret);
  79 }
  80
  81 static void floor0_free_info(vorbis_info_floor *i){
  82   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
  83   if(info){
  84     memset(info,0,sizeof(*info));
  85     _ogg_free(info);
  86   }
  87 }
  88
  89 static void floor0_free_look(vorbis_look_floor *i){
  90   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
  91   if(look){
  92
  93     /*fprintf(stderr,"floor 0 bit usage %f\n",
  94       (float)look->bits/look->frames);*/
  95
  96     if(look->linearmap)_ogg_free(look->linearmap);
  97     if(look->lsp_look)_ogg_free(look->lsp_look);
  98     lpc_clear(&look->lpclook);
  99     memset(look,0,sizeof(*look));
 100     _ogg_free(look);
 101   }
 102 }
 103
 104 static void floor0_pack (vorbis_info_floor *i,oggpack_buffer *opb){
 105   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 106   int j;
 107   oggpack_write(opb,info->order,8);
 108   oggpack_write(opb,info->rate,16);
 109   oggpack_write(opb,info->barkmap,16);
 110   oggpack_write(opb,info->ampbits,6);
 111   oggpack_write(opb,info->ampdB,8);
 112   oggpack_write(opb,info->numbooks-1,4);
 113   for(j=0;j<info->numbooks;j++)
 114     oggpack_write(opb,info->books[j],8);
 115 }
 116
 117 static vorbis_info_floor *floor0_unpack (vorbis_info *vi,oggpack_buffer *opb){
 118   codec_setup_info     *ci=vi->codec_setup;
 119   int j;
 120
 121   vorbis_info_floor0 *info=_ogg_malloc(sizeof(*info));
 122   info->order=oggpack_read(opb,8);
 123   info->rate=oggpack_read(opb,16);
 124   info->barkmap=oggpack_read(opb,16);
 125   info->ampbits=oggpack_read(opb,6);
 126   info->ampdB=oggpack_read(opb,8);
 127   info->numbooks=oggpack_read(opb,4)+1;
 128
 129   if(info->order<1)goto err_out;
 130   if(info->rate<1)goto err_out;
 131   if(info->barkmap<1)goto err_out;
 132   if(info->numbooks<1)goto err_out;
 133
 134   for(j=0;j<info->numbooks;j++){
 135     info->books[j]=oggpack_read(opb,8);
 136     if(info->books[j]<0 || info->books[j]>=ci->books)goto err_out;
 137   }
 138   return(info);
 139
 140  err_out:
 141   floor0_free_info(info);
 142   return(NULL);
 143 }
 144
 145 /* initialize Bark scale and normalization lookups.  We could do this
 146    with static tables, but Vorbis allows a number of possible
 147    combinations, so it's best to do it computationally.
 148
 149    The below is authoritative in terms of defining scale mapping.
 150    Note that the scale depends on the sampling rate as well as the
 151    linear block and mapping sizes */
 152
 153 static vorbis_look_floor *floor0_look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 154                               vorbis_info_floor *i){
 155   int j;
 156   float scale;
 157   vorbis_info        *vi=vd->vi;
 158   codec_setup_info   *ci=vi->codec_setup;
 159   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 160   vorbis_look_floor0 *look=_ogg_calloc(1,sizeof(*look));
 161   look->m=info->order;
 162   look->n=ci->blocksizes[mi->blockflag]/2;
 163   look->ln=info->barkmap;
 164   look->vi=info;
 165
 166   if(vd->analysisp)
 167     lpc_init(&look->lpclook,look->ln,look->m);
 168
 169   /* we choose a scaling constant so that:
 170      floor(bark(rate/2-1)*C)=mapped-1
 171      floor(bark(rate/2)*C)=mapped */
 172   scale=look->ln/toBARK(info->rate/2.f);
 173
 174   /* the mapping from a linear scale to a smaller bark scale is
 175      straightforward.  We do *not* make sure that the linear mapping
 176      does not skip bark-scale bins; the decoder simply skips them and
 177      the encoder may do what it wishes in filling them.  They're
 178      necessary in some mapping combinations to keep the scale spacing
 179      accurate */
 180   look->linearmap=_ogg_malloc((look->n+1)*sizeof(*look->linearmap));
 181   for(j=0;j<look->n;j++){
 182     int val=floor( toBARK((info->rate/2.f)/look->n*j)
 183                    *scale); /* bark numbers represent band edges */
 184     if(val>=look->ln)val=look->ln; /* guard against the approximation */
 185     look->linearmap[j]=val;
 186   }
 187   look->linearmap[j]=-1;
 188
 189   look->lsp_look=_ogg_malloc(look->ln*sizeof(*look->lsp_look));
 190   for(j=0;j<look->ln;j++)
 191     look->lsp_look[j]=2*cos(M_PI/look->ln*j);
 192
 193   return look;
 194 }
 195
 196 /* less efficient than the decode side (written for clarity).  We're
 197    not bottlenecked here anyway */
 198
 199 float _curve_to_lpc(float *curve,float *lpc,
 200                      vorbis_look_floor0 *l){
 201   /* map the input curve to a bark-scale curve for encoding */
 202
 203   int mapped=l->ln;
 204   float *work=alloca(sizeof(*work)*mapped);
 205   int i,j,last=0;
 206   int bark=0;
 207   static int seq=0;
 208
 209   memset(work,0,sizeof(*work)*mapped);
 210
 211   /* Only the decode side is behavior-specced; for now in the encoder,
 212      we select the maximum value of each band as representative (this
 213      helps make sure peaks don't go out of range.  In error terms,
 214      selecting min would make more sense, but the codebook is trained
 215      numerically, so we don't actually lose.  We'd still want to
 216      use the original curve for error and noise estimation */
 217
 218   for(i=0;i<l->n;i++){
 219     bark=l->linearmap[i];
 220     if(work[bark]<curve[i])work[bark]=curve[i];
 221     if(bark>last+1){
 222       /* If the bark scale is climbing rapidly, some bins may end up
 223          going unused.  This isn't a waste actually; it keeps the
 224          scale resolution even so that the LPC generator has an easy
 225          time.  However, if we leave the bins empty we lose energy.
 226          So, fill 'em in.  The decoder does not do anything with  he
 227          unused bins, so we can fill them anyway we like to end up
 228          with a better spectral curve */
 229
 230       /* we'll always have a bin zero, so we don't need to guard init */
 231       long span=bark-last;
 232       for(j=1;j<span;j++){
 233         float del=(float)j/span;
 234         work[j+last]=work[bark]*del+work[last]*(1.f-del);
 235       }
 236     }
 237     last=bark;
 238   }
 239
 240   /* If we're over-ranged to avoid edge effects, fill in the end of spectrum gap */
 241   for(i=bark+1;i<mapped;i++)
 242     work[i]=work[i-1];
 243
 244
 245   /**********************/
 246
 247   for(i=0;i<l->n;i++)
 248     curve[i]-=150;
 249
 250   _analysis_output("barkfloor",seq,work,bark,0,0);
 251   _analysis_output("barkcurve",seq++,curve,l->n,1,0);
 252
 253   for(i=0;i<l->n;i++)
 254     curve[i]+=150;
 255
 256   /**********************/
 257
 258   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 259 }
 260
 261 static int floor0_forward(vorbis_block *vb,vorbis_look_floor *in,
 262                           float *mdct, const float *logmdct,   /* in */
 263                           const float *logmask, const float *logmax, /* in */
 264                           float *codedflr){          /* out */
 265   long j;
 266   vorbis_look_floor0 *look=(vorbis_look_floor0 *)in;
 267   vorbis_info_floor0 *info=look->vi;
 268   float amp;
 269   long val=0;
 270   static int seq=0;
 271
 272 #ifdef TRAIN_LSP
 273   FILE *of;
 274   FILE *ef;
 275   char buffer[80];
 276
 277 #if 1
 278   sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 279   of=fopen(buffer,"a");
 280 #endif
 281 #endif
 282
 283   seq++;
 284
 285
 286   /* our floor comes in on a [-Inf...0] dB scale.  The curve has to be
 287      positive, so we offset it. */
 288
 289   for(j=0;j<look->n;j++)
 290     codedflr[j]=logmask[j]+info->ampdB;
 291
 292   /* use 'out' as temp storage */
 293   /* Convert our floor to a set of lpc coefficients */
 294   amp=sqrt(_curve_to_lpc(codedflr,codedflr,look));
 295
 296   /* amp is in the range (0. to ampdB].  Encode that range using
 297      ampbits bits */
 298
 299   {
 300     long maxval=(1L<<info->ampbits)-1;
 301
 302     val=rint(amp/info->ampdB*maxval);
 303
 304     if(val<0)val=0;           /* likely */
 305     if(val>maxval)val=maxval; /* not bloody likely */
 306
 307     if(val>0)
 308       amp=(float)val/maxval*info->ampdB;
 309     else
 310       amp=0;
 311   }
 312
 313   if(val){
 314     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 315     _analysis_output("lpc",seq-1,codedflr,look->m,0,0);
 316     if(vorbis_lpc_to_lsp(codedflr,codedflr,look->m))
 317       val=0;
 318
 319   }
 320
 321   oggpack_write(&vb->opb,val,info->ampbits);
 322   look->bits+=info->ampbits+1;
 323   look->frames++;
 324
 325   if(val){
 326     float *lspwork=alloca(look->m*sizeof(*lspwork));
 327
 328     /* the spec supports using one of a number of codebooks.  Right
 329        now, encode using this lib supports only one */
 330     backend_lookup_state *be=vb->vd->backend_state;
 331     codebook *b;
 332     int booknum;
 333
 334     _analysis_output("lsp",seq-1,codedflr,look->m,0,0);
 335
 336     /* which codebook to use? We do it only by range right now. */
 337     if(info->numbooks>1){
 338       float last=0.;
 339       for(j=0;j<look->m;j++){
 340         float val=codedflr[j]-last;
 341         if(val<info->lessthan || val>info->greaterthan)break;
 342         last=codedflr[j];
 343       }
 344       if(j<look->m)
 345         booknum=0;
 346       else
 347         booknum=1;
 348     }else
 349       booknum=0;
 350
 351     b=be->fullbooks+info->books[booknum];
 352     oggpack_write(&vb->opb,booknum,_ilog(info->numbooks));
 353     look->bits+=_ilog(info->numbooks);
 354
 355 #ifdef TRAIN_LSP
 356     {
 357       float last=0.f;
 358       for(j=0;j<look->m;j++){
 359         fprintf(of,"%.12g, ",codedflr[j]-last);
 360         last=codedflr[j];
 361       }
 362     }
 363     fprintf(of,"\n");
 364     fclose(of);
 365
 366     sprintf(buffer,"lsp0ent_m%d_b%d.vqd",vb->mode,booknum);
 367     ef=fopen(buffer,"a");
 368
 369 #endif
 370
 371     /* code the spectral envelope, and keep track of the actual
 372        quantized values; we don't want creeping error as each block is
 373        nailed to the last quantized value of the previous block. */
 374
 375     for(j=0;j<look->m;j+=b->dim){
 376       int entry=_f0_fit(b,codedflr,lspwork,j);
 377       look->bits+=vorbis_book_encode(b,entry,&vb->opb);
 378
 379 #ifdef TRAIN_LSP
 380       fprintf(ef,"%d,\n",entry);
 381 #endif
 382
 383     }
 384
 385 #ifdef TRAIN_LSP
 386     fclose(ef);
 387 #endif
 388
 389     _analysis_output("lsp2",seq-1,lspwork,look->m,0,0);
 390
 391     /* take the coefficients back to a spectral envelope curve */
 392     for(j=0;j<look->n;j++)
 393       codedflr[j]=1.f;
 394     vorbis_lsp_to_curve(codedflr,look->linearmap,look->n,look->ln,
 395                         lspwork,look->m,amp,info->ampdB);
 396
 397     _analysis_output("barklsp",seq-1,codedflr,look->n,1,1);
 398     _analysis_output("lsp3",seq-1,codedflr,look->n,0,1);
 399
 400     return(val);
 401   }
 402
 403 #ifdef TRAIN_LSP
 404     fclose(of);
 405 #endif
 406
 407   memset(codedflr,0,sizeof(*codedflr)*look->n);
 408   memset(mdct,0,sizeof(*mdct)*look->n);
 409   return(val);
 410 }
 411
 412 static void *floor0_inverse1(vorbis_block *vb,vorbis_look_floor *i){
 413   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 414   vorbis_info_floor0 *info=look->vi;
 415   int j,k;
 416
 417   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 418   if(ampraw>0){ /* also handles the -1 out of data case */
 419     long maxval=(1<<info->ampbits)-1;
 420     float amp=(float)ampraw/maxval*info->ampdB;
 421     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 422
 423     if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
 424       backend_lookup_state *be=vb->vd->backend_state;
 425       codebook *b=be->fullbooks+info->books[booknum];
 426       float last=0.f;
 427       float *lsp=_vorbis_block_alloc(vb,sizeof(*lsp)*(look->m+1));
 428
 429       for(j=0;j<look->m;j+=b->dim)
 430         if(vorbis_book_decodev_set(b,lsp+j,&vb->opb,b->dim)==-1)goto eop;
 431       for(j=0;j<look->m;){
 432         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 433         last=lsp[j-1];
 434       }
 435
 436       lsp[look->m]=amp;
 437       return(lsp);
 438     }
 439   }
 440  eop:
 441   return(NULL);
 442 }
 443
 444 static int floor0_inverse2(vorbis_block *vb,vorbis_look_floor *i,
 445                            void *memo,float *out){
 446   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 447   vorbis_info_floor0 *info=look->vi;
 448
 449   if(memo){
 450     float *lsp=(float *)memo;
 451     float amp=lsp[look->m];
 452
 453     /* take the coefficients back to a spectral envelope curve */
 454     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 455                         lsp,look->m,amp,info->ampdB);
 456     return(1);
 457   }
 458   memset(out,0,sizeof(*out)*look->n);
 459   return(0);
 460 }
 461
 462 /* export hooks */
 463 vorbis_func_floor floor0_exportbundle={
 464   &floor0_pack,&floor0_unpack,&floor0_look,&floor0_copy_info,&floor0_free_info,
 465   &floor0_free_look,&floor0_forward,&floor0_inverse1,&floor0_inverse2
 466 };
 467
 468