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