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.9 2000/01/04 09:05:02 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 "scales.h"
  28
  29 /* Set up decibel threshhold slopes on a Bark frequency scale */
  30
  31 static void set_curve(double *ref,double *c,int n, double crate){
  32   int i,j=0;
  33
  34   for(i=0;i<MAX_BARK-1;i++){
  35     int endpos=rint(fromBARK(i+1)*2*n/crate);
  36     double base=ref[i];
  37     double delta=(ref[i+1]-base)/(endpos-j);
  38     for(;j<endpos && j<n;j++){
  39       c[j]=base;
  40       base+=delta;
  41     }
  42   }
  43 }
  44
  45 void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n){
  46   long i;
  47   memset(p,0,sizeof(psy_lookup));
  48   p->maskthresh=malloc(n*sizeof(double));
  49   p->barknum=malloc(n*sizeof(double));
  50   p->vi=vi;
  51   p->n=n;
  52
  53   /* set up the lookups for a given blocksize and sample rate */
  54   /* Vorbis max sample rate is limited by 26 Bark (54kHz) */
  55   set_curve(vi->maskthresh, p->maskthresh, n,vi->rate);
  56
  57   for(i=0;i<n;i++)
  58     p->barknum[i]=toBARK(vi->rate/2.*i/n);
  59
  60 #ifdef ANALYSIS
  61   {
  62     int j;
  63     FILE *out;
  64     char buffer[80];
  65
  66     sprintf(buffer,"mask_threshhold_%d.m",n);
  67     out=fopen(buffer,"w+");
  68     for(j=0;j<n;j++)
  69       fprintf(out,"%g\n",p->maskthresh[j]);
  70     fclose(out);
  71   }
  72 #endif
  73
  74 }
  75
  76 void _vp_psy_clear(psy_lookup *p){
  77   if(p){
  78     if(p->maskthresh)free(p->maskthresh);
  79     if(p->barknum)free(p->barknum);
  80     memset(p,0,sizeof(psy_lookup));
  81   }
  82 }
  83
  84 /* Masking curve: linear rolloff on a Bark/dB scale, attenuated by
  85    maskthresh */
  86
  87 void _vp_mask_floor(psy_lookup *p,double *f, double *m){
  88   int n=p->n;
  89   double hroll=p->vi->hrolldB;
  90   double lroll=p->vi->lrolldB;
  91   double curmask=todB(f[0])+p->maskthresh[0];
  92   double curoc=0.;
  93   long i;
  94
  95   /* run mask forward then backward */
  96   for(i=0;i<n;i++){
  97     double newmask=todB(f[i])+p->maskthresh[i];
  98     double newoc=p->barknum[i];
  99     double roll=curmask-(newoc-curoc)*hroll;
 100     double troll;
 101     if(newmask>roll){
 102       roll=curmask=newmask;
 103       curoc=newoc;
 104     }
 105     troll=fromdB(roll);
 106     if(m[i]<troll)m[i]=troll;
 107   }
 108
 109   curmask=todB(f[n-1])+p->maskthresh[n-1];
 110   curoc=p->barknum[n-1];
 111   for(i=n-1;i>=0;i--){
 112     double newmask=todB(f[i])+p->maskthresh[i];
 113     double newoc=p->barknum[i];
 114     double roll=curmask-(curoc-newoc)*lroll;
 115     double troll;
 116     if(newmask>roll){
 117       roll=curmask=newmask;
 118       curoc=newoc;
 119     }
 120     troll=fromdB(roll);
 121     if(m[i]<troll)m[i]=troll;
 122   }
 123 }
 124
 125 /* s must be padded at the end with m-1 zeroes */
 126 static void time_convolve(double *s,double *r,int n,int m){
 127   int i;
 128
 129   for(i=0;i<n;i++){
 130     int j;
 131     double acc=0;
 132
 133     for(j=0;j<m;j++)
 134       acc+=s[i+j]*r[m-j-1];
 135
 136     s[i]=acc;
 137   }
 138 }
 139