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.24 2000/08/31 08:01:34 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.;
  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 LSP 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     work[j]=todB(in[j])+info->ampdB;
 289
 290   /* use 'out' as temp storage */
 291   /* Convert our floor to a set of lpc coefficients */
 292   amp=sqrt(_curve_to_lpc(work,out,look,seq));
 293
 294   /* amp is in the range (0. to ampdB].  Encode that range using
 295      ampbits bits */
 296
 297   {
 298     long maxval=(1L<<info->ampbits)-1;
 299
 300     long val=rint(amp/info->ampdB*maxval);
 301
 302     if(val<0)val=0;           /* likely */
 303     if(val>maxval)val=maxval; /* not bloody likely */
 304
 305     _oggpack_write(&vb->opb,val,info->ampbits);
 306     if(val>0)
 307       amp=(float)val/maxval*info->ampdB;
 308     else
 309       amp=0;
 310   }
 311
 312   if(amp>0){
 313
 314     /* the spec supports using one of a number of codebooks.  Right
 315        now, encode using this lib supports only one */
 316     codebook *b=vb->vd->fullbooks+info->books[0];
 317     _oggpack_write(&vb->opb,0,_ilog(info->numbooks));
 318
 319     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 320     vorbis_lpc_to_lsp(out,out,look->m);
 321
 322 #if 1
 323 #ifdef TRAIN_LSP
 324     {
 325       double last=0.;
 326       for(j=0;j<look->m;j++){
 327         fprintf(of,"%.12g, ",out[j]-last);
 328         last=out[j];
 329       }
 330     }
 331     fprintf(of,"\n");
 332     fclose(of);
 333 #endif
 334 #endif
 335
 336     /* code the spectral envelope, and keep track of the actual
 337        quantized values; we don't want creeping error as each block is
 338        nailed to the last quantized value of the previous block. */
 339
 340     for(j=0;j<look->m;j+=b->dim){
 341       int entry=_f0_fit(b,out,work,j);
 342       bits+=vorbis_book_encode(b,entry,&vb->opb);
 343
 344 #ifdef TRAIN_LSP
 345       fprintf(ef,"%d,\n",entry);
 346 #endif
 347
 348     }
 349
 350 #ifdef ANALYSIS
 351     {
 352       double last=0;
 353       for(j=0;j<look->m;j++){
 354         out[j]=work[j]-last;
 355         last=work[j];
 356       }
 357     }
 358     _analysis_output("lsp",seq,out,look->m,0,0);
 359
 360 #endif
 361
 362 #ifdef TRAIN_LSP
 363     fclose(ef);
 364 #endif
 365
 366     /* take the coefficients back to a spectral envelope curve */
 367     _lsp_to_curve(out,work,amp,look,"Ffloor",seq++);
 368     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
 369     return(1);
 370   }
 371
 372   memset(out,0,sizeof(double)*look->n);
 373   seq++;
 374   return(0);
 375 }
 376
 377 static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
 378   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 379   vorbis_info_floor0 *info=look->vi;
 380   int j,k;
 381
 382   int ampraw=_oggpack_read(&vb->opb,info->ampbits);
 383   if(ampraw>0){ /* also handles the -1 out of data case */
 384     long maxval=(1<<info->ampbits)-1;
 385     double amp=(float)ampraw/maxval*info->ampdB;
 386     int booknum=_oggpack_read(&vb->opb,_ilog(info->numbooks));
 387
 388     if(booknum!=-1){
 389       codebook *b=vb->vd->fullbooks+info->books[booknum];
 390       double last=0.;
 391
 392       memset(out,0,sizeof(double)*look->m);
 393
 394       for(j=0;j<look->m;j+=b->dim)
 395         if(vorbis_book_decodevs(b,out+j,&vb->opb,1,-1)==-1)goto eop;
 396       for(j=0;j<look->m;){
 397         for(k=0;k<b->dim;k++,j++)out[j]+=last;
 398         last=out[j-1];
 399       }
 400
 401       /* take the coefficients back to a spectral envelope curve */
 402       _lsp_to_curve(out,out,amp,look,"",0);
 403
 404       for(j=0;j<look->n;j++)out[j]=fromdB(out[j]-info->ampdB);
 405       return(1);
 406     }
 407   }
 408
 409  eop:
 410   memset(out,0,sizeof(double)*look->n);
 411   return(0);
 412 }
 413
 414 /* export hooks */
 415 vorbis_func_floor floor0_exportbundle={
 416   &pack,&unpack,&look,&free_info,&free_look,&forward,&inverse
 417 };
 418
 419