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-2000             *
   9  * by 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  last mod: $Id: envelope.c,v 1.19 2000/06/18 12:33:47 xiphmont Exp $
  16
  17  Preecho calculation.
  18
  19  ********************************************************************/
  20
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include <stdio.h>
  24 #include <math.h>
  25 #include "vorbis/codec.h"
  26
  27 #include "os.h"
  28 #include "scales.h"
  29 #include "smallft.h"
  30 #include "envelope.h"
  31 #include "bitwise.h"
  32 #include "window.h"
  33 #include "misc.h"
  34
  35 void _ve_envelope_init(envelope_lookup *e,int samples_per){
  36   int i;
  37   e->winlen=samples_per*2;
  38   e->window=malloc(e->winlen*sizeof(double));
  39
  40   e->fft=calloc(1,sizeof(drft_lookup));
  41   drft_init(e->fft,samples_per*2);
  42
  43   /* We just use a straight sin(x) window for this */
  44   for(i=0;i<e->winlen;i++)
  45     e->window[i]=sin((i+.5)/e->winlen*M_PI);
  46 }
  47
  48 void _ve_envelope_clear(envelope_lookup *e){
  49   drft_clear(e->fft);
  50   free(e->fft);
  51   if(e->window)free(e->window);
  52   memset(e,0,sizeof(envelope_lookup));
  53 }
  54
  55 static void smooth_noise(long n,double *f,double *noise){
  56   long i;
  57   long lo=0,hi=0;
  58   double acc=0.;
  59
  60   for(i=0;i<n;i++){
  61     /* not exactly correct, (the center frequency should be centered
  62        on a *log* scale), but not worth quibbling */
  63     long newhi=i*1.0442718740+5;
  64     long newlo=i*.8781245150-5;
  65     if(newhi>n)newhi=n;
  66
  67     for(;lo<newlo;lo++)
  68       acc-=todB(f[lo]); /* yeah, this ain't RMS */
  69     for(;hi<newhi;hi++)
  70       acc+=todB(f[hi]);
  71     noise[i]=acc/(hi-lo);
  72   }
  73 }
  74
  75 /* use FFT for spectral power estimation */
  76 static int frameno=0;
  77 static int frameno2=0;
  78
  79 static void _ve_deltas(double *deltas,double *pcm,int n,double *window,
  80                        int samples_per,drft_lookup *l){
  81   int i,j;
  82   double *out=alloca(sizeof(double)*samples_per*2);
  83   double *cache=alloca(sizeof(double)*samples_per*2);
  84
  85   for(j=-1;j<n;j++){
  86
  87     memcpy(out,pcm+(j+1)*samples_per,samples_per*2*sizeof(double));
  88     for(i=0;i<samples_per*2;i++)
  89       out[i]*=window[i];
  90
  91     _analysis_output("Dpcm",frameno*1000+frameno2,out,samples_per*2,0,0);
  92
  93
  94     drft_forward(l,out);
  95     for(i=1;i<samples_per;i++)
  96       out[i]=hypot(out[i*2],out[i*2-1]);
  97     _analysis_output("Dfft",frameno*1000+frameno2,out,samples_per,0,1);
  98     smooth_noise(samples_per,out,out+samples_per);
  99
 100     if(j==-1){
 101       for(i=samples_per/10;i<samples_per;i++)
 102         cache[i]=out[i+samples_per];
 103     }else{
 104       double max=0;
 105       _analysis_output("Dcache",frameno*1000+frameno2,cache,samples_per,0,0);
 106       for(i=samples_per/10;i<samples_per;i++){
 107         double val=out[i+samples_per]-cache[i];
 108         cache[i]=out[i+samples_per];
 109         if(val>0)max+=val;
 110       }
 111       max/=samples_per;
 112       if(deltas[j]<max)deltas[j]=max;
 113     }
 114     _analysis_output("Dnoise",frameno*1000+frameno2++,out+samples_per,samples_per,0,0);
 115   }
 116 }
 117
 118 void _ve_envelope_deltas(vorbis_dsp_state *v){
 119   vorbis_info *vi=v->vi;
 120   int step=vi->envelopesa;
 121
 122   int dtotal=v->pcm_current/vi->envelopesa-1;
 123   int dcurr=v->envelope_current;
 124   int pch;
 125
 126   if(dtotal>dcurr){
 127     double *mult=v->multipliers+dcurr;
 128     memset(mult,0,sizeof(double)*(dtotal-dcurr));
 129
 130     for(pch=0;pch<vi->channels;pch++){
 131       double *pcm=v->pcm[pch]+(dcurr-1)*step;
 132       _ve_deltas(mult,pcm,dtotal-dcurr,v->ve.window,step,v->ve.fft);
 133
 134       {
 135         double *multexp=alloca(sizeof(double)*v->pcm_current);
 136         int i,j,k;
 137
 138         memset(multexp,0,sizeof(double)*v->pcm_current);
 139         j=0;
 140         for(i=0;i<dtotal;i++)
 141           for(k=0;k<step;k++)
 142             multexp[j++]=v->multipliers[i];
 143
 144         _analysis_output("Apcm",frameno,v->pcm[pch],v->pcm_current,0,0);
 145         _analysis_output("Amult",frameno++,multexp,v->pcm_current,0,0);
 146       }
 147
 148     }
 149     v->envelope_current=dtotal;
 150   }
 151 }
 152
 153
 154
 155
 156
 157