lib/floor0.c

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