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