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