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.18 2000/07/12 09:36:17 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 /* error relationship from coefficient->curve is nonlinear, so fit
  37    using curve lookups, not the coefficients */
  38 typedef struct {
  39   int usage;
  40   long *entrymap;
  41   double **fits;
  42   int fitsize;
  43   codebook *book;
  44   lpc_lookup lpclook;
  45 } LSP_fit_lookup;
  46
  47 typedef struct {
  48   long n;
  49   int ln;
  50   int  m;
  51   int *linearmap;
  52
  53   vorbis_info_floor0 *vi;
  54   lpc_lookup lpclook;
  55
  56   LSP_fit_lookup *encodefit;
  57 } vorbis_look_floor0;
  58
  59 /* infrastructure for setting up, finding and encoding fit */
  60 static void _f0_fit_init(LSP_fit_lookup *f,codebook *b){
  61   int i,j,usage=0;
  62   double max=0;
  63   int dim=b->dim;
  64   double *work;
  65   double *lsp;
  66
  67   memset(f,0,sizeof(LSP_fit_lookup));
  68
  69   /* count actual codeword usage in case the book is sparse */
  70   for(i=0;i<b->entries;i++)
  71     if(b->c->lengthlist[i]>0){
  72       if(b->valuelist[(i+1)*b->dim-1]>max)
  73         max=b->valuelist[(i+1)*b->dim-1];
  74       usage++;
  75     }
  76
  77   /* allocate memory */
  78   f->usage=usage;
  79   f->entrymap=malloc(usage*sizeof(long));
  80   f->fits=malloc(usage*sizeof(double *));
  81   f->fitsize=16;
  82   f->book=b;
  83   lpc_init(&f->lpclook,f->fitsize,dim-1);
  84
  85   usage=0;
  86   work=alloca(f->fitsize*2*sizeof(double));
  87   lsp=alloca(dim*sizeof(double));
  88   for(i=0;i<b->entries;i++)
  89     if(b->c->lengthlist[i]>0){
  90       double *orig=b->valuelist+i*b->dim;
  91       double norm=orig[b->dim-1];
  92       f->fits[usage]=malloc(f->fitsize*sizeof(double));
  93
  94       for(j=0;j<dim-1;j++)
  95         lsp[j]=orig[j]/norm*M_PI;
  96       vorbis_lsp_to_lpc(lsp,lsp,dim-1);
  97       vorbis_lpc_to_curve(work,lsp,norm*norm,&f->lpclook);
  98       memcpy(f->fits[usage],work,f->fitsize*sizeof(double));
  99
 100       f->entrymap[usage]=i;
 101
 102       usage++;
 103     }
 104 }
 105
 106 static void _f0_fit_clear(LSP_fit_lookup *f){
 107   if(f){
 108     int i;
 109     for(i=0;i<f->usage;i++)
 110       free(f->fits[i]);
 111     free(f->fits);
 112     free(f->entrymap);
 113     lpc_clear(&f->lpclook);
 114     memset(f,0,sizeof(LSP_fit_lookup));
 115   }
 116 }
 117
 118 double _curve_error1(double *curve1,double *curve2,long n){
 119   double acc=0.;
 120   long i;
 121   for(i=0;i<n;i++){
 122     double val=curve1[i]-curve2[i];
 123     acc+=val*val;
 124   }
 125   return(acc);
 126 }
 127
 128 double _curve_error2(double *curve1,double *curve2,long n,double max){
 129   double acc=0.;
 130   long i;
 131   for(i=0;i<n;i++){
 132     double val=curve1[i]-curve2[i];
 133     acc+=val*val;
 134     if(acc>max)return(acc);
 135   }
 136   return(acc);
 137 }
 138
 139 static long _f0_fit(LSP_fit_lookup *f,
 140                     codebook *book,
 141                     double *orig,
 142                     double *workfit,
 143                     int cursor){
 144   int dim=book->dim;
 145   double norm,base=0.,err=0.;
 146   int i,best=0;
 147   double *lsp=alloca(dim*sizeof(double));
 148   double *ref=alloca(f->fitsize*2*sizeof(double));
 149
 150   if(f->usage<=0)return(-1);
 151
 152   /* gen a curve for fitting */
 153   if(cursor)base=workfit[cursor-1];
 154   norm=orig[cursor+dim-1]-base;
 155   for(i=0;i<dim-1;i++)
 156     lsp[i]=(orig[i+cursor]-base)/norm*M_PI;
 157   vorbis_lsp_to_lpc(lsp,lsp,dim-1);
 158   vorbis_lpc_to_curve(ref,lsp,norm*norm,&f->lpclook);
 159
 160   err=_curve_error1(f->fits[0],ref,f->fitsize);
 161   for(i=1;i<f->usage;i++){
 162     double this=_curve_error2(f->fits[i],ref,f->fitsize,err);
 163     if(this<err){
 164       err=this;
 165       best=i;
 166     }
 167   }
 168
 169   /* choose the best fit from the tests and set workfit to it */
 170   best=f->entrymap[best];
 171   memcpy(workfit+cursor,book->valuelist+best*dim,dim*sizeof(double));
 172   for(i=0;i<dim;i++)
 173     workfit[i+cursor]+=base;
 174   return(best);
 175 }
 176
 177 /***********************************************/
 178
 179 static void free_info(vorbis_info_floor *i){
 180   if(i){
 181     memset(i,0,sizeof(vorbis_info_floor0));
 182     free(i);
 183   }
 184 }
 185
 186 static void free_look(vorbis_look_floor *i){
 187   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 188   if(i){
 189     int j;
 190     for(j=0;j<look->vi->numbooks;j++)
 191       _f0_fit_clear(look->encodefit+j);
 192     free(look->encodefit);
 193     if(look->linearmap)free(look->linearmap);
 194     lpc_clear(&look->lpclook);
 195     memset(look,0,sizeof(vorbis_look_floor0));
 196     free(look);
 197   }
 198 }
 199
 200 static void pack (vorbis_info_floor *i,oggpack_buffer *opb){
 201   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 202   int j;
 203   _oggpack_write(opb,info->order,8);
 204   _oggpack_write(opb,info->rate,16);
 205   _oggpack_write(opb,info->barkmap,16);
 206   _oggpack_write(opb,info->ampbits,6);
 207   _oggpack_write(opb,info->ampdB,8);
 208   _oggpack_write(opb,info->numbooks-1,4);
 209   for(j=0;j<info->numbooks;j++)
 210     _oggpack_write(opb,info->books[j],8);
 211 }
 212
 213 static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
 214   int j;
 215   vorbis_info_floor0 *info=malloc(sizeof(vorbis_info_floor0));
 216   info->order=_oggpack_read(opb,8);
 217   info->rate=_oggpack_read(opb,16);
 218   info->barkmap=_oggpack_read(opb,16);
 219   info->ampbits=_oggpack_read(opb,6);
 220   info->ampdB=_oggpack_read(opb,8);
 221   info->numbooks=_oggpack_read(opb,4)+1;
 222
 223   if(info->order<1)goto err_out;
 224   if(info->rate<1)goto err_out;
 225   if(info->barkmap<1)goto err_out;
 226   if(info->numbooks<1)goto err_out;
 227
 228   for(j=0;j<info->numbooks;j++){
 229     info->books[j]=_oggpack_read(opb,8);
 230     if(info->books[j]<0 || info->books[j]>=vi->books)goto err_out;
 231   }
 232   return(info);
 233  err_out:
 234   free_info(info);
 235   return(NULL);
 236 }
 237
 238 /* initialize Bark scale and normalization lookups.  We could do this
 239    with static tables, but Vorbis allows a number of possible
 240    combinations, so it's best to do it computationally.
 241
 242    The below is authoritative in terms of defining scale mapping.
 243    Note that the scale depends on the sampling rate as well as the
 244    linear block and mapping sizes */
 245
 246 static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 247                               vorbis_info_floor *i){
 248   int j;
 249   double scale;
 250   vorbis_info        *vi=vd->vi;
 251   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 252   vorbis_look_floor0 *look=malloc(sizeof(vorbis_look_floor0));
 253   look->m=info->order;
 254   look->n=vi->blocksizes[mi->blockflag]/2;
 255   look->ln=info->barkmap;
 256   look->vi=info;
 257   lpc_init(&look->lpclook,look->ln,look->m);
 258
 259   /* we choose a scaling constant so that:
 260      floor(bark(rate/2-1)*C)=mapped-1
 261      floor(bark(rate/2)*C)=mapped */
 262   scale=look->ln/toBARK(info->rate/2.);
 263
 264   /* the mapping from a linear scale to a smaller bark scale is
 265      straightforward.  We do *not* make sure that the linear mapping
 266      does not skip bark-scale bins; the decoder simply skips them and
 267      the encoder may do what it wishes in filling them.  They're
 268      necessary in some mapping combinations to keep the scale spacing
 269      accurate */
 270   look->linearmap=malloc(look->n*sizeof(int));
 271   for(j=0;j<look->n;j++){
 272     int val=floor( toBARK((info->rate/2.)/look->n*j)
 273                    *scale); /* bark numbers represent band edges */
 274     if(val>look->ln)val=look->ln; /* guard against the approximation */
 275     look->linearmap[j]=val;
 276   }
 277
 278   if(vd->analysisp){
 279     look->encodefit=calloc(info->numbooks,sizeof(LSP_fit_lookup));
 280     for(j=0;j<info->numbooks;j++){
 281       codebook *b=vd->fullbooks+info->books[j];
 282       _f0_fit_init(look->encodefit+j,b);
 283     }
 284   }
 285
 286   return look;
 287 }
 288
 289 /* less efficient than the decode side (written for clarity).  We're
 290    not bottlenecked here anyway */
 291
 292 double _curve_to_lpc(double *curve,double *lpc,
 293                      vorbis_look_floor0 *l,long frameno){
 294   /* map the input curve to a bark-scale curve for encoding */
 295
 296   int mapped=l->ln;
 297   double *work=alloca(sizeof(double)*mapped);
 298   int i,j,last=0;
 299   int bark=0;
 300
 301   memset(work,0,sizeof(double)*mapped);
 302
 303   /* Only the decode side is behavior-specced; for now in the encoder,
 304      we select the maximum value of each band as representative (this
 305      helps make sure peaks don't go out of range.  In error terms,
 306      selecting min would make more sense, but the codebook is trained
 307      numerically, so we don't actually lose.  We'd still want to
 308      use the original curve for error and noise estimation */
 309
 310   for(i=0;i<l->n;i++){
 311     bark=l->linearmap[i];
 312     if(work[bark]<curve[i])work[bark]=curve[i];
 313     if(bark>last+1){
 314       /* If the bark scale is climbing rapidly, some bins may end up
 315          going unused.  This isn't a waste actually; it keeps the
 316          scale resolution even so that the LPC generator has an easy
 317          time.  However, if we leave the bins empty we lose energy.
 318          So, fill 'em in.  The decoder does not do anything with  he
 319          unused bins, so we can fill them anyway we like to end up
 320          with a better spectral curve */
 321
 322       /* we'll always have a bin zero, so we don't need to guard init */
 323       long span=bark-last;
 324       for(j=1;j<span;j++){
 325         double del=(double)j/span;
 326         work[j+last]=work[bark]*del+work[last]*(1.-del);
 327       }
 328     }
 329     last=bark;
 330   }
 331
 332   /* If we're over-ranged to avoid edge effects, fill in the end of spectrum gap */
 333   for(i=bark+1;i<mapped;i++)
 334     work[i]=work[i-1];
 335
 336 #if 0
 337     { /******************/
 338       FILE *of;
 339       char buffer[80];
 340       int i;
 341
 342       sprintf(buffer,"Fmask_%d.m",frameno);
 343       of=fopen(buffer,"w");
 344       for(i=0;i<mapped;i++)
 345         fprintf(of,"%g\n",work[i]);
 346       fclose(of);
 347     }
 348 #endif
 349
 350   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 351 }
 352
 353 /* generate the whole freq response curve of an LPC IIR filter */
 354
 355 void _lpc_to_curve(double *curve,double *lpc,double amp,
 356                           vorbis_look_floor0 *l,char *name,long frameno){
 357   /* l->m+1 must be less than l->ln, but guard in case we get a bad stream */
 358   double *lcurve=alloca(sizeof(double)*max(l->ln*2,l->m*2+2));
 359   int i;
 360
 361   if(amp==0){
 362     memset(curve,0,sizeof(double)*l->n);
 363     return;
 364   }
 365   vorbis_lpc_to_curve(lcurve,lpc,amp,&(l->lpclook));
 366
 367 #if 0
 368     { /******************/
 369       FILE *of;
 370       char buffer[80];
 371       int i;
 372
 373       sprintf(buffer,"%s_%d.m",name,frameno);
 374       of=fopen(buffer,"w");
 375       for(i=0;i<l->ln;i++)
 376         fprintf(of,"%g\n",lcurve[i]);
 377       fclose(of);
 378     }
 379 #endif
 380
 381   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 382
 383 }
 384
 385 static long seq=0;
 386 static int forward(vorbis_block *vb,vorbis_look_floor *i,
 387                     double *in,double *out){
 388   long j;
 389   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 390   vorbis_info_floor0 *info=look->vi;
 391   double *work=alloca((look->ln+look->n)*sizeof(double));
 392   double amp;
 393   long bits=0;
 394
 395 #ifdef TRAIN_LSP
 396   FILE *of;
 397   FILE *ef;
 398   char buffer[80];
 399
 400 #if 0
 401   sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 402   of=fopen(buffer,"a");
 403 #endif
 404
 405   sprintf(buffer,"lsp0ent_%d.vqd",vb->mode);
 406   ef=fopen(buffer,"a");
 407 #endif
 408
 409   /* our floor comes in on a linear scale; go to a [-Inf...0] dB
 410      scale.  The curve has to be positive, so we offset it. */
 411
 412   for(j=0;j<look->n;j++){
 413     double val=todB(in[j])+info->ampdB;
 414     if(val<1.)
 415       work[j]=1.;
 416     else
 417       work[j]=val;
 418   }
 419
 420   /* use 'out' as temp storage */
 421   /* Convert our floor to a set of lpc coefficients */
 422   amp=sqrt(_curve_to_lpc(work,out,look,seq));
 423
 424   /* amp is in the range (0. to ampdB].  Encode that range using
 425      ampbits bits */
 426
 427   {
 428     long maxval=(1L<<info->ampbits)-1;
 429
 430     long val=rint(amp/info->ampdB*maxval);
 431
 432     if(val<0)val=0;           /* likely */
 433     if(val>maxval)val=maxval; /* not bloody likely */
 434
 435     _oggpack_write(&vb->opb,val,info->ampbits);
 436     if(val>0)
 437       amp=(float)val/maxval*info->ampdB;
 438     else
 439       amp=0;
 440   }
 441
 442   if(amp>0){
 443
 444     /* the spec supports using one of a number of codebooks.  Right
 445        now, encode using this lib supports only one */
 446     codebook *b=vb->vd->fullbooks+info->books[0];
 447     LSP_fit_lookup *f=look->encodefit;
 448     _oggpack_write(&vb->opb,0,_ilog(info->numbooks));
 449
 450     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 451     vorbis_lpc_to_lsp(out,out,look->m);
 452
 453 #ifdef ANALYSIS
 454     if(vb->W==0){fprintf(stderr,"%d ",seq);}
 455     vorbis_lsp_to_lpc(out,work,look->m);
 456     _lpc_to_curve(work,work,amp,look,"Ffloor",seq);
 457     for(j=0;j<look->n;j++)work[j]-=info->ampdB;
 458     _analysis_output("rawfloor",seq,work,look->n,0,0);
 459     {
 460       double last=0;
 461       for(j=0;j<look->m;j++){
 462         work[j]=out[j]-last;
 463         last=out[j];
 464       }
 465     }
 466     _analysis_output("rawlsp",seq,work,look->m,0,0);
 467
 468 #endif
 469
 470 #if 0
 471 #ifdef TRAIN_LSP
 472     for(j=0;j<look->m;j++)
 473       fprintf(of,"%.12g, ",out[j]);
 474     fprintf(of,"\n");
 475     fclose(of);
 476 #endif
 477 #endif
 478
 479     /* code the spectral envelope, and keep track of the actual
 480        quantized values; we don't want creeping error as each block is
 481        nailed to the last quantized value of the previous block. */
 482
 483     for(j=0;j<look->m;j+=b->dim){
 484       int entry=_f0_fit(f,b,out,work,j);
 485       bits+=vorbis_book_encode(b,entry,&vb->opb);
 486 #ifdef ANALYSIS
 487       if(entry==2921)fprintf(stderr,"\n*************found\n");
 488 #endif
 489
 490 #ifdef TRAIN_LSP
 491       fprintf(ef,"%d,\n",entry);
 492 #endif
 493
 494     }
 495
 496 #ifdef ANALYSIS
 497     {
 498       double last=0;
 499       for(j=0;j<look->m;j++){
 500         out[j]=work[j]-last;
 501         last=work[j];
 502       }
 503     }
 504     _analysis_output("lsp",seq,out,look->m,0,0);
 505
 506 #endif
 507
 508 #ifdef TRAIN_LSP
 509     fclose(ef);
 510 #endif
 511
 512     /* take the coefficients back to a spectral envelope curve */
 513     vorbis_lsp_to_lpc(work,out,look->m);
 514     _lpc_to_curve(out,out,amp,look,"Ffloor",seq++);
 515     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
 516     return(1);
 517   }
 518
 519   memset(out,0,sizeof(double)*look->n);
 520   seq++;
 521   return(0);
 522 }
 523
 524 static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
 525   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 526   vorbis_info_floor0 *info=look->vi;
 527   int j,k;
 528
 529   int ampraw=_oggpack_read(&vb->opb,info->ampbits);
 530   if(ampraw>0){ /* also handles the -1 out of data case */
 531     long maxval=(1<<info->ampbits)-1;
 532     double amp=(float)ampraw/maxval*info->ampdB;
 533     int booknum=_oggpack_read(&vb->opb,_ilog(info->numbooks));
 534
 535     if(booknum!=-1){
 536       codebook *b=vb->vd->fullbooks+info->books[booknum];
 537       double last=0.;
 538
 539       memset(out,0,sizeof(double)*look->m);
 540
 541       for(j=0;j<look->m;j+=b->dim)
 542         if(vorbis_book_decodevs(b,out+j,&vb->opb,1,-1)==-1)goto eop;
 543       for(j=0;j<look->m;){
 544         for(k=0;k<b->dim;k++,j++)out[j]+=last;
 545         last=out[j-1];
 546       }
 547
 548       /* take the coefficients back to a spectral envelope curve */
 549       vorbis_lsp_to_lpc(out,out,look->m);
 550       _lpc_to_curve(out,out,amp,look,"",0);
 551
 552       for(j=0;j<look->n;j++)out[j]=fromdB(out[j]-info->ampdB);
 553       return(1);
 554     }
 555   }
 556
 557  eop:
 558   memset(out,0,sizeof(double)*look->n);
 559   return(0);
 560 }
 561
 562 /* export hooks */
 563 vorbis_func_floor floor0_exportbundle={
 564   &pack,&unpack,&look,&free_info,&free_look,&forward,&inverse
 565 };
 566
 567