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-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: random psychoacoustics (not including preecho)
  15  last mod: $Id: psy.c,v 1.8 1999/12/31 12:35:16 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include <stdio.h>
  23 #include "codec.h"
  24 #include "psy.h"
  25 #include "lpc.h"
  26 #include "smallft.h"
  27
  28 /* Set up decibel threshhold 'curves'.  Actually, just set a level at
  29    log frequency intervals, interpolate, and call it a curve. */
  30
  31 static void set_curve(double *points,
  32                       double *ref,int rn,double rrate,
  33                       double *c,int n, double crate){
  34   int i,j=0;
  35
  36   for(i=0;i<rn-1;i++){
  37     int endpos=points[i+1]*n*rrate/crate;
  38     double base=ref[i];
  39     double delta=(ref[i+1]-base)/(endpos-j);
  40     for(;j<endpos && j<n;j++){
  41       c[j]=base;
  42       base+=delta;
  43     }
  44   }
  45 }
  46
  47 void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n){
  48   memset(p,0,sizeof(psy_lookup));
  49   p->noisethresh=malloc(n*sizeof(double));
  50   p->maskthresh=malloc(n*sizeof(double));
  51   p->vi=vi;
  52   p->n=n;
  53
  54   /* set up the curves for a given blocksize and sample rate */
  55
  56   set_curve(vi->threshhold_points,
  57             vi->noisethresh,THRESH_POINTS,48000,p->noisethresh,n,vi->rate);
  58   set_curve(vi->threshhold_points,
  59             vi->maskthresh, THRESH_POINTS,48000,p->maskthresh, n,vi->rate);
  60
  61 #ifdef ANALYSIS
  62   {
  63     int j;
  64     FILE *out;
  65     char buffer[80];
  66
  67     sprintf(buffer,"noise_threshhold_%d.m",n);
  68     out=fopen(buffer,"w+");
  69     for(j=0;j<n;j++)
  70       fprintf(out,"%g\n",p->noisethresh[j]);
  71     fclose(out);
  72     sprintf(buffer,"mask_threshhold_%d.m",n);
  73     out=fopen(buffer,"w+");
  74     for(j=0;j<n;j++)
  75       fprintf(out,"%g\n",p->maskthresh[j]);
  76     fclose(out);
  77   }
  78 #endif
  79
  80 }
  81
  82 void _vp_psy_clear(psy_lookup *p){
  83   if(p){
  84     if(p->noisethresh)free(p->noisethresh);
  85     if(p->maskthresh)free(p->maskthresh);
  86     memset(p,0,sizeof(psy_lookup));
  87   }
  88 }
  89
  90 /* Find the mean log energy of a given 'band'; used to evaluate tones
  91    against background noise */
  92
  93 /* This is faster than a real convolution, gives us roughly the log f
  94    scale we seek, and gives OK results.  So, that means it's a good
  95    hack */
  96
  97 void _vp_noise_floor(psy_lookup *p, double *f, double *m){
  98   int n=p->n;
  99   vorbis_info *vi=p->vi;
 100
 101   long lo=0,hi=0;
 102   double acc=0.;
 103   double div=0.;
 104   int i,j;
 105
 106   double bias=n*vi->noisebias;
 107
 108   for(i=0;i<n;i++){
 109     long newlo=i*vi->lnoise-bias;
 110     long newhi=i*vi->hnoise+bias;
 111     double temp;
 112
 113     if(newhi>n)newhi=n;
 114     if(newlo<0)newlo=0;
 115
 116     for(j=hi;j<newhi;j++){
 117       acc+=todB(f[j]);
 118       div+=1.;
 119     }
 120     for(j=lo;j<newlo;j++){
 121       acc-=todB(f[j]);
 122       div-=1.;
 123     }
 124
 125     hi=newhi;
 126     lo=newlo;
 127
 128     /* attenuate by the noise threshhold curve */
 129     temp=fromdB(acc/div+p->noisethresh[i]);
 130     if(m[i]<temp)m[i]=temp;
 131   }
 132 }
 133
 134 /* Masking curve: linear rolloff on a dB scale, adjusted by octave,
 135    attenuated by maskthresh */
 136
 137 void _vp_mask_floor(psy_lookup *p,double *f, double *m){
 138   int n=p->n;
 139   double hroll=p->vi->hroll;
 140   double lroll=p->vi->lroll;
 141   double ocSCALE=1./log(2);
 142   double curmask=-9.e40;
 143   double maskbias=n*p->vi->maskbias;
 144   double curoc=log(maskbias)*ocSCALE;
 145   long i;
 146
 147   /* run mask forward then backward */
 148   for(i=0;i<n;i++){
 149     double newmask=todB(f[i])+p->maskthresh[i];
 150     double newoc=log(i+maskbias)*ocSCALE;
 151     double roll=curmask-(newoc-curoc)*hroll;
 152     double troll;
 153     if(newmask>roll){
 154       roll=curmask=newmask;
 155       curoc=newoc;
 156     }
 157     troll=fromdB(roll);
 158     if(m[i]<troll)m[i]=troll;
 159   }
 160
 161   curmask=-9.e40;
 162   curoc=log(n+maskbias)*ocSCALE;
 163   for(i=n-1;i>=0;i--){
 164     double newmask=todB(f[i])+p->maskthresh[i];
 165     double newoc=log(i+maskbias)*ocSCALE;
 166     double roll=curmask-(curoc-newoc)*lroll;
 167     double troll;
 168     if(newmask>roll){
 169       roll=curmask=newmask;
 170       curoc=newoc;
 171     }
 172     troll=fromdB(roll);
 173     if(m[i]<troll)m[i]=troll;
 174   }
 175 }
 176
 177 /* s must be padded at the end with m-1 zeroes */
 178 static void time_convolve(double *s,double *r,int n,int m){
 179   int i;
 180
 181   for(i=0;i<n;i++){
 182     int j;
 183     double acc=0;
 184
 185     for(j=0;j<m;j++)
 186       acc+=s[i+j]*r[m-j-1];
 187
 188     s[i]=acc;
 189   }
 190 }
 191