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