lib/envelope.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 envelope analysis and manipulation
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Oct 17 1999
  18
  19  Vorbis manipulates the dynamic range of the incoming PCM data
  20  envelope to minimise time-domain energy leakage from percussive and
  21  plosive waveforms being quantized in the MDCT domain.
  22
  23  ********************************************************************/
  24
  25 #include <stdlib.h>
  26 #include <string.h>
  27 #include <stdio.h>
  28 #include <math.h>
  29
  30 #include "codec.h"
  31 #include "mdct.h"
  32 #include "envelope.h"
  33 #include "bitwise.h"
  34
  35 void _ve_envelope_init(envelope_lookup *e,int samples_per){
  36   int i;
  37
  38   e->winlen=samples_per*2;
  39   e->window=malloc(e->winlen*sizeof(double));
  40
  41   /* We just use a straight sin^2(x) window for this */
  42   for(i=0;i<e->winlen;i++){
  43     double temp=sin((i+.5)/e->winlen*M_PI);
  44     e->window[i]=temp*temp;
  45   }
  46 }
  47
  48 void _ve_envelope_clear(envelope_lookup *e){
  49   if(e->window)free(e->window);
  50   memset(e,0,sizeof(envelope_lookup));
  51 }
  52
  53 /* initial and final blocks are special cases. Eg:
  54    ______
  55          `--_
  56    |_______|_`-.___|_______|_______|
  57
  58               ___
  59           _--'   `--_
  60    |___.-'_|_______|_`-.___|_______|
  61
  62                       ___
  63                   _--'   `--_
  64    |_______|___.-'_|_______|_`-.___|
  65
  66                                _____
  67                            _--'
  68    |_______|_______|____.-'|_______|
  69
  70    as we go block by block, we watch the collective metrics span. If we
  71    span the threshhold (assuming the threshhold is active), we use an
  72    abbreviated vector */
  73
  74 static int _ve_envelope_generate(double *mult,double *env,double *look,
  75                                  int n,int step){
  76   int i,j,p;
  77   double m,mo;
  78
  79   for(j=0;j<n;j++)if(mult[j]!=1)break;
  80   if(j==n)return(0);
  81
  82   n*=step;
  83   /* first multiplier special case */
  84   m=ldexp(1,mult[0]-1);
  85   for(p=0;p<step/2;p++)env[p]=m;
  86
  87   /* mid multipliers normal case */
  88   for(j=1;p<n-step/2;j++){
  89     mo=m;
  90     m=ldexp(1,mult[j]-1);
  91     if(mo==m)
  92       for(i=0;i<step;i++,p++)env[p]=m;
  93     else
  94       for(i=0;i<step;i++,p++)env[p]=m*look[i]+mo*look[i+step];
  95   }
  96
  97   /* last multiplier special case */
  98   for(;p<n;p++)env[p]=m;
  99   return(1);
 100 }
 101
 102 /* use MDCT for spectral power estimation */
 103
 104 static void _ve_deltas(double *deltas,double *pcm,int n,double *win,
 105                        int winsize){
 106   int i,j;
 107   mdct_lookup m;
 108   double out[winsize/2];
 109
 110   mdct_init(&m,winsize);
 111
 112   for(j=0;j<n;j++){
 113     double acc=0.;
 114
 115     mdct_forward(&m,pcm+j*winsize/2,out,win);
 116     for(i=0;i<winsize/2;i++)
 117       acc+=fabs(out[i]);
 118     if(deltas[j]<acc)deltas[j]=acc;
 119   }
 120
 121   mdct_clear(&m);
 122
 123 #ifdef ANALYSIS
 124   {
 125     static int frameno=0;
 126     FILE *out;
 127     char buffer[80];
 128
 129     sprintf(buffer,"delta%d.m",frameno++);
 130     out=fopen(buffer,"w+");
 131     for(j=0;j<n;j++)
 132       fprintf(out,"%g\n",deltas[j]);
 133     fclose(out);
 134   }
 135 #endif
 136 }
 137
 138 void _ve_envelope_multipliers(vorbis_dsp_state *v){
 139   vorbis_info *vi=v->vi;
 140   int step=vi->envelopesa;
 141   static int frame=0;
 142
 143   /* we need a 1-1/4 envelope window overlap begin and 1/4 end */
 144   int dtotal=(v->pcm_current-step/2)/vi->envelopesa;
 145   int dcurr=v->envelope_current;
 146   double *window=v->ve.window;
 147   int winlen=v->ve.winlen;
 148   int pch,ech;
 149
 150   if(dtotal>dcurr){
 151     for(ech=0;ech<vi->envelopech;ech++){
 152       double *mult=v->multipliers[ech]+dcurr;
 153       memset(mult,0,sizeof(double)*(dtotal-dcurr));
 154
 155       for(pch=0;pch<vi->channels;pch++){
 156
 157         /* does this channel contribute to the envelope analysis */
 158         /*if(vi->envelopemap[pch]==ech){ not mapping yet */
 159         if(pch==ech){
 160
 161           /* we need a 1/4 envelope window overlap front and back */
 162           double *pcm=v->pcm[pch]+dcurr*step-step/2;
 163           _ve_deltas(mult,pcm,dtotal-dcurr,window,winlen);
 164
 165         }
 166       }
 167     }
 168     v->envelope_current=dtotal;
 169     frame++;
 170   }
 171 }
 172
 173 /* This readies the multiplier vector for use/coding.  Clamp/adjust
 174    the multipliers to the allowed range and eliminate unneeded
 175    coefficients */
 176
 177 void _ve_envelope_sparsify(vorbis_block *vb){
 178   vorbis_info *vi=vb->vd->vi;
 179   int ch;
 180   for(ch=0;ch<vi->envelopech;ch++){
 181     int flag=0;
 182     double *mult=vb->mult[ch];
 183     int n=vb->multend;
 184     double first=mult[0];
 185     double last=first;
 186     double clamp;
 187     int i;
 188
 189 #ifdef ANALYSIS
 190     {
 191       static int frameno=0.;
 192       FILE *out;
 193       int j;
 194       char buffer[80];
 195
 196       sprintf(buffer,"del%d.m",frameno++);
 197       out=fopen(buffer,"w+");
 198       for(j=0;j<n;j++)
 199         fprintf(out,"%g\n",mult[j]);
 200       fclose(out);
 201     }
 202 #endif
 203
 204     /* are we going to multiply anything? */
 205
 206     for(i=1;i<n;i++){
 207       if(mult[i]>=last*vi->preecho_thresh){
 208         flag=1;
 209         break;
 210       }
 211       if(i<n-1 && mult[i+1]>=last*vi->preecho_thresh){
 212         flag=1;
 213         break;
 214       }
 215       last=mult[i];
 216     }
 217
 218     if(flag){
 219       /* we need to adjust, so we might as well go nuts */
 220
 221       int begin=i;
 222       clamp=last?last:1;
 223
 224       for(i=0;i<begin;i++)mult[i]=0;
 225
 226       for(;i<n;i++){
 227         if(mult[i]>=last*vi->preecho_thresh){
 228           last=mult[i];
 229
 230           mult[i]=floor(log(mult[i]/clamp)/log(2));
 231           if(mult[i]>15)mult[i]=15;
 232           if(mult[i]<1)mult[i]=1;
 233
 234         }else{
 235           mult[i]=0;
 236         }
 237       }
 238     }else
 239       memset(mult,0,sizeof(double)*n);
 240
 241 #ifdef ANALYSIS
 242     {
 243       static int frameno=0.;
 244       FILE *out;
 245       int j;
 246       char buffer[80];
 247
 248       sprintf(buffer,"sparse%d.m",frameno++);
 249       out=fopen(buffer,"w+");
 250       for(j=0;j<n;j++)
 251         fprintf(out,"%g\n",mult[j]);
 252       fclose(out);
 253     }
 254 #endif
 255
 256
 257   }
 258 }
 259
 260 void _ve_envelope_apply(vorbis_block *vb,int multp){
 261   static int frameno=0;
 262   vorbis_info *vi=vb->vd->vi;
 263   double env[vb->multend*vi->envelopesa];
 264   envelope_lookup *look=&vb->vd->ve;
 265   int i,j,k;
 266
 267   for(i=0;i<vi->envelopech;i++){
 268     double *mult=vb->mult[i];
 269     double last=1.;
 270
 271     /* fill in the multiplier placeholders */
 272
 273     for(j=0;j<vb->multend;j++){
 274       if(mult[j]){
 275         last=mult[j];
 276       }else{
 277         mult[j]=last;
 278       }
 279     }
 280
 281     /* compute the envelope curve */
 282     if(_ve_envelope_generate(mult,env,look->window,vb->multend,
 283                              vi->envelopesa)){
 284 #ifdef ANALYSIS
 285       {
 286         FILE *out;
 287         char buffer[80];
 288
 289         sprintf(buffer,"env%d.m",frameno);
 290         out=fopen(buffer,"w+");
 291         for(j=0;j<vb->multend*vi->envelopesa;j++)
 292           fprintf(out,"%g\n",env[j]);
 293         fclose(out);
 294         sprintf(buffer,"prepcm%d.m",frameno);
 295         out=fopen(buffer,"w+");
 296         for(j=0;j<vb->multend*vi->envelopesa;j++)
 297           fprintf(out,"%g\n",vb->pcm[i][j]);
 298         fclose(out);
 299       }
 300 #endif
 301
 302       for(k=0;k<vi->channels;k++){
 303
 304         if(multp)
 305           for(j=0;j<vb->multend*vi->envelopesa;j++)
 306             vb->pcm[k][j]*=env[j];
 307         else
 308           for(j=0;j<vb->multend*vi->envelopesa;j++)
 309             vb->pcm[k][j]/=env[j];
 310
 311       }
 312
 313 #ifdef ANALYSIS
 314       {
 315         FILE *out;
 316         char buffer[80];
 317
 318         sprintf(buffer,"pcm%d.m",frameno);
 319         out=fopen(buffer,"w+");
 320         for(j=0;j<vb->multend*vi->envelopesa;j++)
 321           fprintf(out,"%g\n",vb->pcm[i][j]);
 322         fclose(out);
 323       }
 324 #endif
 325     }
 326     frameno++;
 327   }
 328 }
 329
 330 int _ve_envelope_encode(vorbis_block *vb){
 331   /* Not huffcoded yet. */
 332
 333   vorbis_info *vi=vb->vd->vi;
 334   int scale=vb->W;
 335   int n=vb->multend;
 336   int i,j;
 337
 338   /* range is 0-15 */
 339
 340   for(i=0;i<vi->envelopech;i++){
 341     double *mult=vb->mult[i];
 342     for(j=0;j<n;j++){
 343       _oggpack_write(&vb->opb,(int)(mult[j]),4);
 344     }
 345   }
 346   return(0);
 347 }
 348
 349 /* synthesis expands the buffers in vb if needed.  We can assume we
 350    have enough storage handed to us. */
 351 int _ve_envelope_decode(vorbis_block *vb){
 352   vorbis_info *vi=vb->vd->vi;
 353   int scale=vb->W;
 354   int n=vb->multend;
 355   int i,j;
 356
 357   /* range is 0-15 */
 358
 359   for(i=0;i<vi->envelopech;i++){
 360     double *mult=vb->mult[i];
 361     for(j=0;j<n;j++)
 362       mult[j]=_oggpack_read(&vb->opb,4);
 363   }
 364   return(0);
 365 }
 366
 367
 368
 369
 370