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.46 2001/10/02 00:14:30 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 bits=0;
 270   long val=0;
 271   static int seq=0;
 272
 273 #ifdef TRAIN_LSP
 274   FILE *of;
 275   FILE *ef;
 276   char buffer[80];
 277
 278 #if 1
 279   sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 280   of=fopen(buffer,"a");
 281 #endif
 282 #endif
 283
 284   seq++;
 285
 286
 287   /* our floor comes in on a [-Inf...0] dB scale.  The curve has to be
 288      positive, so we offset it. */
 289
 290   for(j=0;j<look->n;j++)
 291     codedflr[j]=logmask[j]+info->ampdB;
 292
 293   /* use 'out' as temp storage */
 294   /* Convert our floor to a set of lpc coefficients */
 295   amp=sqrt(_curve_to_lpc(codedflr,codedflr,look));
 296
 297   /* amp is in the range (0. to ampdB].  Encode that range using
 298      ampbits bits */
 299
 300   {
 301     long maxval=(1L<<info->ampbits)-1;
 302
 303     val=rint(amp/info->ampdB*maxval);
 304
 305     if(val<0)val=0;           /* likely */
 306     if(val>maxval)val=maxval; /* not bloody likely */
 307
 308     if(val>0)
 309       amp=(float)val/maxval*info->ampdB;
 310     else
 311       amp=0;
 312   }
 313
 314   if(val){
 315     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 316     _analysis_output("lpc",seq-1,codedflr,look->m,0,0);
 317     if(vorbis_lpc_to_lsp(codedflr,codedflr,look->m))
 318       val=0;
 319
 320   }
 321
 322   oggpack_write(&vb->opb,val,info->ampbits);
 323   look->bits+=info->ampbits+1;
 324   look->frames++;
 325
 326   if(val){
 327     float *lspwork=alloca(look->m*sizeof(*lspwork));
 328
 329     /* the spec supports using one of a number of codebooks.  Right
 330        now, encode using this lib supports only one */
 331     backend_lookup_state *be=vb->vd->backend_state;
 332     codebook *b;
 333     int booknum;
 334
 335     _analysis_output("lsp",seq-1,codedflr,look->m,0,0);
 336
 337     /* which codebook to use? We do it only by range right now. */
 338     if(info->numbooks>1){
 339       float last=0.;
 340       for(j=0;j<look->m;j++){
 341         float val=codedflr[j]-last;
 342         if(val<info->lessthan || val>info->greaterthan)break;
 343         last=codedflr[j];
 344       }
 345       if(j<look->m)
 346         booknum=0;
 347       else
 348         booknum=1;
 349     }else
 350       booknum=0;
 351
 352     b=be->fullbooks+info->books[booknum];
 353     oggpack_write(&vb->opb,booknum,_ilog(info->numbooks));
 354     look->bits+=_ilog(info->numbooks);
 355
 356 #ifdef TRAIN_LSP
 357     {
 358       float last=0.f;
 359       for(j=0;j<look->m;j++){
 360         fprintf(of,"%.12g, ",codedflr[j]-last);
 361         last=codedflr[j];
 362       }
 363     }
 364     fprintf(of,"\n");
 365     fclose(of);
 366
 367     sprintf(buffer,"lsp0ent_m%d_b%d.vqd",vb->mode,booknum);
 368     ef=fopen(buffer,"a");
 369
 370 #endif
 371
 372     /* code the spectral envelope, and keep track of the actual
 373        quantized values; we don't want creeping error as each block is
 374        nailed to the last quantized value of the previous block. */
 375
 376     for(j=0;j<look->m;j+=b->dim){
 377       int entry=_f0_fit(b,codedflr,lspwork,j);
 378       look->bits+=vorbis_book_encode(b,entry,&vb->opb);
 379
 380 #ifdef TRAIN_LSP
 381       fprintf(ef,"%d,\n",entry);
 382 #endif
 383
 384     }
 385
 386 #ifdef TRAIN_LSP
 387     fclose(ef);
 388 #endif
 389
 390     _analysis_output("lsp2",seq-1,lspwork,look->m,0,0);
 391
 392     /* take the coefficients back to a spectral envelope curve */
 393     for(j=0;j<look->n;j++)
 394       codedflr[j]=1.f;
 395     vorbis_lsp_to_curve(codedflr,look->linearmap,look->n,look->ln,
 396                         lspwork,look->m,amp,info->ampdB);
 397
 398     _analysis_output("barklsp",seq-1,codedflr,look->n,1,1);
 399     _analysis_output("lsp3",seq-1,codedflr,look->n,0,1);
 400
 401     return(val);
 402   }
 403
 404 #ifdef TRAIN_LSP
 405     fclose(of);
 406 #endif
 407
 408   memset(codedflr,0,sizeof(*codedflr)*look->n);
 409   memset(mdct,0,sizeof(*mdct)*look->n);
 410   return(val);
 411 }
 412
 413 static void *floor0_inverse1(vorbis_block *vb,vorbis_look_floor *i){
 414   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 415   vorbis_info_floor0 *info=look->vi;
 416   int j,k;
 417
 418   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 419   if(ampraw>0){ /* also handles the -1 out of data case */
 420     long maxval=(1<<info->ampbits)-1;
 421     float amp=(float)ampraw/maxval*info->ampdB;
 422     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 423
 424     if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
 425       backend_lookup_state *be=vb->vd->backend_state;
 426       codebook *b=be->fullbooks+info->books[booknum];
 427       float last=0.f;
 428       float *lsp=_vorbis_block_alloc(vb,sizeof(*lsp)*(look->m+1));
 429
 430       for(j=0;j<look->m;j+=b->dim)
 431         if(vorbis_book_decodev_set(b,lsp+j,&vb->opb,b->dim)==-1)goto eop;
 432       for(j=0;j<look->m;){
 433         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 434         last=lsp[j-1];
 435       }
 436
 437       lsp[look->m]=amp;
 438       return(lsp);
 439     }
 440   }
 441  eop:
 442   return(NULL);
 443 }
 444
 445 static int floor0_inverse2(vorbis_block *vb,vorbis_look_floor *i,
 446                            void *memo,float *out){
 447   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 448   vorbis_info_floor0 *info=look->vi;
 449
 450   if(memo){
 451     float *lsp=(float *)memo;
 452     float amp=lsp[look->m];
 453
 454     /* take the coefficients back to a spectral envelope curve */
 455     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 456                         lsp,look->m,amp,info->ampdB);
 457     return(1);
 458   }
 459   memset(out,0,sizeof(*out)*look->n);
 460   return(0);
 461 }
 462
 463 /* export hooks */
 464 vorbis_func_floor floor0_exportbundle={
 465   &floor0_pack,&floor0_unpack,&floor0_look,&floor0_copy_info,&floor0_free_info,
 466   &floor0_free_look,&floor0_forward,&floor0_inverse1,&floor0_inverse2
 467 };
 468
 469