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.28 2000/11/14 00:05:31 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     _ogg_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)_ogg_free(look->linearmap);
  91     if(look->lsp_look)_ogg_free(look->lsp_look);
  92     lpc_clear(&look->lpclook);
  93     memset(look,0,sizeof(vorbis_look_floor0));
  94     _ogg_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   static int seq=0;
 250
 251 #ifdef TRAIN_LSP
 252   FILE *of;
 253   FILE *ef;
 254   char buffer[80];
 255
 256 #if 1
 257   sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 258   of=fopen(buffer,"a");
 259 #endif
 260
 261   sprintf(buffer,"lsp0ent_%d.vqd",vb->mode);
 262   ef=fopen(buffer,"a");
 263 #endif
 264
 265   /* our floor comes in on a linear scale; go to a [-Inf...0] dB
 266      scale.  The curve has to be positive, so we offset it. */
 267
 268   for(j=0;j<look->n;j++)
 269     work[j]=todB(in[j])+info->ampdB;
 270
 271   /* use 'out' as temp storage */
 272   /* Convert our floor to a set of lpc coefficients */
 273   amp=sqrt(_curve_to_lpc(work,out,look));
 274
 275   /* amp is in the range (0. to ampdB].  Encode that range using
 276      ampbits bits */
 277
 278   {
 279     long maxval=(1L<<info->ampbits)-1;
 280
 281     long val=rint(amp/info->ampdB*maxval);
 282
 283     if(val<0)val=0;           /* likely */
 284     if(val>maxval)val=maxval; /* not bloody likely */
 285
 286     oggpack_write(&vb->opb,val,info->ampbits);
 287     if(val>0)
 288       amp=(float)val/maxval*info->ampdB;
 289     else
 290       amp=0;
 291   }
 292
 293   if(amp>0){
 294
 295     /* the spec supports using one of a number of codebooks.  Right
 296        now, encode using this lib supports only one */
 297     backend_lookup_state *be=vb->vd->backend_state;
 298     codebook *b=be->fullbooks+info->books[0];
 299     oggpack_write(&vb->opb,0,_ilog(info->numbooks));
 300
 301     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 302     vorbis_lpc_to_lsp(out,out,look->m);
 303
 304 #ifdef ANALYSIS
 305     {
 306       float *lspwork=alloca(look->m*sizeof(float));
 307       memcpy(lspwork,out,look->m*sizeof(float));
 308       vorbis_lsp_to_curve(lspwork,look->linearmap,look->n,look->ln,
 309                           work,look->m,amp,info->ampdB);
 310       _analysis_output("prefit",seq,work,look->n,0,1);
 311
 312     }
 313
 314 #endif
 315
 316
 317 #if 1
 318 #ifdef TRAIN_LSP
 319     {
 320       float last=0.;
 321       for(j=0;j<look->m;j++){
 322         fprintf(of,"%.12g, ",out[j]-last);
 323         last=out[j];
 324       }
 325     }
 326     fprintf(of,"\n");
 327     fclose(of);
 328 #endif
 329 #endif
 330
 331     /* code the spectral envelope, and keep track of the actual
 332        quantized values; we don't want creeping error as each block is
 333        nailed to the last quantized value of the previous block. */
 334
 335     for(j=0;j<look->m;j+=b->dim){
 336       int entry=_f0_fit(b,out,work,j);
 337       bits+=vorbis_book_encode(b,entry,&vb->opb);
 338
 339 #ifdef TRAIN_LSP
 340       fprintf(ef,"%d,\n",entry);
 341 #endif
 342
 343     }
 344
 345 #ifdef ANALYSIS
 346     {
 347       float last=0;
 348       for(j=0;j<look->m;j++){
 349         out[j]=work[j]-last;
 350         last=work[j];
 351       }
 352     }
 353
 354 #endif
 355
 356 #ifdef TRAIN_LSP
 357     fclose(ef);
 358 #endif
 359
 360     /* take the coefficients back to a spectral envelope curve */
 361     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 362                         work,look->m,amp,info->ampdB);
 363     return(1);
 364   }
 365
 366   memset(out,0,sizeof(float)*look->n);
 367   seq++;
 368   return(0);
 369 }
 370
 371 static int floor0_inverse(vorbis_block *vb,vorbis_look_floor *i,float *out){
 372   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 373   vorbis_info_floor0 *info=look->vi;
 374   int j,k;
 375
 376   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 377   if(ampraw>0){ /* also handles the -1 out of data case */
 378     long maxval=(1<<info->ampbits)-1;
 379     float amp=(float)ampraw/maxval*info->ampdB;
 380     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 381     float *lsp=alloca(sizeof(float)*look->m);
 382
 383     if(booknum!=-1){
 384       backend_lookup_state *be=vb->vd->backend_state;
 385       codebook *b=be->fullbooks+info->books[booknum];
 386       float last=0.;
 387
 388       memset(out,0,sizeof(float)*look->m);
 389
 390       for(j=0;j<look->m;j+=b->dim)
 391         if(vorbis_book_decodevs(b,lsp+j,&vb->opb,1,-1)==-1)goto eop;
 392       for(j=0;j<look->m;){
 393         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 394         last=lsp[j-1];
 395       }
 396
 397       /* take the coefficients back to a spectral envelope curve */
 398       vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 399                           lsp,look->m,amp,info->ampdB);
 400       return(1);
 401     }
 402   }
 403
 404  eop:
 405   memset(out,0,sizeof(float)*look->n);
 406   return(0);
 407 }
 408
 409 /* export hooks */
 410 vorbis_func_floor floor0_exportbundle={
 411   &floor0_pack,&floor0_unpack,&floor0_look,&floor0_copy_info,&floor0_free_info,
 412   &floor0_free_look,&floor0_forward,&floor0_inverse
 413 };
 414
 415