lib/floor0.c

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