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 10
  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 void _vp_psy_quantize(double *f, double *m,int n){
 120   int i;
 121   for(i=0;i<n;i++){
 122     int val=rint(f[i]/m[i]);
 123     if(val>16)val=16;
 124     if(val<-16)val=-16;
 125     if(val==0 || val==2){
 126       if(f[i]<0){
 127         f[i]=-1;
 128       }else{
 129         f[i]=1;
 130       }
 131     }
 132
 133     f[i]=val;
 134   }
 135 }
 136 void _vp_psy_unquantize(double *f, double *m,int n){
 137   int i;
 138   for(i=0;i<n;i++)
 139     f[i]*=m[i];
 140 }
 141
 142 void _vp_psy_sparsify(double *f, double *m,int n){
 143   int i;
 144   for(i=0;i<n;i++)
 145     if(fabs(f[i])<m[i]*.5)f[i]=0;
 146 }
 147
 148 /* s must be padded at the end with m-1 zeroes */
 149 static void time_convolve(double *s,double *r,int n,int m){
 150   int i;
 151
 152   for(i=0;i<n;i++){
 153     int j;
 154     double acc=0;
 155
 156     for(j=0;j<m;j++)
 157       acc+=s[i+j]*r[m-j-1];
 158
 159     s[i]=acc;
 160   }
 161 }
 162
 163 /*************************************************************************/
 164 /* Continuous balance analysis/synthesis *********************************/
 165
 166
 167 /* Compute the average continual spectral balance of the given vectors
 168    (in radians); encode into LPC coefficients */
 169
 170 double _vp_balance_compute(double *A, double *B, double *lpc,lpc_lookup *vb){
 171   /* correlate in time (the response function is small).  Log
 172      frequency scale, small mapping */
 173
 174   int n=vb->n;
 175   int mapped=vb->ln;
 176
 177   /* 256/15 are arbitrary but unimportant to decoding */
 178   int resp=15;
 179   int i;
 180
 181   /* This is encode side. Don't think too hard about it */
 182
 183   double workA[mapped+resp];
 184   double workB[mapped+resp];
 185   double p[mapped];
 186   double workC[resp];
 187
 188   memset(workA,0,sizeof(workA));
 189   memset(workB,0,sizeof(workB));
 190   memset(workC,0,sizeof(workC));
 191
 192   for(i=0;i<n;i++){
 193     int j=vb->bscale[i]+(resp>>1);
 194     double mag_sq=A[i]*A[i]+B[i]*B[i];
 195     double phi;
 196
 197     if(B[i]==0)
 198       phi=M_PI/2.;
 199     else{
 200       phi=atan(A[i]/B[i]);
 201       if((A[i]<0 && B[i]>0)||
 202          (A[i]>0 && B[i]<0)){
 203         /* rotate II and IV into the first quadrant.  III is already there */
 204         phi+=M_PI/2;
 205       }
 206     }
 207
 208     workA[j]+=mag_sq*sin(phi);
 209     workB[j]+=mag_sq*cos(phi);
 210   }
 211
 212   /* prepare convolution vector.  Play with a few different shapes */
 213
 214   for(i=0;i<resp;i++){
 215     workC[i]=sin(M_PI*(i+1)/(float)(resp+1));
 216     workC[i]*=workC[i];
 217   }
 218
 219   time_convolve(workA,workC,mapped,resp);
 220   time_convolve(workB,workC,mapped,resp);
 221
 222   {
 223     double amp1;
 224
 225     for(i=0;i<mapped;i++){
 226       p[i]=atan2(workA[i],workB[i]);
 227     }
 228
 229     amp1=sqrt(vorbis_gen_lpc(p,lpc,vb));
 230
 231     return(amp1);
 232   }
 233
 234 }
 235
 236 void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
 237                      lpc_lookup *vb,int divp){
 238   int i;
 239   for(i=0;i<vb->n;i++){
 240     double mag=sqrt(A[i]*A[i]+B[i]*B[i]);
 241     double del=vorbis_lpc_magnitude(vb->dscale[i],lpc,vb->m)*amp;
 242     double phi=atan2(A[i],B[i]);
 243
 244     if(divp)
 245       phi-=del;
 246     else
 247       phi+=del;
 248
 249     A[i]=mag*sin(phi);
 250     B[i]=mag*cos(phi);
 251   }
 252 }