lib/envelope.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH    *
   6  * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.        *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: PCM data envelope analysis and manipulation
  14  last mod: $Id: envelope.c,v 1.32 2001/02/15 19:05:45 xiphmont Exp $
  15
  16  Preecho calculation.
  17
  18  ********************************************************************/
  19
  20 #include <stdlib.h>
  21 #include <string.h>
  22 #include <stdio.h>
  23 #include <math.h>
  24 #include <ogg/ogg.h>
  25 #include "vorbis/codec.h"
  26 #include "codec_internal.h"
  27
  28 #include "os.h"
  29 #include "scales.h"
  30 #include "envelope.h"
  31 #include "misc.h"
  32
  33 /* Digital filter designed by mkfilter/mkshape/gencode A.J. Fisher */
  34
  35
  36
  37 static int   cheb_highpass_stages=6;
  38 static float cheb_highpass_B[]={1.f,-6.f,15.f,-20.f,15.f,-6.f,1.f};
  39
  40 static int   cheb_bandpass_stages=6;
  41 static float cheb_bandpass_B[]={-1.f,0.f,3.f,0.f,-3.f,0.f,1.f};
  42
  43
  44 /* 10kHz Chebyshev highpass */
  45 static float cheb_highpass10k_gain= 54.34519586f;
  46 static float cheb_highpass10k_A[]={
  47   -0.2064797169f,
  48   -0.5609713214f,
  49   -1.1352465327f,
  50   -1.4495555418f,
  51   -1.7938140760f,
  52   -0.9473564683f};
  53
  54 /* 6kHz-10kHz Chebyshev bandpass */
  55 static float cheb_bandpass6k_gain=113.4643935f;
  56 static float cheb_bandpass6k_A[]={
  57    -0.5712621337f,
  58     1.5626130710f,
  59    -3.3348854983f,
  60     4.0471340821f,
  61    -4.0051680331f,
  62    2.2786325610f};
  63
  64 /* 3kHz-6kHz Chebyshev bandpass */
  65 static float cheb_bandpass3k_gain= 248.8359377f;
  66 static float cheb_bandpass3k_A[]={
  67      -0.6564230022f,
  68       3.3747911257f,
  69      -8.0098635981f,
  70      11.0040876874f,
  71      -9.2250963484f,
  72       4.4760355389f};
  73
  74 /* 1.5kHz-3kHz Chebyshev bandpass */
  75 static float cheb_bandpass1k_gain= 1798.537183f;
  76 static float cheb_bandpass1k_A[]={
  77      -0.8097527363f,
  78       4.7725742682f,
  79     -11.9800219408f,
  80      16.3770336223f,
  81     -12.8553129536f,
  82      5.4948074309f};
  83
  84
  85 void _ve_envelope_init(envelope_lookup *e,vorbis_info *vi){
  86   codec_setup_info *ci=vi->codec_setup;
  87   int ch=vi->channels;
  88   int window=ci->envelopesa;
  89   int i;
  90   e->winlength=window;
  91   e->minenergy=fromdB(ci->preecho_minenergy);
  92   e->iir=_ogg_calloc(ch*4,sizeof(IIR_state));
  93   e->filtered=_ogg_calloc(ch*4,sizeof(float *));
  94   e->ch=ch;
  95   e->storage=128;
  96   for(i=0;i<ch*4;i+=4){
  97
  98     IIR_init(e->iir+i,cheb_highpass_stages,cheb_highpass10k_gain,
  99              cheb_highpass10k_A,cheb_highpass_B);
 100     IIR_init(e->iir+i+1,cheb_bandpass_stages,cheb_bandpass6k_gain,
 101              cheb_bandpass6k_A,cheb_bandpass_B);
 102     IIR_init(e->iir+i+2,cheb_bandpass_stages,cheb_bandpass3k_gain,
 103              cheb_bandpass3k_A,cheb_bandpass_B);
 104     IIR_init(e->iir+i+3,cheb_bandpass_stages,cheb_bandpass1k_gain,
 105              cheb_bandpass1k_A,cheb_bandpass_B);
 106
 107     e->filtered[i]=_ogg_calloc(e->storage,sizeof(float));
 108     e->filtered[i+1]=_ogg_calloc(e->storage,sizeof(float));
 109     e->filtered[i+2]=_ogg_calloc(e->storage,sizeof(float));
 110     e->filtered[i+3]=_ogg_calloc(e->storage,sizeof(float));
 111   }
 112
 113 }
 114
 115 void _ve_envelope_clear(envelope_lookup *e){
 116   int i;
 117   for(i=0;i<e->ch*4;i++){
 118     IIR_clear((e->iir+i));
 119     _ogg_free(e->filtered[i]);
 120   }
 121   _ogg_free(e->filtered);
 122   _ogg_free(e->iir);
 123   memset(e,0,sizeof(envelope_lookup));
 124 }
 125
 126 /* straight threshhold based until we find something that works better
 127    and isn't patented */
 128 static float _ve_deltai(envelope_lookup *ve,float *pre,float *post){
 129   long n=ve->winlength;
 130
 131   long i;
 132
 133   /* we want to have a 'minimum bar' for energy, else we're just
 134      basing blocks on quantization noise that outweighs the signal
 135      itself (for low power signals) */
 136
 137   float min=ve->minenergy;
 138   float A=min*min*n;
 139   float B=A;
 140
 141   /*_analysis_output("A",granulepos,pre,n,0,0);
 142     _analysis_output("B",granulepos,post,n,0,0);*/
 143
 144   for(i=0;i<n;i++){
 145     A+=pre[i]*pre[i];
 146     B+=post[i]*post[i];
 147   }
 148
 149   A=todB(A);
 150   B=todB(B);
 151
 152   return(B-A);
 153 }
 154
 155 long _ve_envelope_search(vorbis_dsp_state *v,long searchpoint){
 156   vorbis_info *vi=v->vi;
 157   codec_setup_info *ci=vi->codec_setup;
 158   envelope_lookup *ve=((backend_lookup_state *)(v->backend_state))->ve;
 159   long i,j,k,l;
 160   float *work=alloca(sizeof(float)*ve->winlength*2);
 161   static int seq=0;
 162
 163   /* make sure we have enough storage to match the PCM */
 164   if(v->pcm_storage>ve->storage){
 165     ve->storage=v->pcm_storage;
 166     for(i=0;i<ve->ch*4;i++)
 167       ve->filtered[i]=_ogg_realloc(ve->filtered[i],ve->storage*sizeof(float));
 168   }
 169
 170   /* catch up the highpass to match the pcm */
 171   for(i=0;i<ve->ch;i++){
 172     float *pcm=v->pcm[i];
 173     float *filtered0=ve->filtered[i*4];
 174     float *filtered1=ve->filtered[i*4+1];
 175     float *filtered2=ve->filtered[i*4+2];
 176     float *filtered3=ve->filtered[i*4+3];
 177     IIR_state *iir0=ve->iir+i*4;
 178     IIR_state *iir1=ve->iir+i*4+1;
 179     IIR_state *iir2=ve->iir+i*4+2;
 180     IIR_state *iir3=ve->iir+i*4+3;
 181     int flag=1;
 182
 183     for(j=ve->current;j<v->pcm_current;j++){
 184       filtered0[j]=IIR_filter(iir0,pcm[j]);
 185       filtered1[j]=IIR_filter(iir1,pcm[j]);
 186       filtered2[j]=IIR_filter(iir2,pcm[j]);
 187       filtered3[j]=IIR_filter(iir3,pcm[j]);
 188       if(pcm[j])flag=0;
 189     }
 190     if(flag && ve->current+64<v->pcm_current){
 191       IIR_reset(iir0);
 192       IIR_reset(iir1);
 193       IIR_reset(iir2);
 194       IIR_reset(iir3);
 195     }
 196
 197     _analysis_output("pcm",seq,pcm+v->centerW,v->pcm_current-v->centerW,0,0);
 198     _analysis_output("f0",seq,filtered0+v->centerW,v->pcm_current-v->centerW,
 199                      0,0);
 200     _analysis_output("f1",seq,filtered1+v->centerW,v->pcm_current-v->centerW,
 201                      0,0);
 202     _analysis_output("f2",seq,filtered2+v->centerW,v->pcm_current-v->centerW,
 203                      0,0);
 204     _analysis_output("f3",seq++,filtered3+v->centerW,v->pcm_current-v->centerW,
 205                      0,0);
 206
 207   }
 208
 209   ve->current=v->pcm_current;
 210
 211   /* Now search through our cached highpass data for breaking points */
 212   /* starting point */
 213   if(v->W)
 214     j=v->centerW+ci->blocksizes[1]/4-ci->blocksizes[0]/4;
 215   else
 216     j=v->centerW;
 217
 218   while(j+ve->winlength<=v->pcm_current){
 219     for(i=0;i<ve->ch;i++){
 220       for(k=0;k<4;k++){
 221         float *filtered=ve->filtered[i*4+k]+j;
 222         float m=_ve_deltai(ve,filtered-ve->winlength,filtered);
 223
 224         if(m>ci->preecho_thresh[k]){
 225           /*granulepos++;*/
 226           return(0);
 227         }
 228         if(m<ci->postecho_thresh[k]){
 229           /*granulepos++;*/
 230           return(0);
 231         }
 232         /*granulepos++;*/
 233       }
 234     }
 235
 236     j+=min(ci->blocksizes[0],ve->winlength)/2;
 237
 238     if(j>=searchpoint){
 239       return(1);
 240     }
 241   }
 242
 243   return(-1);
 244 }
 245
 246 void _ve_envelope_shift(envelope_lookup *e,long shift){
 247   int i;
 248   for(i=0;i<e->ch*4;i++)
 249     memmove(e->filtered[i],e->filtered[i]+shift,(e->current-shift)*
 250             sizeof(float));
 251   e->current-=shift;
 252 }
 253
 254