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