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.43 2001/06/15 23:31:00 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   static int seq=0;
 201
 202   memset(work,0,sizeof(float)*mapped);
 203
 204   /* Only the decode side is behavior-specced; for now in the encoder,
 205      we select the maximum value of each band as representative (this
 206      helps make sure peaks don't go out of range.  In error terms,
 207      selecting min would make more sense, but the codebook is trained
 208      numerically, so we don't actually lose.  We'd still want to
 209      use the original curve for error and noise estimation */
 210
 211   for(i=0;i<l->n;i++){
 212     bark=l->linearmap[i];
 213     if(work[bark]<curve[i])work[bark]=curve[i];
 214     if(bark>last+1){
 215       /* If the bark scale is climbing rapidly, some bins may end up
 216          going unused.  This isn't a waste actually; it keeps the
 217          scale resolution even so that the LPC generator has an easy
 218          time.  However, if we leave the bins empty we lose energy.
 219          So, fill 'em in.  The decoder does not do anything with  he
 220          unused bins, so we can fill them anyway we like to end up
 221          with a better spectral curve */
 222
 223       /* we'll always have a bin zero, so we don't need to guard init */
 224       long span=bark-last;
 225       for(j=1;j<span;j++){
 226         float del=(float)j/span;
 227         work[j+last]=work[bark]*del+work[last]*(1.f-del);
 228       }
 229     }
 230     last=bark;
 231   }
 232
 233   /* If we're over-ranged to avoid edge effects, fill in the end of spectrum gap */
 234   for(i=bark+1;i<mapped;i++)
 235     work[i]=work[i-1];
 236
 237
 238   /**********************/
 239
 240   for(i=0;i<l->n;i++)
 241     curve[i]-=150;
 242
 243   _analysis_output("barkfloor",seq,work,bark,0,0);
 244   _analysis_output("barkcurve",seq++,curve,l->n,1,0);
 245
 246   for(i=0;i<l->n;i++)
 247     curve[i]+=150;
 248
 249   /**********************/
 250
 251   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 252 }
 253
 254 static int floor0_forward(vorbis_block *vb,vorbis_look_floor *in,
 255                           const float *mdct, const float *logmdct,   /* in */
 256                           const float *logmask, const float *logmax, /* in */
 257                           float *residue, float *codedflr){          /* out */
 258   long j;
 259   vorbis_look_floor0 *look=(vorbis_look_floor0 *)in;
 260   vorbis_info_floor0 *info=look->vi;
 261   float amp;
 262   long bits=0;
 263   long val=0;
 264   static int seq=0;
 265
 266 #ifdef TRAIN_LSP
 267   FILE *of;
 268   FILE *ef;
 269   char buffer[80];
 270
 271 #if 1
 272   sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 273   of=fopen(buffer,"a");
 274 #endif
 275 #endif
 276
 277   seq++;
 278
 279
 280   /* our floor comes in on a [-Inf...0] dB scale.  The curve has to be
 281      positive, so we offset it. */
 282
 283   for(j=0;j<look->n;j++)
 284     codedflr[j]=logmask[j]+info->ampdB;
 285
 286   /* use 'out' as temp storage */
 287   /* Convert our floor to a set of lpc coefficients */
 288   amp=sqrt(_curve_to_lpc(codedflr,codedflr,look));
 289
 290   /* amp is in the range (0. to ampdB].  Encode that range using
 291      ampbits bits */
 292
 293   {
 294     long maxval=(1L<<info->ampbits)-1;
 295
 296     val=rint(amp/info->ampdB*maxval);
 297
 298     if(val<0)val=0;           /* likely */
 299     if(val>maxval)val=maxval; /* not bloody likely */
 300
 301     if(val>0)
 302       amp=(float)val/maxval*info->ampdB;
 303     else
 304       amp=0;
 305   }
 306
 307   if(val){
 308     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 309     _analysis_output("lpc",seq-1,codedflr,look->m,0,0);
 310     if(vorbis_lpc_to_lsp(codedflr,codedflr,look->m))
 311       val=0;
 312
 313   }
 314
 315   oggpack_write(&vb->opb,val,info->ampbits);
 316
 317   if(val){
 318     float *lspwork=alloca(look->m*sizeof(float));
 319
 320     /* the spec supports using one of a number of codebooks.  Right
 321        now, encode using this lib supports only one */
 322     backend_lookup_state *be=vb->vd->backend_state;
 323     codebook *b;
 324     int booknum;
 325
 326     _analysis_output("lsp",seq-1,codedflr,look->m,0,0);
 327
 328     /* which codebook to use? We do it only by range right now. */
 329     if(info->numbooks>1){
 330       float last=0.;
 331       for(j=0;j<look->m;j++){
 332         float val=codedflr[j]-last;
 333         if(val<info->lessthan || val>info->greaterthan)break;
 334         last=codedflr[j];
 335       }
 336       if(j<look->m)
 337         booknum=0;
 338       else
 339         booknum=1;
 340     }else
 341       booknum=0;
 342
 343     b=be->fullbooks+info->books[booknum];
 344     oggpack_write(&vb->opb,booknum,_ilog(info->numbooks));
 345
 346
 347 #ifdef TRAIN_LSP
 348     {
 349       float last=0.f;
 350       for(j=0;j<look->m;j++){
 351         fprintf(of,"%.12g, ",codedflr[j]-last);
 352         last=codedflr[j];
 353       }
 354     }
 355     fprintf(of,"\n");
 356     fclose(of);
 357
 358     sprintf(buffer,"lsp0ent_m%d_b%d.vqd",vb->mode,booknum);
 359     ef=fopen(buffer,"a");
 360
 361 #endif
 362
 363     /* code the spectral envelope, and keep track of the actual
 364        quantized values; we don't want creeping error as each block is
 365        nailed to the last quantized value of the previous block. */
 366
 367     for(j=0;j<look->m;j+=b->dim){
 368       int entry=_f0_fit(b,codedflr,lspwork,j);
 369       bits+=vorbis_book_encode(b,entry,&vb->opb);
 370
 371 #ifdef TRAIN_LSP
 372       fprintf(ef,"%d,\n",entry);
 373 #endif
 374
 375     }
 376
 377 #ifdef TRAIN_LSP
 378     fclose(ef);
 379 #endif
 380
 381     _analysis_output("lsp2",seq-1,lspwork,look->m,0,0);
 382
 383     /* take the coefficients back to a spectral envelope curve */
 384     vorbis_lsp_to_curve(codedflr,look->linearmap,look->n,look->ln,
 385                         lspwork,look->m,amp,info->ampdB);
 386
 387     _analysis_output("barklsp",seq-1,codedflr,look->n,1,1);
 388     _analysis_output("lsp3",seq-1,codedflr,look->n,0,1);
 389
 390     /* generate residue output */
 391     for(j=0;j<look->n;j++)
 392       residue[j]=mdct[j]/codedflr[j];
 393
 394     return(val);
 395   }
 396
 397 #ifdef TRAIN_LSP
 398     fclose(of);
 399 #endif
 400
 401   memset(codedflr,0,sizeof(float)*look->n);
 402   memset(residue,0,sizeof(float)*look->n);
 403   return(val);
 404 }
 405
 406 static void *floor0_inverse1(vorbis_block *vb,vorbis_look_floor *i){
 407   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 408   vorbis_info_floor0 *info=look->vi;
 409   int j,k;
 410
 411   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 412   if(ampraw>0){ /* also handles the -1 out of data case */
 413     long maxval=(1<<info->ampbits)-1;
 414     float amp=(float)ampraw/maxval*info->ampdB;
 415     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 416
 417     if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
 418       backend_lookup_state *be=vb->vd->backend_state;
 419       codebook *b=be->fullbooks+info->books[booknum];
 420       float last=0.f;
 421       float *lsp=_vorbis_block_alloc(vb,sizeof(float)*(look->m+1));
 422
 423       for(j=0;j<look->m;j+=b->dim)
 424         if(vorbis_book_decodev_set(b,lsp+j,&vb->opb,b->dim)==-1)goto eop;
 425       for(j=0;j<look->m;){
 426         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 427         last=lsp[j-1];
 428       }
 429
 430       lsp[look->m]=amp;
 431       return(lsp);
 432     }
 433   }
 434  eop:
 435   return(NULL);
 436 }
 437
 438 static int floor0_inverse2(vorbis_block *vb,vorbis_look_floor *i,
 439                            void *memo,float *out){
 440   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 441   vorbis_info_floor0 *info=look->vi;
 442
 443   if(memo){
 444     float *lsp=(float *)memo;
 445     float amp=lsp[look->m];
 446
 447     /* take the coefficients back to a spectral envelope curve */
 448     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 449                         lsp,look->m,amp,info->ampdB);
 450     return(1);
 451   }
 452   memset(out,0,sizeof(float)*look->n);
 453   return(0);
 454 }
 455
 456 /* export hooks */
 457 vorbis_func_floor floor0_exportbundle={
 458   &floor0_pack,&floor0_unpack,&floor0_look,&floor0_copy_info,&floor0_free_info,
 459   &floor0_free_look,&floor0_forward,&floor0_inverse1,&floor0_inverse2
 460 };
 461
 462