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