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.33 2000/12/12 08:54:27 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   }
 181   look->linearmap[j]=-1;
 182
 183   look->lsp_look=_ogg_malloc(look->ln*sizeof(float));
 184   for(j=0;j<look->ln;j++)
 185     look->lsp_look[j]=2*cos(M_PI/look->ln*j);
 186
 187   return look;
 188 }
 189
 190 /* less efficient than the decode side (written for clarity).  We're
 191    not bottlenecked here anyway */
 192
 193 float _curve_to_lpc(float *curve,float *lpc,
 194                      vorbis_look_floor0 *l){
 195   /* map the input curve to a bark-scale curve for encoding */
 196
 197   int mapped=l->ln;
 198   float *work=alloca(sizeof(float)*mapped);
 199   int i,j,last=0;
 200   int bark=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.-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   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 238 }
 239
 240 /* generate the whole freq response curve of an LSP IIR filter */
 241 static int floor0_forward(vorbis_block *vb,vorbis_look_floor *i,
 242                     float *in,float *out,vorbis_bitbuffer *vbb){
 243   long j;
 244   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 245   vorbis_info_floor0 *info=look->vi;
 246   float *work=alloca((look->ln+look->n)*sizeof(float));
 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
 262   sprintf(buffer,"lsp0ent_%d.vqd",vb->mode);
 263   ef=fopen(buffer,"a");
 264 #endif
 265
 266   /* our floor comes in on a linear scale; go to a [-Inf...0] dB
 267      scale.  The curve has to be positive, so we offset it. */
 268
 269   for(j=0;j<look->n;j++)
 270     work[j]=todB(in[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(work,out,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     /*oggpack_write(&vb->opb,val,info->ampbits);*/
 288     if(val>0)
 289       amp=(float)val/maxval*info->ampdB;
 290     else
 291       amp=0;
 292   }
 293
 294   if(val){
 295
 296     /* the spec supports using one of a number of codebooks.  Right
 297        now, encode using this lib supports only one */
 298     backend_lookup_state *be=vb->vd->backend_state;
 299     codebook *b=be->fullbooks+info->books[0];
 300     bitbuf_write(vbb,0,_ilog(info->numbooks));
 301
 302     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 303     vorbis_lpc_to_lsp(out,out,look->m);
 304
 305 #ifdef ANALYSIS
 306     {
 307       float *lspwork=alloca(look->m*sizeof(float));
 308       memcpy(lspwork,out,look->m*sizeof(float));
 309       vorbis_lsp_to_curve(work,look->linearmap,look->n,look->ln,
 310                           lspwork,look->m,amp,info->ampdB);
 311       _analysis_output("prefit",seq,work,look->n,0,1);
 312
 313     }
 314
 315 #endif
 316
 317
 318 #if 1
 319 #ifdef TRAIN_LSP
 320     {
 321       float last=0.;
 322       for(j=0;j<look->m;j++){
 323         fprintf(of,"%.12g, ",out[j]-last);
 324         last=out[j];
 325       }
 326     }
 327     fprintf(of,"\n");
 328     fclose(of);
 329 #endif
 330 #endif
 331
 332     /* code the spectral envelope, and keep track of the actual
 333        quantized values; we don't want creeping error as each block is
 334        nailed to the last quantized value of the previous block. */
 335
 336     for(j=0;j<look->m;j+=b->dim){
 337       int entry=_f0_fit(b,out,work,j);
 338       bits+=vorbis_book_bufencode(b,entry,vbb);
 339
 340 #ifdef TRAIN_LSP
 341       fprintf(ef,"%d,\n",entry);
 342 #endif
 343
 344     }
 345
 346 #ifdef ANALYSIS
 347     {
 348       float last=0;
 349       for(j=0;j<look->m;j++){
 350         out[j]=work[j]-last;
 351         last=work[j];
 352       }
 353     }
 354
 355 #endif
 356
 357 #ifdef TRAIN_LSP
 358     fclose(ef);
 359 #endif
 360
 361     /* take the coefficients back to a spectral envelope curve */
 362     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 363                         work,look->m,amp,info->ampdB);
 364     return(val);
 365   }
 366
 367   memset(out,0,sizeof(float)*look->n);
 368   seq++;
 369   return(val);
 370 }
 371
 372 static float floor0_forward2(vorbis_block *vb,vorbis_look_floor *i,
 373                           long amp,float error,
 374                           vorbis_bitbuffer *vbb){
 375
 376   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 377   vorbis_info_floor0 *info=look->vi;
 378   if(amp){
 379     long maxval=(1L<<info->ampbits)-1;
 380     long adj=rint(todB(error)/info->ampdB*maxval/2);
 381
 382     amp+=adj;
 383     if(amp<1)amp=1;
 384
 385     oggpack_write(&vb->opb,amp,info->ampbits);
 386     bitbuf_pack(&vb->opb,vbb);
 387     return(fromdB((float)adj/maxval*info->ampdB));
 388   }else{
 389     oggpack_write(&vb->opb,0,info->ampbits);
 390     bitbuf_pack(&vb->opb,vbb);
 391   }
 392   return(0.);
 393 }
 394
 395
 396 static int floor0_inverse(vorbis_block *vb,vorbis_look_floor *i,float *out){
 397   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 398   vorbis_info_floor0 *info=look->vi;
 399   int j,k;
 400
 401   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 402   if(ampraw>0){ /* also handles the -1 out of data case */
 403     long maxval=(1<<info->ampbits)-1;
 404     float amp=(float)ampraw/maxval*info->ampdB;
 405     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 406     float *lsp=alloca(sizeof(float)*look->m);
 407
 408     if(booknum!=-1){
 409       backend_lookup_state *be=vb->vd->backend_state;
 410       codebook *b=be->fullbooks+info->books[booknum];
 411       float last=0.;
 412
 413       memset(out,0,sizeof(float)*look->m);
 414
 415       for(j=0;j<look->m;j+=b->dim)
 416         if(vorbis_book_decodevs(b,lsp+j,&vb->opb,1,-1)==-1)goto eop;
 417       for(j=0;j<look->m;){
 418         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 419         last=lsp[j-1];
 420       }
 421
 422       /* take the coefficients back to a spectral envelope curve */
 423       vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 424                           lsp,look->m,amp,info->ampdB);
 425       return(1);
 426     }
 427   }
 428
 429  eop:
 430   memset(out,0,sizeof(float)*look->n);
 431   return(0);
 432 }
 433
 434 /* export hooks */
 435 vorbis_func_floor floor0_exportbundle={
 436   &floor0_pack,&floor0_unpack,&floor0_look,&floor0_copy_info,&floor0_free_info,
 437   &floor0_free_look,&floor0_forward,&floor0_forward2,&floor0_inverse
 438 };
 439
 440