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