lib/floor0.c

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