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