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.25 2000/10/12 03:12:52 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 "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   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 void floor0_free_info(vorbis_info_floor *i){
  74   if(i){
  75     memset(i,0,sizeof(vorbis_info_floor0));
  76     free(i);
  77   }
  78 }
  79
  80 static void floor0_free_look(vorbis_look_floor *i){
  81   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
  82   if(i){
  83     if(look->linearmap)free(look->linearmap);
  84     if(look->lsp_look)free(look->lsp_look);
  85     lpc_clear(&look->lpclook);
  86     memset(look,0,sizeof(vorbis_look_floor0));
  87     free(look);
  88   }
  89 }
  90
  91 static void floor0_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 *floor0_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   floor0_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 *floor0_look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 138                               vorbis_info_floor *i){
 139   int j;
 140   float scale;
 141   vorbis_info        *vi=vd->vi;
 142   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 143   vorbis_look_floor0 *look=calloc(1,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
 149   if(vd->analysisp)
 150     lpc_init(&look->lpclook,look->ln,look->m);
 151
 152   /* we choose a scaling constant so that:
 153      floor(bark(rate/2-1)*C)=mapped-1
 154      floor(bark(rate/2)*C)=mapped */
 155   scale=look->ln/toBARK(info->rate/2.);
 156
 157   /* the mapping from a linear scale to a smaller bark scale is
 158      straightforward.  We do *not* make sure that the linear mapping
 159      does not skip bark-scale bins; the decoder simply skips them and
 160      the encoder may do what it wishes in filling them.  They're
 161      necessary in some mapping combinations to keep the scale spacing
 162      accurate */
 163   look->linearmap=malloc((look->n+1)*sizeof(int));
 164   for(j=0;j<look->n;j++){
 165     int val=floor( toBARK((info->rate/2.)/look->n*j)
 166                    *scale); /* bark numbers represent band edges */
 167     if(val>look->ln)val=look->ln; /* guard against the approximation */
 168     look->linearmap[j]=val;
 169   }
 170   look->linearmap[j]=-1;
 171
 172   look->lsp_look=malloc(look->ln*sizeof(float));
 173   for(j=0;j<look->ln;j++)
 174     look->lsp_look[j]=2*cos(M_PI/look->ln*j);
 175
 176   return look;
 177 }
 178
 179 /* less efficient than the decode side (written for clarity).  We're
 180    not bottlenecked here anyway */
 181
 182 float _curve_to_lpc(float *curve,float *lpc,
 183                      vorbis_look_floor0 *l){
 184   /* map the input curve to a bark-scale curve for encoding */
 185
 186   int mapped=l->ln;
 187   float *work=alloca(sizeof(float)*mapped);
 188   int i,j,last=0;
 189   int bark=0;
 190
 191   memset(work,0,sizeof(float)*mapped);
 192
 193   /* Only the decode side is behavior-specced; for now in the encoder,
 194      we select the maximum value of each band as representative (this
 195      helps make sure peaks don't go out of range.  In error terms,
 196      selecting min would make more sense, but the codebook is trained
 197      numerically, so we don't actually lose.  We'd still want to
 198      use the original curve for error and noise estimation */
 199
 200   for(i=0;i<l->n;i++){
 201     bark=l->linearmap[i];
 202     if(work[bark]<curve[i])work[bark]=curve[i];
 203     if(bark>last+1){
 204       /* If the bark scale is climbing rapidly, some bins may end up
 205          going unused.  This isn't a waste actually; it keeps the
 206          scale resolution even so that the LPC generator has an easy
 207          time.  However, if we leave the bins empty we lose energy.
 208          So, fill 'em in.  The decoder does not do anything with  he
 209          unused bins, so we can fill them anyway we like to end up
 210          with a better spectral curve */
 211
 212       /* we'll always have a bin zero, so we don't need to guard init */
 213       long span=bark-last;
 214       for(j=1;j<span;j++){
 215         float del=(float)j/span;
 216         work[j+last]=work[bark]*del+work[last]*(1.-del);
 217       }
 218     }
 219     last=bark;
 220   }
 221
 222   /* If we're over-ranged to avoid edge effects, fill in the end of spectrum gap */
 223   for(i=bark+1;i<mapped;i++)
 224     work[i]=work[i-1];
 225
 226   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 227 }
 228
 229 /* generate the whole freq response curve of an LSP IIR filter */
 230 static int floor0_forward(vorbis_block *vb,vorbis_look_floor *i,
 231                     float *in,float *out){
 232   long j;
 233   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 234   vorbis_info_floor0 *info=look->vi;
 235   float *work=alloca((look->ln+look->n)*sizeof(float));
 236   float amp;
 237   long bits=0;
 238
 239 #ifdef TRAIN_LSP
 240   FILE *of;
 241   FILE *ef;
 242   char buffer[80];
 243
 244 #if 1
 245   sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 246   of=fopen(buffer,"a");
 247 #endif
 248
 249   sprintf(buffer,"lsp0ent_%d.vqd",vb->mode);
 250   ef=fopen(buffer,"a");
 251 #endif
 252
 253   /* our floor comes in on a linear scale; go to a [-Inf...0] dB
 254      scale.  The curve has to be positive, so we offset it. */
 255
 256   for(j=0;j<look->n;j++)
 257     work[j]=todB(in[j])+info->ampdB;
 258
 259   /* use 'out' as temp storage */
 260   /* Convert our floor to a set of lpc coefficients */
 261   amp=sqrt(_curve_to_lpc(work,out,look));
 262
 263   /* amp is in the range (0. to ampdB].  Encode that range using
 264      ampbits bits */
 265
 266   {
 267     long maxval=(1L<<info->ampbits)-1;
 268
 269     long val=rint(amp/info->ampdB*maxval);
 270
 271     if(val<0)val=0;           /* likely */
 272     if(val>maxval)val=maxval; /* not bloody likely */
 273
 274     oggpack_write(&vb->opb,val,info->ampbits);
 275     if(val>0)
 276       amp=(float)val/maxval*info->ampdB;
 277     else
 278       amp=0;
 279   }
 280
 281   if(amp>0){
 282
 283     /* the spec supports using one of a number of codebooks.  Right
 284        now, encode using this lib supports only one */
 285     codebook *b=vb->vd->fullbooks+info->books[0];
 286     oggpack_write(&vb->opb,0,_ilog(info->numbooks));
 287
 288     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 289     vorbis_lpc_to_lsp(out,out,look->m);
 290
 291 #if 1
 292 #ifdef TRAIN_LSP
 293     {
 294       float last=0.;
 295       for(j=0;j<look->m;j++){
 296         fprintf(of,"%.12g, ",out[j]-last);
 297         last=out[j];
 298       }
 299     }
 300     fprintf(of,"\n");
 301     fclose(of);
 302 #endif
 303 #endif
 304
 305     /* code the spectral envelope, and keep track of the actual
 306        quantized values; we don't want creeping error as each block is
 307        nailed to the last quantized value of the previous block. */
 308
 309     for(j=0;j<look->m;j+=b->dim){
 310       int entry=_f0_fit(b,out,work,j);
 311       bits+=vorbis_book_encode(b,entry,&vb->opb);
 312
 313 #ifdef TRAIN_LSP
 314       fprintf(ef,"%d,\n",entry);
 315 #endif
 316
 317     }
 318
 319 #ifdef ANALYSIS
 320     {
 321       float last=0;
 322       for(j=0;j<look->m;j++){
 323         out[j]=work[j]-last;
 324         last=work[j];
 325       }
 326     }
 327
 328 #endif
 329
 330 #ifdef TRAIN_LSP
 331     fclose(ef);
 332 #endif
 333
 334     /* take the coefficients back to a spectral envelope curve */
 335     vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 336                         work,look->m,amp,info->ampdB);
 337     return(1);
 338   }
 339
 340   memset(out,0,sizeof(float)*look->n);
 341   return(0);
 342 }
 343
 344 static int floor0_inverse(vorbis_block *vb,vorbis_look_floor *i,float *out){
 345   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 346   vorbis_info_floor0 *info=look->vi;
 347   int j,k;
 348
 349   int ampraw=oggpack_read(&vb->opb,info->ampbits);
 350   if(ampraw>0){ /* also handles the -1 out of data case */
 351     long maxval=(1<<info->ampbits)-1;
 352     float amp=(float)ampraw/maxval*info->ampdB;
 353     int booknum=oggpack_read(&vb->opb,_ilog(info->numbooks));
 354     float *lsp=alloca(sizeof(float)*look->m);
 355
 356     if(booknum!=-1){
 357       codebook *b=vb->vd->fullbooks+info->books[booknum];
 358       float last=0.;
 359
 360       memset(out,0,sizeof(double)*look->m);
 361
 362       for(j=0;j<look->m;j+=b->dim)
 363         if(vorbis_book_decodevs(b,lsp+j,&vb->opb,1,-1)==-1)goto eop;
 364       for(j=0;j<look->m;){
 365         for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
 366         last=lsp[j-1];
 367       }
 368
 369       /* take the coefficients back to a spectral envelope curve */
 370       vorbis_lsp_to_curve(out,look->linearmap,look->n,look->ln,
 371                           lsp,look->m,amp,info->ampdB);
 372       return(1);
 373     }
 374   }
 375
 376  eop:
 377   memset(out,0,sizeof(float)*look->n);
 378   return(0);
 379 }
 380
 381 /* export hooks */
 382 vorbis_func_floor floor0_exportbundle={
 383   &floor0_pack,&floor0_unpack,&floor0_look,&floor0_free_info,
 384   &floor0_free_look,&floor0_forward,&floor0_inverse
 385 };
 386
 387