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