lib/block.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-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: PCM data vector blocking, windowing and dis/reassembly
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Oct 2 1999
  18
  19  Handle windowing, overlap-add, etc of the PCM vectors.  This is made
  20  more amusing by Vorbis' current two allowed block sizes.
  21
  22  Vorbis manipulates the dynamic range of the incoming PCM data
  23  envelope to minimise time-domain energy leakage from percussive and
  24  plosive waveforms being quantized in the MDCT domain.
  25
  26  ********************************************************************/
  27
  28 #include <stdlib.h>
  29 #include <string.h>
  30 #include "codec.h"
  31 #include "window.h"
  32 #include "envelope.h"
  33 #include "mdct.h"
  34 #include "lpc.h"
  35
  36 /* pcm accumulator examples (not exhaustive):
  37
  38  <-------------- lW ---------------->
  39                    <--------------- W ---------------->
  40 :            .....|.....       _______________         |
  41 :        .'''     |     '''_---      |       |\        |
  42 :.....'''         |_____--- '''......|       | \_______|
  43 :.................|__________________|_______|__|______|
  44                   |<------ Sl ------>|      > Sr <     |endW
  45                   |beginSl           |endSl  |  |endSr
  46                   |beginW            |endlW  |beginSr
  47
  48
  49                       |< lW >|
  50                    <--------------- W ---------------->
  51                   |   |  ..  ______________            |
  52                   |   | '  `/        |     ---_        |
  53                   |___.'___/`.       |         ---_____|
  54                   |_______|__|_______|_________________|
  55                   |      >|Sl|<      |<------ Sr ----->|endW
  56                   |       |  |endSl  |beginSr          |endSr
  57                   |beginW |  |endlW
  58                   mult[0] |beginSl                     mult[n]
  59
  60  <-------------- lW ----------------->
  61                           |<--W-->|
  62 :            ..............  ___  |   |
  63 :        .'''             |`/   \ |   |
  64 :.....'''                 |/`....\|...|
  65 :.........................|___|___|___|
  66                           |Sl |Sr |endW
  67                           |   |   |endSr
  68                           |   |beginSr
  69                           |   |endSl
  70                           |beginSl
  71                           |beginW
  72 */
  73
  74 /* block abstraction setup *********************************************/
  75
  76 int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb){
  77   int i;
  78   vorbis_info *vi=v->vi;
  79   memset(vb,0,sizeof(vorbis_block));
  80   vb->vd=v;
  81
  82   vb->pcm_storage=vi->blocksize[1];
  83   vb->pcm_channels=vi->channels;
  84   vb->mult_storage=vi->blocksize[1]/vi->envelopesa;
  85   vb->mult_channels=vi->envelopech;
  86   vb->floor_channels=vi->floorch;
  87   vb->floor_storage=max(vi->floororder[0],vi->floororder[1]);
  88
  89   vb->pcm=malloc(vb->pcm_channels*sizeof(double *));
  90   for(i=0;i<vb->pcm_channels;i++)
  91     vb->pcm[i]=malloc(vb->pcm_storage*sizeof(double));
  92
  93   vb->mult=malloc(vb->mult_channels*sizeof(double *));
  94   for(i=0;i<vb->mult_channels;i++)
  95     vb->mult[i]=malloc(vb->mult_storage*sizeof(double));
  96
  97   vb->lsp=malloc(vb->floor_channels*sizeof(double *));
  98   vb->lpc=malloc(vb->floor_channels*sizeof(double *));
  99   vb->amp=malloc(vb->floor_channels*sizeof(double));
 100   for(i=0;i<vb->floor_channels;i++){
 101     vb->lsp[i]=malloc(vb->floor_storage*sizeof(double));
 102     vb->lpc[i]=malloc(vb->floor_storage*sizeof(double));
 103   }
 104   if(v->analysisp)
 105     _oggpack_writeinit(&vb->opb);
 106
 107   return(0);
 108 }
 109
 110 int vorbis_block_clear(vorbis_block *vb){
 111   int i;
 112   if(vb->pcm){
 113     for(i=0;i<vb->pcm_channels;i++)
 114       free(vb->pcm[i]);
 115     free(vb->pcm);
 116   }
 117   if(vb->mult){
 118     for(i=0;i<vb->mult_channels;i++)
 119       free(vb->mult[i]);
 120     free(vb->mult);
 121   }
 122   if(vb->vd->analysisp)
 123     _oggpack_writeclear(&vb->opb);
 124
 125   memset(vb,0,sizeof(vorbis_block));
 126   return(0);
 127 }
 128
 129 /* Analysis side code, but directly related to blocking.  Thus it's
 130    here and not in analysis.c (which is for analysis transforms only).
 131    The init is here because some of it is shared */
 132
 133 static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
 134   memset(v,0,sizeof(vorbis_dsp_state));
 135
 136   v->vi=vi;
 137   _ve_envelope_init(&v->ve,vi->envelopesa);
 138   mdct_init(&v->vm[0],vi->blocksize[0]);
 139   mdct_init(&v->vm[1],vi->blocksize[1]);
 140
 141   v->window[0][0][0]=_vorbis_window(vi->blocksize[0],
 142                                    vi->blocksize[0]/2,vi->blocksize[0]/2);
 143   v->window[0][0][1]=v->window[0][0][0];
 144   v->window[0][1][0]=v->window[0][0][0];
 145   v->window[0][1][1]=v->window[0][0][0];
 146
 147   v->window[1][0][0]=_vorbis_window(vi->blocksize[1],
 148                                    vi->blocksize[0]/2,vi->blocksize[0]/2);
 149   v->window[1][0][1]=_vorbis_window(vi->blocksize[1],
 150                                    vi->blocksize[0]/2,vi->blocksize[1]/2);
 151   v->window[1][1][0]=_vorbis_window(vi->blocksize[1],
 152                                    vi->blocksize[1]/2,vi->blocksize[0]/2);
 153   v->window[1][1][1]=_vorbis_window(vi->blocksize[1],
 154                                     vi->blocksize[1]/2,vi->blocksize[1]/2);
 155
 156   /* initialize the storage vectors to a decent size greater than the
 157      minimum */
 158
 159   v->pcm_storage=8192; /* we'll assume later that we have
 160                           a minimum of twice the blocksize of
 161                           accumulated samples in analysis */
 162   v->pcm=malloc(vi->channels*sizeof(double *));
 163   v->pcmret=malloc(vi->channels*sizeof(double *));
 164   {
 165     int i;
 166     for(i=0;i<vi->channels;i++)
 167       v->pcm[i]=calloc(v->pcm_storage,sizeof(double));
 168   }
 169
 170   /* Initialize the envelope multiplier storage */
 171
 172   if(vi->envelopech){
 173     v->envelope_storage=v->pcm_storage/vi->envelopesa;
 174     v->multipliers=calloc(vi->envelopech,sizeof(double *));
 175     {
 176       int i;
 177       for(i=0;i<vi->envelopech;i++){
 178         v->multipliers[i]=calloc(v->envelope_storage,sizeof(double));
 179       }
 180     }
 181   }
 182
 183   /* all 1 (large block) or 0 (small block) */
 184   /* explicitly set for the sake of clarity */
 185   v->lW=0; /* previous window size */
 186   v->W=0;  /* current window size */
 187
 188   /* all vector indexes; multiples of samples_per_envelope_step */
 189   v->centerW=vi->blocksize[1]/2;
 190
 191   v->pcm_current=v->centerW;
 192   v->envelope_current=v->centerW/vi->envelopesa;
 193   return(0);
 194 }
 195
 196 /* arbitrary settings and spec-mandated numbers get filled in here */
 197 int vorbis_analysis_init(vorbis_dsp_state *v,vorbis_info *vi){
 198   _vds_shared_init(v,vi);
 199
 200   /* the coder init is different for read/write */
 201   v->analysisp=1;
 202
 203   /* Yes, wasteful to have four lookups.  This will get collapsed once
 204      things crystallize */
 205   lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
 206            vi->floororder[0],vi->flooroctaves[0],1);
 207   lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/2,
 208            vi->floororder[0],vi->flooroctaves[0],1);
 209
 210   /*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
 211            vi->balanceorder,vi->balanceoctaves,1);
 212   lpc_init(&v->vbal[1],vi->blocksize[1]/2,256,
 213            vi->balanceorder,vi->balanceoctaves,1);*/
 214
 215   return(0);
 216 }
 217
 218 void vorbis_dsp_clear(vorbis_dsp_state *v){
 219   int i,j,k;
 220   if(v){
 221     vorbis_info *vi=v->vi;
 222
 223     if(v->window[0][0][0])free(v->window[0][0][0]);
 224     for(j=0;j<2;j++)
 225       for(k=0;k<2;k++)
 226         if(v->window[1][j][k])free(v->window[1][j][k]);
 227     if(v->pcm){
 228       for(i=0;i<vi->channels;i++)
 229         if(v->pcm[i])free(v->pcm[i]);
 230       free(v->pcm);
 231       if(v->pcmret)free(v->pcmret);
 232     }
 233     if(v->multipliers){
 234       for(i=0;i<vi->envelopech;i++)
 235         if(v->multipliers[i])free(v->multipliers[i]);
 236       free(v->multipliers);
 237     }
 238     _ve_envelope_clear(&v->ve);
 239     mdct_clear(&v->vm[0]);
 240     mdct_clear(&v->vm[1]);
 241     lpc_clear(&v->vl[0]);
 242     lpc_clear(&v->vl[1]);
 243     /*lpc_clear(&v->vbal[0]);
 244       lpc_clear(&v->vbal[1]);*/
 245     memset(v,0,sizeof(vorbis_dsp_state));
 246   }
 247 }
 248
 249 double **vorbis_analysis_buffer(vorbis_dsp_state *v, int vals){
 250   int i;
 251   vorbis_info *vi=v->vi;
 252
 253   /* Do we have enough storage space for the requested buffer? If not,
 254      expand the PCM (and envelope) storage */
 255
 256   if(v->pcm_current+vals>=v->pcm_storage){
 257     v->pcm_storage=v->pcm_current+vals*2;
 258     v->envelope_storage=v->pcm_storage/v->vi->envelopesa;
 259
 260     for(i=0;i<vi->channels;i++){
 261       v->pcm[i]=realloc(v->pcm[i],v->pcm_storage*sizeof(double));
 262     }
 263     for(i=0;i<vi->envelopech;i++){
 264       v->multipliers[i]=realloc(v->multipliers[i],
 265                                 v->envelope_storage*sizeof(double));
 266     }
 267   }
 268
 269   for(i=0;i<vi->channels;i++)
 270     v->pcmret[i]=v->pcm[i]+v->pcm_current;
 271
 272   return(v->pcmret);
 273 }
 274
 275 /* call with val<=0 to set eof */
 276
 277 int vorbis_analysis_wrote(vorbis_dsp_state *v, int vals){
 278   vorbis_info *vi=v->vi;
 279   if(vals<=0){
 280     /* We're encoding the end of the stream.  Just make sure we have
 281        [at least] a full block of zeroes at the end. */
 282
 283     int i;
 284     vorbis_analysis_buffer(v,v->vi->blocksize[1]*2);
 285     v->eofflag=v->pcm_current;
 286     v->pcm_current+=v->vi->blocksize[1]*2;
 287     for(i=0;i<vi->channels;i++)
 288       memset(v->pcm[i]+v->eofflag,0,
 289              (v->pcm_current-v->eofflag)*sizeof(double));
 290   }else{
 291
 292     if(v->pcm_current+vals>v->pcm_storage)
 293       return(-1);
 294
 295     v->pcm_current+=vals;
 296   }
 297   return(0);
 298 }
 299
 300 /* do the deltas, envelope shaping, pre-echo and determine the size of
 301    the next block on which to continue analysis */
 302 int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
 303   int i,j;
 304   vorbis_info *vi=v->vi;
 305   long beginW=v->centerW-vi->blocksize[v->W]/2,centerNext;
 306   long beginM=beginW/vi->envelopesa;
 307
 308   /* check to see if we're done... */
 309   if(v->eofflag==-1)return(0);
 310
 311   /* if we have any unfilled envelope blocks for which we have PCM
 312      data, fill them up in before proceeding. */
 313
 314   if(v->pcm_current/vi->envelopesa>v->envelope_current){
 315     /* This generates the multipliers, but does not sparsify the vector.
 316        That's done by block before coding */
 317     _ve_envelope_multipliers(v);
 318   }
 319
 320   /* By our invariant, we have lW, W and centerW set.  Search for
 321      the next boundary so we can determine nW (the next window size)
 322      which lets us compute the shape of the current block's window */
 323
 324   /* overconserve for now; any block with a non-placeholder multiplier
 325      should be minimal size. We can be greedy and only look at nW size */
 326
 327   if(vi->blocksize[0]<vi->blocksize[1]){
 328
 329     if(v->W)
 330       /* this is a long window; we start the search forward of centerW
 331          because that's the fastest we could react anyway */
 332       i=v->centerW+vi->blocksize[1]/4-vi->blocksize[0]/4;
 333     else
 334       /* short window.  Search from centerW */
 335       i=v->centerW;
 336     i/=vi->envelopesa;
 337
 338     for(;i<v->envelope_current;i++){
 339       for(j=0;j<vi->envelopech;j++)
 340         if(v->multipliers[j][i-1]*vi->preecho_thresh<
 341            v->multipliers[j][i])break;
 342       if(j<vi->envelopech)break;
 343     }
 344
 345     if(i<v->envelope_current){
 346       /* Ooo, we hit a multiplier. Is it beyond the boundary to make the
 347          upcoming block large ? */
 348       long largebound;
 349       if(v->W)
 350         largebound=v->centerW+vi->blocksize[1];
 351       else
 352         largebound=v->centerW+vi->blocksize[0]/4+vi->blocksize[1]*3/4;
 353       largebound/=vi->envelopesa;
 354
 355       if(i>=largebound)
 356         v->nW=1;
 357       else
 358         v->nW=0;
 359
 360     }else{
 361       /* Assume maximum; if the block is incomplete given current
 362          buffered data, this will be detected below */
 363       v->nW=1;
 364     }
 365   }else
 366     v->nW=1;
 367     v->nW=1;
 368
 369   /* Do we actually have enough data *now* for the next block? The
 370      reason to check is that if we had no multipliers, that could
 371      simply been due to running out of data.  In that case, we don;t
 372      know the size of the next block for sure and we need that now to
 373      figure out the window shape of this block */
 374
 375   centerNext=v->centerW+vi->blocksize[v->W]/4+vi->blocksize[v->nW]/4;
 376
 377   {
 378     long blockbound=centerNext+vi->blocksize[v->nW]/2;
 379     if(v->pcm_current<blockbound)return(0); /* not enough data yet */
 380   }
 381
 382   /* fill in the block */
 383   vb->lW=v->lW;
 384   vb->W=v->W;
 385   vb->nW=v->nW;
 386   vb->vd=v;
 387
 388   vb->pcmend=vi->blocksize[v->W];
 389   vb->multend=vb->pcmend / vi->envelopesa;
 390
 391   /* copy the vectors */
 392   for(i=0;i<vi->channels;i++)
 393     memcpy(vb->pcm[i],v->pcm[i]+beginW,vi->blocksize[v->W]*sizeof(double));
 394   for(i=0;i<vi->envelopech;i++)
 395     memcpy(vb->mult[i],v->multipliers[i]+beginM,vi->blocksize[v->W]/
 396            vi->envelopesa*sizeof(double));
 397
 398   vb->frameno=v->frame;
 399
 400   /* handle eof detection: eof==0 means that we've not yet received EOF
 401                            eof>0  marks the last 'real' sample in pcm[]
 402                            eof<0  'no more to do'; doesn't get here */
 403
 404   if(v->eofflag){
 405     if(v->centerW>=v->eofflag){
 406       v->eofflag=-1;
 407       vb->eofflag=1;
 408     }
 409   }
 410
 411   /* advance storage vectors and clean up */
 412   {
 413     int new_centerNext=vi->blocksize[1]/2;
 414     int movementW=centerNext-new_centerNext;
 415     int movementM=movementW/vi->envelopesa;
 416
 417     /* the multipliers and pcm stay synced up because the blocksizes
 418        must be multiples of samples_per_envelope_step (minimum
 419        multiple is 2) */
 420
 421     for(i=0;i<vi->channels;i++)
 422       memmove(v->pcm[i],v->pcm[i]+movementW,
 423               (v->pcm_current-movementW)*sizeof(double));
 424
 425     for(i=0;i<vi->envelopech;i++){
 426       memmove(v->multipliers[i],v->multipliers[i]+movementM,
 427               (v->envelope_current-movementM)*sizeof(double));
 428     }
 429
 430     v->pcm_current-=movementW;
 431     v->envelope_current-=movementM;
 432
 433     v->lW=v->W;
 434     v->W=v->nW;
 435     v->centerW=new_centerNext;
 436
 437     v->frame++;
 438     v->samples+=movementW;
 439
 440     if(v->eofflag)
 441       v->eofflag-=movementW;
 442   }
 443
 444   /* done */
 445   return(1);
 446 }
 447
 448 int vorbis_synthesis_init(vorbis_dsp_state *v,vorbis_info *vi){
 449   int temp=vi->envelopech;
 450   vi->envelopech=0; /* we don't need multiplier buffering in syn */
 451   _vds_shared_init(v,vi);
 452   vi->envelopech=temp;
 453
 454   /* Yes, wasteful to have four lookups.  This will get collapsed once
 455      things crystallize */
 456   lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
 457            vi->floororder[0],vi->flooroctaves[0],0);
 458   lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/2,
 459            vi->floororder[1],vi->flooroctaves[1],0);
 460   /*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
 461            vi->balanceorder,vi->balanceoctaves,0);
 462   lpc_init(&v->vbal[1],vi->blocksize[1]/2,256,
 463            vi->balanceorder,vi->balanceoctaves,0);*/
 464
 465
 466   /* Adjust centerW to allow an easier mechanism for determining output */
 467   v->pcm_returned=v->centerW;
 468   v->centerW-= vi->blocksize[v->W]/4+vi->blocksize[v->lW]/4;
 469   return(0);
 470 }
 471
 472 /* Unike in analysis, the window is only partially applied.  Envelope
 473    is previously applied (the whole envelope, if any, is shipped in
 474    each block) */
 475
 476 int vorbis_synthesis_blockin(vorbis_dsp_state *v,vorbis_block *vb){
 477   vorbis_info *vi=v->vi;
 478
 479   /* Shift out any PCM that we returned previously */
 480
 481   if(v->pcm_returned  && v->centerW>vi->blocksize[1]/2){
 482
 483     /* don't shift too much; we need to have a minimum PCM buffer of
 484        1/2 long block */
 485
 486     int shift=v->centerW-vi->blocksize[1]/2;
 487     shift=(v->pcm_returned<shift?v->pcm_returned:shift);
 488
 489     v->pcm_current-=shift;
 490     v->centerW-=shift;
 491     v->pcm_returned-=shift;
 492
 493     if(shift){
 494       int i;
 495       for(i=0;i<vi->channels;i++)
 496         memmove(v->pcm[i],v->pcm[i]+shift,
 497                 v->pcm_current*sizeof(double));
 498     }
 499   }
 500
 501   v->lW=v->W;
 502   v->W=vb->W;
 503
 504   {
 505     int sizeW=vi->blocksize[v->W];
 506     int centerW=v->centerW+vi->blocksize[v->lW]/4+sizeW/4;
 507     int beginW=centerW-sizeW/2;
 508     int endW=beginW+sizeW;
 509     int beginSl;
 510     int endSl;
 511     int i,j;
 512     double *window;
 513
 514     /* Do we have enough PCM storage for the block? */
 515     if(endW>v->pcm_storage){
 516       /* expand the PCM storage */
 517
 518       v->pcm_storage=endW+vi->blocksize[1];
 519
 520       for(i=0;i<vi->channels;i++)
 521         v->pcm[i]=realloc(v->pcm[i],v->pcm_storage*sizeof(double));
 522     }
 523
 524     /* Overlap/add */
 525     switch(v->W){
 526     case 0:
 527       beginSl=0;
 528       endSl=vi->blocksize[0]/2;
 529       break;
 530     case 1:
 531       beginSl=vi->blocksize[1]/4-vi->blocksize[v->lW]/4;
 532       endSl=beginSl+vi->blocksize[v->lW]/2;
 533       break;
 534     }
 535
 536     window=v->window[v->W][0][v->lW]+vi->blocksize[v->W]/2;
 537
 538     for(j=0;j<vi->channels;j++){
 539       double *pcm=v->pcm[j]+beginW;
 540
 541       /* the add section */
 542       for(i=beginSl;i<endSl;i++)
 543         pcm[i]=pcm[i]*window[i]+vb->pcm[j][i];
 544       /* the remaining section */
 545       for(;i<sizeW;i++)
 546         pcm[i]=vb->pcm[j][i];
 547     }
 548
 549     /* Update, cleanup */
 550
 551     v->centerW=centerW;
 552     v->pcm_current=endW;
 553
 554     if(vb->eofflag)v->eofflag=1;
 555   }
 556   return(0);
 557 }
 558
 559 int vorbis_synthesis_pcmout(vorbis_dsp_state *v,double ***pcm){
 560   vorbis_info *vi=v->vi;
 561   if(v->pcm_returned<v->centerW){
 562     int i;
 563     for(i=0;i<vi->channels;i++)
 564       v->pcmret[i]=v->pcm[i]+v->pcm_returned;
 565     *pcm=v->pcmret;
 566     return(v->centerW-v->pcm_returned);
 567   }
 568   return(0);
 569 }
 570
 571 int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes){
 572   if(bytes && v->pcm_returned+bytes>v->centerW)return(-1);
 573   v->pcm_returned+=bytes;
 574   return(0);
 575 }
 576