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