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