lib/psy.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: random psychoacoustics (not including preecho)
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Aug 26 1999
  18
  19  ********************************************************************/
  20
  21 #include <math.h>
  22 #include <string.h>
  23 #include "stdio.h"
  24 #include "codec.h"
  25 #include "psy.h"
  26 #include "lpc.h"
  27 #include "smallft.h"
  28 #include "xlogmap.h"
  29
  30 #define NOISEdB 0
  31
  32 #define MASKdB  18
  33 #define HROLL   60
  34 #define LROLL   90
  35 #define MASKBIAS  10
  36
  37 #define LNOISE  .95
  38 #define HNOISE  1.01
  39 #define NOISEBIAS  20
  40
  41 /* Find the mean log energy of a given 'band'; used to evaluate tones
  42    against background noise */
  43
  44 /* This is faster than a real convolution, gives us roughly the log f
  45    scale we seek, and gives OK results.  So, that means it's a good
  46    hack */
  47
  48 /* To add: f scale noise attenuation curve */
  49
  50 void _vp_noise_floor(double *f, double *m,int n){
  51   long lo=0,hi=0;
  52   double acc=0,div=0;
  53   int i,j;
  54
  55   for(i=100;i<n;i++){
  56     long newlo=i*LNOISE-NOISEBIAS;
  57     long newhi=i*HNOISE+NOISEBIAS;
  58     double temp;
  59
  60     if(newhi>n)newhi=n;
  61     if(newlo<0)newlo=0;
  62
  63     for(j=hi;j<newhi;j++){
  64       acc+=todB(f[j]);
  65       div++;
  66     }
  67     for(j=lo;j<newlo;j++){
  68       acc-=todB(f[j]);
  69       div--;
  70     }
  71
  72     hi=newhi;
  73     lo=newlo;
  74
  75     temp=fromdB(acc/div+NOISEdB); /* The NOISEdB constant should be an
  76                                      attenuation curve */
  77     if(m[i]<temp)m[i]=temp;
  78   }
  79 }
  80
  81 /* figure the masking curve.  linear rolloff on a dB scale, adjusted
  82    by octave */
  83 void _vp_mask_floor(double *f, double *m,int n){
  84   double ocSCALE=1./log(2);
  85   double curmask=-9.e40;
  86   double curoc=log(MASKBIAS)*ocSCALE;
  87   long i;
  88
  89   /* run mask forward then backward */
  90   for(i=0;i<n;i++){
  91     double newmask=todB(f[i])-MASKdB;
  92     double newoc=log(i+MASKBIAS)*ocSCALE;
  93     double roll=curmask-(newoc-curoc)*HROLL;
  94     double lroll;
  95     if(newmask>roll){
  96       roll=curmask=newmask;
  97       curoc=newoc;
  98     }
  99     lroll=fromdB(roll);
 100     if(m[i]<lroll)m[i]=lroll;
 101   }
 102
 103   curmask=-9.e40;
 104   curoc=log(n+MASKBIAS)*ocSCALE;
 105   for(i=n-1;i>=0;i--){
 106     double newmask=todB(f[i])-MASKdB;
 107     double newoc=log(i+MASKBIAS)*ocSCALE;
 108     double roll=curmask-(curoc-newoc)*LROLL;
 109     double lroll;
 110     if(newmask>roll){
 111       roll=curmask=newmask;
 112       curoc=newoc;
 113     }
 114     lroll=fromdB(roll);
 115     if(m[i]<lroll)m[i]=lroll;
 116   }
 117 }
 118
 119 /* s must be padded at the end with m-1 zeroes */
 120 static void time_convolve(double *s,double *r,int n,int m){
 121   int i;
 122
 123   for(i=0;i<n;i++){
 124     int j;
 125     double acc=0;
 126
 127     for(j=0;j<m;j++)
 128       acc+=s[i+j]*r[m-j-1];
 129
 130     s[i]=acc;
 131   }
 132 }
 133
 134 /*************************************************************************/
 135 /* Continuous balance analysis/synthesis *********************************/
 136
 137
 138 /* Compute the average continual spectral balance of the given vectors
 139    (in radians); encode into LPC coefficients */
 140
 141 double _vp_balance_compute(double *A, double *B, double *lpc,lpc_lookup *vb){
 142   /* correlate in time (the response function is small).  Log
 143      frequency scale, small mapping */
 144
 145   int n=vb->n;
 146   int mapped=vb->ln;
 147
 148   /* 256/15 are arbitrary but unimportant to decoding */
 149   int resp=15;
 150   int i;
 151
 152   /* This is encode side. Don't think too hard about it */
 153
 154   double workA[mapped+resp];
 155   double workB[mapped+resp];
 156   double p[mapped];
 157   double workC[resp];
 158
 159   memset(workA,0,sizeof(workA));
 160   memset(workB,0,sizeof(workB));
 161   memset(workC,0,sizeof(workC));
 162
 163   for(i=0;i<n;i++){
 164     int j=vb->bscale[i]+(resp>>1);
 165     double mag_sq=A[i]*A[i]+B[i]*B[i];
 166     double phi;
 167
 168     if(B[i]==0)
 169       phi=M_PI/2.;
 170     else{
 171       phi=atan(A[i]/B[i]);
 172       if((A[i]<0 && B[i]>0)||
 173          (A[i]>0 && B[i]<0)){
 174         /* rotate II and IV into the first quadrant.  III is already there */
 175         phi+=M_PI/2;
 176       }
 177     }
 178
 179     workA[j]+=mag_sq*sin(phi);
 180     workB[j]+=mag_sq*cos(phi);
 181   }
 182
 183   /* prepare convolution vector.  Play with a few different shapes */
 184
 185   for(i=0;i<resp;i++){
 186     workC[i]=sin(M_PI*(i+1)/(float)(resp+1));
 187     workC[i]*=workC[i];
 188   }
 189
 190   time_convolve(workA,workC,mapped,resp);
 191   time_convolve(workB,workC,mapped,resp);
 192
 193   {
 194     double amp1;
 195
 196     for(i=0;i<mapped;i++){
 197       p[i]=atan2(workA[i],workB[i]);
 198     }
 199
 200     amp1=sqrt(vorbis_gen_lpc(p,lpc,vb));
 201
 202     return(amp1);
 203   }
 204
 205 }
 206
 207 void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
 208                      lpc_lookup *vb,int divp){
 209   int i;
 210   for(i=0;i<vb->n;i++){
 211     double mag=sqrt(A[i]*A[i]+B[i]*B[i]);
 212     double del=vorbis_lpc_magnitude(vb->dscale[i],lpc,vb->m)*amp;
 213     double phi=atan2(A[i],B[i]);
 214
 215     if(divp)
 216       phi-=del;
 217     else
 218       phi+=del;
 219
 220     A[i]=mag*sin(phi);
 221     B[i]=mag*cos(phi);
 222   }
 223 }