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.14 2000/05/08 20:49:48 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 static void free_info(vorbis_info_floor *i){
  44   if(i){
  45     memset(i,0,sizeof(vorbis_info_floor0));
  46     free(i);
  47   }
  48 }
  49
  50 static void free_look(vorbis_look_floor *i){
  51   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
  52   if(i){
  53     if(look->linearmap)free(look->linearmap);
  54     lpc_clear(&look->lpclook);
  55     memset(look,0,sizeof(vorbis_look_floor0));
  56     free(look);
  57   }
  58 }
  59
  60 static void pack (vorbis_info_floor *i,oggpack_buffer *opb){
  61   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
  62   int j;
  63   _oggpack_write(opb,info->order,8);
  64   _oggpack_write(opb,info->rate,16);
  65   _oggpack_write(opb,info->barkmap,16);
  66   _oggpack_write(opb,info->ampbits,6);
  67   _oggpack_write(opb,info->ampdB,8);
  68   _oggpack_write(opb,info->numbooks-1,4);
  69   for(j=0;j<info->numbooks;j++)
  70     _oggpack_write(opb,info->books[j],8);
  71 }
  72
  73 static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
  74   int j;
  75   vorbis_info_floor0 *info=malloc(sizeof(vorbis_info_floor0));
  76   info->order=_oggpack_read(opb,8);
  77   info->rate=_oggpack_read(opb,16);
  78   info->barkmap=_oggpack_read(opb,16);
  79   info->ampbits=_oggpack_read(opb,6);
  80   info->ampdB=_oggpack_read(opb,8);
  81   info->numbooks=_oggpack_read(opb,4)+1;
  82
  83   if(info->order<1)goto err_out;
  84   if(info->rate<1)goto err_out;
  85   if(info->barkmap<1)goto err_out;
  86   if(info->numbooks<1)goto err_out;
  87
  88   for(j=0;j<info->numbooks;j++){
  89     info->books[j]=_oggpack_read(opb,8);
  90     if(info->books[j]<0 || info->books[j]>=vi->books)goto err_out;
  91   }
  92   return(info);
  93  err_out:
  94   free_info(info);
  95   return(NULL);
  96 }
  97
  98 /* initialize Bark scale and normalization lookups.  We could do this
  99    with static tables, but Vorbis allows a number of possible
 100    combinations, so it's best to do it computationally.
 101
 102    The below is authoritative in terms of defining scale mapping.
 103    Note that the scale depends on the sampling rate as well as the
 104    linear block and mapping sizes */
 105
 106 static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
 107                               vorbis_info_floor *i){
 108   int j;
 109   double scale;
 110   vorbis_info        *vi=vd->vi;
 111   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
 112   vorbis_look_floor0 *look=malloc(sizeof(vorbis_look_floor0));
 113   look->m=info->order;
 114   look->n=vi->blocksizes[mi->blockflag]/2;
 115   look->ln=info->barkmap;
 116   look->vi=info;
 117   lpc_init(&look->lpclook,look->ln,look->m);
 118
 119   /* we choose a scaling constant so that:
 120      floor(bark(rate/2-1)*C)=mapped-1
 121      floor(bark(rate/2)*C)=mapped */
 122   scale=look->ln/toBARK(info->rate/2.);
 123
 124   /* the mapping from a linear scale to a smaller bark scale is
 125      straightforward.  We do *not* make sure that the linear mapping
 126      does not skip bark-scale bins; the decoder simply skips them and
 127      the encoder may do what it wishes in filling them.  They're
 128      necessary in some mapping combinations to keep the scale spacing
 129      accurate */
 130   look->linearmap=malloc(look->n*sizeof(int));
 131   for(j=0;j<look->n;j++){
 132     int val=floor( toBARK((info->rate/2.)/look->n*j)
 133                    *scale); /* bark numbers represent band edges */
 134     if(val>look->ln)val=look->ln; /* guard against the approximation */
 135     look->linearmap[j]=val;
 136   }
 137   return look;
 138 }
 139
 140 #include <stdio.h>
 141
 142
 143 /* less efficient than the decode side (written for clarity).  We're
 144    not bottlenecked here anyway */
 145 double _curve_to_lpc(double *curve,double *lpc,vorbis_look_floor0 *l,
 146                             long frameno){
 147   /* map the input curve to a bark-scale curve for encoding */
 148
 149   int mapped=l->ln;
 150   double *work=alloca(sizeof(double)*mapped);
 151   int i,j,last=0;
 152
 153   memset(work,0,sizeof(double)*mapped);
 154
 155   /* Only the decode side is behavior-specced; for now in the encoder,
 156      we select the maximum value of each band as representative (this
 157      helps make sure peaks don't go out of range.  In error terms,
 158      selecting min would make more sense, but the codebook is trained
 159      numerically, so we don't actually lose.  We'd still want to
 160      use the original curve for error and noise estimation */
 161
 162   for(i=0;i<l->n;i++){
 163     int bark=l->linearmap[i];
 164     if(work[bark]<curve[i])work[bark]=curve[i];
 165     if(bark>last+1){
 166       /* If the bark scale is climbing rapidly, some bins may end up
 167          going unused.  This isn't a waste actually; it keeps the
 168          scale resolution even so that the LPC generator has an easy
 169          time.  However, if we leave the bins empty we lose energy.
 170          So, fill 'em in.  The decoder does not do anything with  he
 171          unused bins, so we can fill them anyway we like to end up
 172          with a better spectral curve */
 173
 174       /* we'll always have a bin zero, so we don't need to guard init */
 175       long span=bark-last;
 176       for(j=1;j<span;j++){
 177         double del=(double)j/span;
 178         work[j+last]=work[bark]*del+work[last]*(1.-del);
 179       }
 180     }
 181     last=bark;
 182   }
 183
 184 #if 0
 185     { /******************/
 186       FILE *of;
 187       char buffer[80];
 188       int i;
 189
 190       sprintf(buffer,"Fmask_%d.m",frameno);
 191       of=fopen(buffer,"w");
 192       for(i=0;i<mapped;i++)
 193         fprintf(of,"%g\n",work[i]);
 194       fclose(of);
 195     }
 196 #endif
 197
 198   return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
 199 }
 200
 201 /* generate the whole freq response curve of an LPC IIR filter */
 202
 203 void _lpc_to_curve(double *curve,double *lpc,double amp,
 204                           vorbis_look_floor0 *l,char *name,long frameno){
 205   double *lcurve=alloca(sizeof(double)*(l->ln*2));
 206   int i;
 207
 208   if(amp==0){
 209     memset(curve,0,sizeof(double)*l->n);
 210     return;
 211   }
 212   vorbis_lpc_to_curve(lcurve,lpc,amp,&(l->lpclook));
 213
 214 #if 0
 215     { /******************/
 216       FILE *of;
 217       char buffer[80];
 218       int i;
 219
 220       sprintf(buffer,"%s_%d.m",name,frameno);
 221       of=fopen(buffer,"w");
 222       for(i=0;i<l->ln;i++)
 223         fprintf(of,"%g\n",lcurve[i]);
 224       fclose(of);
 225     }
 226 #endif
 227
 228   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 229
 230 }
 231
 232 static long seq=0;
 233 static int forward(vorbis_block *vb,vorbis_look_floor *i,
 234                     double *in,double *out){
 235   long j,k;
 236   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 237   vorbis_info_floor0 *info=look->vi;
 238   double *work=alloca(look->n*sizeof(double));
 239   double amp;
 240   long bits=0;
 241
 242   /* our floor comes in on a linear scale; go to a [-Inf...0] dB
 243      scale.  The curve has to be positive, so we offset it. */
 244   for(j=0;j<look->n;j++)work[j]=todB(in[j])+info->ampdB;
 245
 246   /* use 'out' as temp storage */
 247   /* Convert our floor to a set of lpc coefficients */
 248   amp=sqrt(_curve_to_lpc(work,out,look,vb->sequence));
 249
 250   /* amp is in the range (0. to ampdB].  Encode that range using
 251      ampbits bits */
 252
 253   {
 254     long maxval=(1<<info->ampbits)-1;
 255
 256     long val=rint(amp/info->ampdB*maxval);
 257
 258     if(val<0)val=0;           /* likely */
 259     if(val>maxval)val=maxval; /* not bloody likely */
 260
 261     _oggpack_write(&vb->opb,val,info->ampbits);
 262     if(val>0)
 263       amp=(float)val/maxval*info->ampdB;
 264     else
 265       amp=0;
 266   }
 267
 268   if(amp>0){
 269
 270     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
 271     vorbis_lpc_to_lsp(out,out,look->m);
 272 #ifdef ANALYSIS
 273     if(vb->mode==0)_analysis_output("lsp",seq++,out,look->m,0,0);
 274 #endif
 275 #ifdef TRAIN
 276     {
 277       int j;
 278       FILE *of;
 279       char buffer[80];
 280       sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
 281       of=fopen(buffer,"a");
 282       for(j=0;j<look->m;j++)
 283         fprintf(of,"%g, ",out[j]);
 284       fprintf(of,"\n");
 285       fclose(of);
 286     }
 287 #endif
 288
 289 #if 0
 290     { /******************/
 291       vorbis_lsp_to_lpc(out,work,look->m);
 292       _lpc_to_curve(work,work,amp,look,"Fprefloor",vb->sequence);
 293     }
 294 #endif
 295
 296     /* code the spectral envelope, and keep track of the actual
 297        quantized values; we don't want creeping error as each block is
 298        nailed to the last quantized value of the previous block. */
 299
 300     /* the spec supports using one of a number of codebooks.  Right
 301        now, encode using this lib supports only one */
 302     _oggpack_write(&vb->opb,0,_ilog(info->numbooks));
 303
 304     {
 305       codebook *b=vb->vd->fullbooks+info->books[0];
 306       double last=0.;
 307       for(j=0;j<look->m;){
 308         for(k=0;k<b->dim;k++)out[j+k]-=last;
 309         bits+=vorbis_book_encodev(b,out+j,&vb->opb);
 310         for(k=0;k<b->dim;k++,j++)out[j]+=last;
 311         last=out[j-1];
 312       }
 313     }
 314
 315     /* take the coefficients back to a spectral envelope curve */
 316     vorbis_lsp_to_lpc(out,out,look->m);
 317     _lpc_to_curve(out,out,amp,look,"Ffloor",vb->sequence);
 318     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
 319     return(1);
 320   }
 321
 322   memset(out,0,sizeof(double)*look->n);
 323   return(0);
 324 }
 325
 326 static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
 327   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
 328   vorbis_info_floor0 *info=look->vi;
 329   int j,k;
 330
 331   int ampraw=_oggpack_read(&vb->opb,info->ampbits);
 332   if(ampraw>0){
 333     long maxval=(1<<info->ampbits)-1;
 334     double amp=(float)ampraw/maxval*info->ampdB;
 335     int booknum=_oggpack_read(&vb->opb,_ilog(info->numbooks));
 336     codebook *b=vb->vd->fullbooks+info->books[booknum];
 337     double last=0.;
 338
 339     memset(out,0,sizeof(double)*look->m);
 340
 341     for(j=0;j<look->m;j+=b->dim)
 342       vorbis_book_decodevs(b,out+j,&vb->opb,1,-1);
 343     for(j=0;j<look->m;){
 344       for(k=0;k<b->dim;k++,j++)out[j]+=last;
 345       last=out[j-1];
 346     }
 347
 348     /* take the coefficients back to a spectral envelope curve */
 349     vorbis_lsp_to_lpc(out,out,look->m);
 350     _lpc_to_curve(out,out,amp,look,"",0);
 351
 352     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
 353     return(1);
 354   }else
 355     memset(out,0,sizeof(double)*look->n);
 356
 357   return(0);
 358 }
 359
 360 /* export hooks */
 361 vorbis_func_floor floor0_exportbundle={
 362   &pack,&unpack,&look,&free_info,&free_look,&forward,&inverse
 363 };
 364
 365