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