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.15 2000/06/14 01:38:31 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 typedef struct {
  34   long n;
  35   int ln;
  36   int  m;
  37   int *linearmap;
  38
  39   vorbis_info_floor0 *vi;
  40   lpc_lookup lpclook;
  41 } vorbis_look_floor0;
  42
  43 typedef struct {
  44   long *codewords;
  45   double *curve;
  46   long frameno;
  47   long codes;
  48 } vorbis_echstate_floor0;
  49
  50 static void free_info(vorbis_info_floor *i){
  51   if(i){
  52     memset(i,0,sizeof(vorbis_info_floor0));
  53     free(i);
  54   }
  55 }
  56
  57 static void free_look(vorbis_look_floor *i){
  58   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
  59   if(i){
  60     if(look->linearmap)free(look->linearmap);
  61     lpc_clear(&look->lpclook);
  62     memset(look,0,sizeof(vorbis_look_floor0));
  63     free(look);
  64   }
  65 }
  66
  67 static void pack (vorbis_info_floor *i,oggpack_buffer *opb){
  68   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
  69   int j;
  70   _oggpack_write(opb,info->order,8);
  71   _oggpack_write(opb,info->rate,16);
  72   _oggpack_write(opb,info->barkmap,16);
  73   _oggpack_write(opb,info->ampbits,6);
  74   _oggpack_write(opb,info->ampdB,8);
  75   _oggpack_write(opb,info->numbooks-1,4);
  76   for(j=0;j<info->numbooks;j++)
  77     _oggpack_write(opb,info->books[j],8);
  78 }
  79
  80 static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
  81   int j;
  82   vorbis_info_floor0 *info=malloc(sizeof(vorbis_info_floor0));
  83   info->order=_oggpack_read(opb,8);
  84   info->rate=_oggpack_read(opb,16);
  85   info->barkmap=_oggpack_read(opb,16);
  86   info->ampbits=_oggpack_read(opb,6);
  87   info->ampdB=_oggpack_read(opb,8);
  88   info->numbooks=_oggpack_read(opb,4)+1;
  89
  90   if(info->order<1)goto err_out;
  91   if(info->rate<1)goto err_out;
  92   if(info->barkmap<1)goto err_out;
  93   if(info->numbooks<1)goto err_out;
  94
  95   for(j=0;j<info->numbooks;j++){
  96     info->books[j]=_oggpack_read(opb,8);
  97     if(info->books[j]<0 || info->books[j]>=vi->books)goto err_out;
  98   }
  99   return(info);
 100  err_out:
 101   free_info(info);
 102   return(NULL);
 103 }
 104
 105 /* initialize Bark scale and normalization lookups.  We could do this
 106    with static tables, but Vorbis allows a number of possible
 107    combinations, so it's best to do it computationally.
 108
 109    The below is authoritative in terms of defining scale mapping.
 110    Note that the scale depends on the sampling rate as well as the
 111    linear block and mapping sizes */
 112
 113 static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 114                               vorbis_info_floor *i){
 115   int j;
 116   double scale;
 117   vorbis_info        *vi=vd->vi;
 118   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 119   vorbis_look_floor0 *look=malloc(sizeof(vorbis_look_floor0));
 120   look->m=info->order;
 121   look->n=vi->blocksizes[mi->blockflag]/2;
 122   look->ln=info->barkmap;
 123   look->vi=info;
 124   lpc_init(&look->lpclook,look->ln,look->m);
 125
 126   /* we choose a scaling constant so that:
 127      floor(bark(rate/2-1)*C)=mapped-1
 128      floor(bark(rate/2)*C)=mapped */
 129   scale=look->ln/toBARK(info->rate/2.);
 130
 131   /* the mapping from a linear scale to a smaller bark scale is
 132      straightforward.  We do *not* make sure that the linear mapping
 133      does not skip bark-scale bins; the decoder simply skips them and
 134      the encoder may do what it wishes in filling them.  They're
 135      necessary in some mapping combinations to keep the scale spacing
 136      accurate */
 137   look->linearmap=malloc(look->n*sizeof(int));
 138   for(j=0;j<look->n;j++){
 139     int val=floor( toBARK((info->rate/2.)/look->n*j)
 140                    *scale); /* bark numbers represent band edges */
 141     if(val>look->ln)val=look->ln; /* guard against the approximation */
 142     look->linearmap[j]=val;
 143   }
 144   return look;
 145 }
 146
 147 static vorbis_echstate_floor *state (vorbis_info_floor *i){
 148   vorbis_echstate_floor0 *state=calloc(1,sizeof(vorbis_echstate_floor0));
 149   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 150
 151   /* a safe size if usually too big (dim==1) */
 152   state->codewords=malloc(info->order*sizeof(long));
 153   state->curve=malloc(info->barkmap*sizeof(double));
 154   state->frameno=-1;
 155   return(state);
 156 }
 157
 158 static void free_state (vorbis_echstate_floor *vs){
 159   vorbis_echstate_floor0 *state=(vorbis_echstate_floor0 *)vs;
 160   if(state){
 161     free(state->codewords);
 162     free(state->curve);
 163     memset(state,0,sizeof(vorbis_echstate_floor0));
 164     free(state);
 165   }
 166 }
 167
 168 #include <stdio.h>
 169
 170 double _curve_error(double *curve1,double *curve2,long n){
 171   double acc=0.;
 172   long i;
 173   for(i=0;i<n;i++){
 174     double val=curve1[i]-curve2[i];
 175     acc+=val*val;
 176   }
 177   return(acc);
 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 _lpc_to_curve(double *curve,double *lpc,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)*max(l->ln*2,l->m*2+2));
 250   int i;
 251
 252   if(amp==0){
 253     memset(curve,0,sizeof(double)*l->n);
 254     return;
 255   }
 256   vorbis_lpc_to_curve(lcurve,lpc,amp,&(l->lpclook));
 257
 258 #if 0
 259     { /******************/
 260       FILE *of;
 261       char buffer[80];
 262       int i;
 263
 264       sprintf(buffer,"%s_%d.m",name,frameno);
 265       of=fopen(buffer,"w");
 266       for(i=0;i<l->ln;i++)
 267         fprintf(of,"%g\n",lcurve[i]);
 268       fclose(of);
 269     }
 270 #endif
 271
 272   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 273
 274 }
 275
 276 static long seq=0;
 277 static int forward(vorbis_block *vb,vorbis_look_floor *i,
 278                     double *in,double *out,vorbis_echstate_floor *vs){
 279   long j,k,l;
 280   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 281   vorbis_info_floor0 *info=look->vi;
 282   vorbis_echstate_floor0 *state=(vorbis_echstate_floor0 *)vs;
 283   double *work=alloca(look->n*sizeof(double));
 284   double amp;
 285   long bits=0;
 286
 287   /* our floor comes in on a linear scale; go to a [-Inf...0] dB
 288      scale.  The curve has to be positive, so we offset it. */
 289
 290   for(j=0;j<look->n;j++)work[j]=todB(in[j])+info->ampdB;
 291
 292   /* use 'out' as temp storage */
 293   /* Convert our floor to a set of lpc coefficients */
 294   amp=sqrt(_curve_to_lpc(work,out,look,seq));
 295
 296   /* amp is in the range (0. to ampdB].  Encode that range using
 297      ampbits bits */
 298
 299   {
 300     long maxval=(1L<<info->ampbits)-1;
 301
 302     long val=rint(amp/info->ampdB*maxval);
 303
 304     if(val<0)val=0;           /* likely */
 305     if(val>maxval)val=maxval; /* not bloody likely */
 306
 307     _oggpack_write(&vb->opb,val,info->ampbits);
 308     if(val>0)
 309       amp=(float)val/maxval*info->ampdB;
 310     else
 311       amp=0;
 312   }
 313
 314   if(amp>0){
 315     double *refcurve=alloca(sizeof(double)*max(look->ln*2,look->m*2+2));
 316     double *newcurve=alloca(sizeof(double)*max(look->ln*2,look->m*2+2));
 317     long *codelist=alloca(sizeof(long)*look->m);
 318     int codes=0;
 319
 320
 321     if(state->frameno == vb->sequence){
 322       /* generate a reference curve for testing */
 323       vorbis_lpc_to_curve(refcurve,out,1,&(look->lpclook));
 324     }
 325
 326     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 327     vorbis_lpc_to_lsp(out,out,look->m);
 328
 329 #ifdef TRAIN
 330     {
 331       int j;
 332       FILE *of;
 333       char buffer[80];
 334       sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 335       of=fopen(buffer,"a");
 336       for(j=0;j<look->m;j++)
 337         fprintf(of,"%.12g, ",out[j]);
 338       fprintf(of,"\n");
 339       fclose(of);
 340     }
 341 #endif
 342
 343     /* code the spectral envelope, and keep track of the actual
 344        quantized values; we don't want creeping error as each block is
 345        nailed to the last quantized value of the previous block. */
 346
 347     /* A new development: the code selection is based on error against
 348        the LSP coefficients and not the curve it produces.  Because
 349        the coefficient error is not linearly related to the curve
 350        error, the fit we select is usually nonoptimal (but
 351        sufficient).  This is fine, but flailing about causes a problem
 352        in generally consistent spectra... so we add hysterisis. */
 353
 354     /* select a new fit, but don't code it.  Just grab it for testing */
 355     {
 356       /* the spec supports using one of a number of codebooks.  Right
 357          now, encode using this lib supports only one */
 358       codebook *b=vb->vd->fullbooks+info->books[0];
 359       double last=0.;
 360
 361       for(j=0;j<look->m;){
 362         for(k=0;k<b->dim;k++)out[j+k]-=last;
 363         codelist[codes++]=vorbis_book_errorv(b,out+j);
 364         for(k=0;k<b->dim;k++,j++)out[j]+=last;
 365         last=out[j-1];
 366       }
 367     }
 368     vorbis_lsp_to_lpc(out,out,look->m);
 369     vorbis_lpc_to_curve(newcurve,out,1,&(look->lpclook));
 370
 371     /* if we're out of sequence, no hysterisis this frame, else check it */
 372     if(state->frameno != vb->sequence ||
 373        _curve_error(refcurve,newcurve,look->ln)<
 374        _curve_error(refcurve,state->curve,look->ln)){
 375       /* new curve is the fit to use.  replace the state */
 376       memcpy(state->curve,newcurve,sizeof(double)*look->ln);
 377       memcpy(state->codewords,codelist,sizeof(long)*codes);
 378       state->codes=codes;
 379     }else{
 380       /* use state */
 381       /*fprintf(stderr,"X");*/
 382       codelist=state->codewords;
 383       codes=state->codes;
 384     }
 385
 386     state->frameno=vb->sequence+1;
 387
 388     /* the spec supports using one of a number of codebooks.  Right
 389        now, encode using this lib supports only one */
 390     _oggpack_write(&vb->opb,0,_ilog(info->numbooks));
 391
 392     {
 393       codebook *b=vb->vd->fullbooks+info->books[0];
 394       double last=0.;
 395       for(l=0,j=0;l<codes;l++){
 396         bits+=vorbis_book_encodev(b,codelist[l],out+j,&vb->opb);
 397         for(k=0;k<b->dim;k++,j++)out[j]+=last;
 398         last=out[j-1];
 399       }
 400     }
 401
 402     /* take the coefficients back to a spectral envelope curve */
 403     vorbis_lsp_to_lpc(out,out,look->m);
 404     _lpc_to_curve(out,out,amp,look,"Ffloor",seq);
 405     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
 406     return(1);
 407   }
 408
 409   memset(out,0,sizeof(double)*look->n);
 410   seq++;
 411   return(0);
 412 }
 413
 414 static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
 415   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 416   vorbis_info_floor0 *info=look->vi;
 417   int j,k;
 418
 419   int ampraw=_oggpack_read(&vb->opb,info->ampbits);
 420   if(ampraw>0){
 421     long maxval=(1<<info->ampbits)-1;
 422     double amp=(float)ampraw/maxval*info->ampdB;
 423     int booknum=_oggpack_read(&vb->opb,_ilog(info->numbooks));
 424     codebook *b=vb->vd->fullbooks+info->books[booknum];
 425     double last=0.;
 426
 427     memset(out,0,sizeof(double)*look->m);
 428
 429     for(j=0;j<look->m;j+=b->dim)
 430       vorbis_book_decodevs(b,out+j,&vb->opb,1,-1);
 431     for(j=0;j<look->m;){
 432       for(k=0;k<b->dim;k++,j++)out[j]+=last;
 433       last=out[j-1];
 434     }
 435
 436     /* take the coefficients back to a spectral envelope curve */
 437     vorbis_lsp_to_lpc(out,out,look->m);
 438     _lpc_to_curve(out,out,amp,look,"",0);
 439
 440     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
 441     return(1);
 442   }else
 443     memset(out,0,sizeof(double)*look->n);
 444
 445   return(0);
 446 }
 447
 448 /* export hooks */
 449 vorbis_func_floor floor0_exportbundle={
 450   &pack,&unpack,&look,&state,&free_info,&free_look,&free_state,&forward,&inverse
 451 };
 452
 453