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