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: psychoacoustics not including preecho
  15  last mod: $Id: psy.c,v 1.18 2000/04/03 08:30:49 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include <stdio.h>
  23 #include "vorbis/codec.h"
  24
  25 #include "psy.h"
  26 #include "os.h"
  27 #include "lpc.h"
  28 #include "smallft.h"
  29 #include "scales.h"
  30 #include "misc.h"
  31
  32 /* Set up decibel threshhold slopes on a Bark frequency scale */
  33
  34 static void set_curve(double *ref,double *c,int n, double crate){
  35   int i,j=0;
  36
  37   for(i=0;i<MAX_BARK-1;i++){
  38     int endpos=rint(fromBARK(i+1)*2*n/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(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate){
  49   long i;
  50   memset(p,0,sizeof(vorbis_look_psy));
  51   p->maskthresh=malloc(n*sizeof(double));
  52   p->barknum=malloc(n*sizeof(double));
  53   p->vi=vi;
  54   p->n=n;
  55
  56   /* set up the lookups for a given blocksize and sample rate */
  57   /* Vorbis max sample rate is limited by 26 Bark (54kHz) */
  58   set_curve(vi->maskthresh, p->maskthresh, n,rate);
  59
  60   for(i=0;i<n;i++)
  61     p->barknum[i]=toBARK(rate/2.*i/n);
  62
  63 #ifdef ANALYSIS
  64   {
  65     int j;
  66     FILE *out;
  67     char buffer[80];
  68
  69     sprintf(buffer,"mask_threshhold_%d.m",n);
  70     out=fopen(buffer,"w+");
  71     for(j=0;j<n;j++)
  72       fprintf(out,"%g\n",p->maskthresh[j]);
  73     fclose(out);
  74   }
  75 #endif
  76
  77 }
  78
  79 void _vp_psy_clear(vorbis_look_psy *p){
  80   if(p){
  81     if(p->maskthresh)free(p->maskthresh);
  82     if(p->barknum)free(p->barknum);
  83     memset(p,0,sizeof(vorbis_look_psy));
  84   }
  85 }
  86
  87 /* Masking curve: linear rolloff on a Bark/dB scale, attenuated by
  88    maskthresh */
  89
  90 void _vp_mask_floor(vorbis_look_psy *p,double *f, double *floor,int attp){
  91   int n=p->n;
  92   double hroll=p->vi->hrolldB;
  93   double lroll=p->vi->lrolldB;
  94   double curmask=todB(f[0])+(attp?p->maskthresh[0]:0);
  95   double curoc=0.;
  96   long i;
  97
  98   /* run mask forward then backward */
  99   for(i=0;i<n;i++){
 100     double newmask=todB(f[i])+(attp?p->maskthresh[i]:0);
 101     double newoc=p->barknum[i];
 102     double roll=curmask-(newoc-curoc)*hroll;
 103     double troll;
 104     if(newmask>roll){
 105       roll=curmask=newmask;
 106       curoc=newoc;
 107     }
 108     troll=fromdB(roll);
 109     if(floor[i]<troll)floor[i]=troll;
 110   }
 111
 112   curmask=todB(f[n-1])+(attp?p->maskthresh[n-1]:0);
 113   curoc=p->barknum[n-1];
 114   for(i=n-1;i>=0;i--){
 115     double newmask=todB(f[i])+(attp?p->maskthresh[i]:0);
 116     double newoc=p->barknum[i];
 117     double roll=curmask-(curoc-newoc)*lroll;
 118     double troll;
 119     if(newmask>roll){
 120       roll=curmask=newmask;
 121       curoc=newoc;
 122     }
 123     troll=fromdB(roll);
 124     if(floor[i]<troll)floor[i]=troll;
 125   }
 126 }
 127
 128 /* s must be padded at the end with m-1 zeroes */
 129 static void time_convolve(double *s,double *r,int n,int m){
 130   int i;
 131
 132   for(i=0;i<n;i++){
 133     int j;
 134     double acc=0;
 135
 136     for(j=0;j<m;j++)
 137       acc+=s[i+j]*r[m-j-1];
 138
 139     s[i]=acc;
 140   }
 141 }
 142
 143 void _vi_psy_free(vorbis_info_psy *i){
 144   if(i){
 145     memset(i,0,sizeof(vorbis_info_psy));
 146     free(i);
 147   }
 148 }