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.15 2000/02/06 13:39:44 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 "lpc.h"
  27 #include "smallft.h"
  28 #include "scales.h"
  29
  30 /* Set up decibel threshhold slopes on a Bark frequency scale */
  31
  32 static void set_curve(double *ref,double *c,int n, double crate){
  33   int i,j=0;
  34
  35   for(i=0;i<MAX_BARK-1;i++){
  36     int endpos=rint(fromBARK(i+1)*2*n/crate);
  37     double base=ref[i];
  38     double delta=(ref[i+1]-base)/(endpos-j);
  39     for(;j<endpos && j<n;j++){
  40       c[j]=base;
  41       base+=delta;
  42     }
  43   }
  44 }
  45
  46 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate){
  47   long i;
  48   memset(p,0,sizeof(vorbis_look_psy));
  49   p->maskthresh=malloc(n*sizeof(double));
  50   p->barknum=malloc(n*sizeof(double));
  51   p->vi=vi;
  52   p->n=n;
  53
  54   /* set up the lookups for a given blocksize and sample rate */
  55   /* Vorbis max sample rate is limited by 26 Bark (54kHz) */
  56   set_curve(vi->maskthresh, p->maskthresh, n,rate);
  57
  58   for(i=0;i<n;i++)
  59     p->barknum[i]=toBARK(rate/2.*i/n);
  60
  61 #ifdef ANALYSIS
  62   {
  63     int j;
  64     FILE *out;
  65     char buffer[80];
  66
  67     sprintf(buffer,"mask_threshhold_%d.m",n);
  68     out=fopen(buffer,"w+");
  69     for(j=0;j<n;j++)
  70       fprintf(out,"%g\n",p->maskthresh[j]);
  71     fclose(out);
  72   }
  73 #endif
  74
  75 }
  76
  77 void _vp_psy_clear(vorbis_look_psy *p){
  78   if(p){
  79     if(p->maskthresh)free(p->maskthresh);
  80     if(p->barknum)free(p->barknum);
  81     memset(p,0,sizeof(vorbis_look_psy));
  82   }
  83 }
  84
  85 /* Masking curve: linear rolloff on a Bark/dB scale, attenuated by
  86    maskthresh */
  87 /* right now, floor==mask */
  88 void _vp_mask_floor(vorbis_look_psy *p,double *f, double *mask,double *floor){
  89   int n=p->n;
  90   double hroll=p->vi->hrolldB;
  91   double lroll=p->vi->lrolldB;
  92   double curmask=todB(f[0])+p->maskthresh[0];
  93   double curoc=0.;
  94   long i;
  95
  96   /* run mask forward then backward */
  97   for(i=0;i<n;i++){
  98     double newmask=todB(f[i])+p->maskthresh[i];
  99     double newoc=p->barknum[i];
 100     double roll=curmask-(newoc-curoc)*hroll;
 101     double troll;
 102     if(newmask>roll){
 103       roll=curmask=newmask;
 104       curoc=newoc;
 105     }
 106     troll=fromdB(roll);
 107     if(mask[i]<troll)mask[i]=troll;
 108     if(floor[i]<troll)floor[i]=troll;
 109   }
 110
 111   curmask=todB(f[n-1])+p->maskthresh[n-1];
 112   curoc=p->barknum[n-1];
 113   for(i=n-1;i>=0;i--){
 114     double newmask=todB(f[i])+p->maskthresh[i];
 115     double newoc=p->barknum[i];
 116     double roll=curmask-(curoc-newoc)*lroll;
 117     double troll;
 118     if(newmask>roll){
 119       roll=curmask=newmask;
 120       curoc=newoc;
 121     }
 122     troll=fromdB(roll);
 123     if(mask[i]<troll)mask[i]=troll;
 124     if(floor[i]<troll)floor[i]=troll;
 125   }
 126 }
 127
 128 /* take a masking curve and raw residue; eliminate the inaduble and
 129    quantize to the final form handed to the VQ.  All and any tricks to
 130    squeeze out bits given knowledge of the encoding mode should go
 131    here too */
 132
 133 /* modifies the pcm vector, returns book membership in aux */
 134
 135 void _vp_quantize(vorbis_look_psy *p, double *pcm, double *mask,
 136                   double *floor,int *aux,long begin, long n,long subn){
 137   long i,j;
 138     /* for now, we're not worrying about subvector, but the idea is
 139        that we normal blocks not have zeroes; zeroes only exist as an
 140        all-zero block */
 141
 142   for(i=begin;i<n;i++){
 143     double value=rint(pcm[i]/floor[i]);
 144     if(value>15)value=15;
 145     if(value<-15)value=-15;
 146     pcm[i]=value;
 147   }
 148
 149 }
 150
 151 /* s must be padded at the end with m-1 zeroes */
 152 static void time_convolve(double *s,double *r,int n,int m){
 153   int i;
 154
 155   for(i=0;i<n;i++){
 156     int j;
 157     double acc=0;
 158
 159     for(j=0;j<m;j++)
 160       acc+=s[i+j]*r[m-j-1];
 161
 162     s[i]=acc;
 163   }
 164 }
 165
 166 void _vi_psy_free(vorbis_info_psy *i){
 167   if(i){
 168     memset(i,0,sizeof(vorbis_info_psy));
 169     free(i);
 170   }
 171 }