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.45 2001/08/13 01:36:56 xiphmont 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(float));
  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(vorbis_info_floor0));
  77   memcpy(ret,info,sizeof(vorbis_info_floor0));
  78   return(ret);
  79 }
  80
  81 static void floor0_free_info(vorbis_info_floor *i){
  82   if(i){
  83     memset(i,0,sizeof(vorbis_info_floor0));
  84     _ogg_free(i);
  85   }
  86 }
  87
  88 static void floor0_free_look(vorbis_look_floor *i){
  89   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
  90   if(i){
  91
  92     /*fprintf(stderr,"floor 0 bit usage %f\n",
  93       (float)look->bits/look->frames);*/
  94
  95     if(look->linearmap)_ogg_free(look->linearmap);
  96     if(look->lsp_look)_ogg_free(look->lsp_look);
  97     lpc_clear(&look->lpclook);
  98     memset(look,0,sizeof(vorbis_look_floor0));
  99     _ogg_free(look);
 100   }
 101 }
 102
 103 static void floor0_pack (vorbis_info_floor *i,oggpack_buffer *opb){
 104   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 105   int j;
 106   oggpack_write(opb,info->order,8);
 107   oggpack_write(opb,info->rate,16);
 108   oggpack_write(opb,info->barkmap,16);
 109   oggpack_write(opb,info->ampbits,6);
 110   oggpack_write(opb,info->ampdB,8);
 111   oggpack_write(opb,info->numbooks-1,4);
 112   for(j=0;j<info->numbooks;j++)
 113     oggpack_write(opb,info->books[j],8);
 114 }
 115
 116 static vorbis_info_floor *floor0_unpack (vorbis_info *vi,oggpack_buffer *opb){
 117   codec_setup_info     *ci=vi->codec_setup;
 118   int j;
 119
 120   vorbis_info_floor0 *info=_ogg_malloc(sizeof(vorbis_info_floor0));
 121   info->order=oggpack_read(opb,8);
 122   info->rate=oggpack_read(opb,16);
 123   info->barkmap=oggpack_read(opb,16);
 124   info->ampbits=oggpack_read(opb,6);
 125   info->ampdB=oggpack_read(opb,8);
 126   info->numbooks=oggpack_read(opb,4)+1;
 127
 128   if(info->order<1)goto err_out;
 129   if(info->rate<1)goto err_out;
 130   if(info->barkmap<1)goto err_out;
 131   if(info->numbooks<1)goto err_out;
 132
 133   for(j=0;j<info->numbooks;j++){
 134     info->books[j]=oggpack_read(opb,8);
 135     if(info->books[j]<0 || info->books[j]>=ci->books)goto err_out;
 136   }
 137   return(info);
 138
 139  err_out:
 140   floor0_free_info(info);
 141   return(NULL);
 142 }
 143
 144 /* initialize Bark scale and normalization lookups.  We could do this
 145    with static tables, but Vorbis allows a number of possible
 146    combinations, so it's best to do it computationally.
 147
 148    The below is authoritative in terms of defining scale mapping.
 149    Note that the scale depends on the sampling rate as well as the
 150    linear block and mapping sizes */
 151
 152 static vorbis_look_floor *floor0_look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 153                               vorbis_info_floor *i){
 154   int j;
 155   float scale;
 156   vorbis_info        *vi=vd->vi;
 157   codec_setup_info   *ci=vi->codec_setup;
 158   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 159   vorbis_look_floor0 *look=_ogg_calloc(1,sizeof(vorbis_look_floor0));
 160   look->m=info->order;
 161   look->n=ci->blocksizes[mi->blockflag]/2;
 162   look->ln=info->barkmap;
 163   look->vi=info;
 164
 165   if(vd->analysisp)
 166     lpc_init(&look->lpclook,look->ln,look->m);
 167
 168   /* we choose a scaling constant so that:
 169      floor(bark(rate/2-1)*C)=mapped-1
 170      floor(bark(rate/2)*C)=mapped */
 171   scale=look->ln/toBARK(info->rate/2.f);
 172
 173   /* the mapping from a linear scale to a smaller bark scale is
 174      straightforward.  We do *not* make sure that the linear mapping
 175      does not skip bark-scale bins; the decoder simply skips them and
 176      the encoder may do what it wishes in filling them.  They're
 177      necessary in some mapping combinations to keep the scale spacing
 178      accurate */
 179   look->linearmap=_ogg_malloc((look->n+1)*sizeof(int));
 180   for(j=0;j<look->n;j++){
 181     int val=floor( toBARK((info->rate/2.f)/look->n*j)
 182                    *scale); /* bark numbers represent band edges */
 183     if(val>=look->ln)val=look->ln; /* guard against the approximation */
 184     look->linearmap[j]=val;
 185   }
 186   look->linearmap[j]=-1;
 187
 188   look->lsp_look=_ogg_malloc(look->ln*sizeof(float));
 189   for(j=0;j<look->ln;j++)
 190     look->lsp_look[j]=2*cos(M_PI/look->ln*j);
 191
 192   return look;
 193 }
 194
 195 /* less efficient than the decode side (written for clarity).  We're
 196    not bottlenecked here anyway */
 197
 198 float _curve_to_lpc(float *curve,float *lpc,
 199                      vorbis_look_floor0 *l){
 200   /* map the input curve to a bark-scale curve for encoding */
 201
 202   int mapped=l->ln;
 203   float *work=alloca(sizeof(float)*mapped);
 204   int i,j,last=0;
 205   int bark=0;
 206   static int seq=0;
 207
 208   memset(work,0,sizeof(float)*mapped);
 209
 210   /* Only the decode side is behavior-specced; for now in the encoder,
 211      we select the maximum value of each band as representative (this
 212      helps make sure peaks don't go out of range.  In error terms,
 213      selecting min would make more sense, but the codebook is trained
 214      numerically, so we don't actually lose.  We'd still want to
 215      use the original curve for error and noise estimation */
 216
 217   for(i=0;i<l->n;i++){
 218     bark=l->linearmap[i];
 219     if(work[bark]<curve[i])work[bark]=curve[i];
 220     if(bark>last+1){
 221       /* If the bark scale is climbing rapidly, some bins may end up
 222          going unused.  This isn't a waste actually; it keeps the
 223          scale resolution even so that the LPC generator has an easy
 224          time.  However, if we leave the bins empty we lose energy.
 225          So, fill 'em in.  The decoder does not do anything with  he
 226          unused bins, so we can fill them anyway we like to end up
 227          with a better spectral curve */
 228
 229       /* we'll always have a bin zero, so we don't need to guard init */
 230       long span=bark-last;
 231       for(j=1;j<span;j++){
 232         float del=(float)j/span;
 233         work[j+last]=work[bark]*del+work[last]*(1.f-del);
 234       }
 235     }
 236     last=bark;
 237   }
 238
 239   /* If we're over-ranged to avoid edge effects, fill in the end of spectrum gap */
 240   for(i=bark+1;i<mapped;i++)
 241     work[i]=work[i-1];
 242
 243
 244   /**********************/
 245
 246   for(i=0;i<l->n;i++)
 247     curve[i]-=150;
 248
 249   _analysis_output("barkfloor",seq,work,bark,0,0);
 250   _analysis_output("barkcurve",seq++,curve,l->n,1,0);
 251
 252   for(i=0;i<l->n;i++)
 253     curve[i]+=150;
 254
 255   /**********************/
 256
 257   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 258 }
 259
 260 static int floor0_forward(vorbis_block *vb,vorbis_look_floor *in,
 261                           float *mdct, const float *logmdct,   /* in */
 262                           const float *logmask, const float *logmax, /* in */
 263                           float *codedflr){          /* out */
 264   long j;
 265   vorbis_look_floor0 *look=(vorbis_look_floor0 *)in;
 266   vorbis_info_floor0 *info=look->vi;
 267   float amp;
 268   long bits=0;
 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(float));
 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(float)*look->n);
 408   memset(mdct,0,sizeof(float)*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(float)*(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(float)*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