psy.c

OGG文件格式音频数解压缩SDK.现在OGG文件格式在一些游戏开发中使用的比较多.

/********************************************************************
 *                                                                  *
 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
 *                                                                  *
 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002             *
 * by the XIPHOPHORUS Company http://www.xiph.org/                  *
 *                                                                  *
 ********************************************************************

 function: psychoacoustics not including preecho
 last mod: $Id: psy.c,v 1.81 2002/10/21 07:00:11 xiphmont Exp $

 ********************************************************************/

#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "vorbis/codec.h"
#include "codec_internal.h"

#include "masking.h"
#include "psy.h"
#include "os.h"
#include "lpc.h"
#include "smallft.h"
#include "scales.h"
#include "misc.h"

#define NEGINF -9999.0f
static float stereo_threshholds[]={0.0f, .5f, 1.0f, 1.5f, 2.5f, 4.5f, 8.5f, 16.5f, (float)9e10};

vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
  codec_setup_info *ci=vi->codec_setup;
  vorbis_info_psy_global *gi=&ci->psy_g_param;
  vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));

  look->channels=vi->channels;

  look->ampmax=-9999.0f;
  look->gi=gi;
  return(look);
}

void _vp_global_free(vorbis_look_psy_global *look){
  if(look){
    memset(look,0,sizeof(*look));
    _ogg_free(look);
  }
}

void _vi_gpsy_free(vorbis_info_psy_global *i){
  if(i){
    memset(i,0,sizeof(*i));
    _ogg_free(i);
  }
}

void _vi_psy_free(vorbis_info_psy *i){
  if(i){
    memset(i,0,sizeof(*i));
    _ogg_free(i);
  }
}

static void min_curve(float *c, 
               float *c2){
  int i;  
  for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
}
static void max_curve(float *c,
		       float *c2){
  int i;  
  for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
}

static void attenuate_curve(float *c,float att){
  int i;
  for(i=0;i<EHMER_MAX;i++)
    c[i]+=att;
}

static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
                  float center_boost, float center_decay_rate){
  int i,j,k,m;
  float ath[EHMER_MAX];
  float workc[P_BANDS][P_LEVELS][EHMER_MAX];
  float athc[P_LEVELS][EHMER_MAX];
  float *brute_buffer=alloca(n*sizeof(*brute_buffer));

  float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS);

  memset(workc,0,sizeof(workc));

  for(i=0;i<P_BANDS;i++){
    /* we add back in the ATH to avoid low level curves falling off to
       -infinity and unnecessarily cutting off high level curves in the
       curve limiting (last step). */

    /* A half-band's settings must be valid over the whole band, and
       it's better to mask too little than too much */  
    int ath_offset=i*4;
    for(j=0;j<EHMER_MAX;j++){
      float min=999.0f;
      for(k=0;k<4;k++)
    if(j+k+ath_offset<MAX_ATH){
      if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
    }else{
      if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
    }
      ath[j]=min;
    }

    /* copy curves into working space, replicate the 50dB curve to 30
       and 40, replicate the 100dB curve to 110 */
    for(j=0;j<6;j++)
      memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j]));
    memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
    memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));

    /* apply centered curve boost/decay */
    for(j=0;j<P_LEVELS;j++){
      for(k=0;k<EHMER_MAX;k++){
    float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate;
    if(adj<0.0f && center_boost>0)adj=0.0f;
    if(adj>0.0f && center_boost<0)adj=0.0f;
    workc[i][j][k]+=adj;
      }
    }

    /* normalize curves so the driving amplitude is 0dB */
    /* make temp curves with the ATH overlayed */
    for(j=0;j<P_LEVELS;j++){
      attenuate_curve(workc[i][j], curveatt_dB[i]+100.0f-(j<2?2:j)*10.0f-P_LEVEL_0);
      memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc));
      attenuate_curve(athc[j],+100.0f-j*10.f-P_LEVEL_0);
      max_curve(athc[j],workc[i][j]);
    }

    /* Now limit the louder curves.

       the idea is this: We don't know what the playback attenuation
       will be; 0dB SL moves every time the user twiddles the volume
       knob. So that means we have to use a single 'most pessimal' curve
       for all masking amplitudes, right?  Wrong.  The *loudest* sound
       can be in (we assume) a range of ...+100dB] SL.  However, sounds
       20dB down will be in a range ...+80], 40dB down is from ...+60],
       etc... */

    for(j=1;j<P_LEVELS;j++){
      min_curve(athc[j],athc[j-1]);
      min_curve(workc[i][j],athc[j]);
    }
  }

  for(i=0;i<P_BANDS;i++){
    int hi_curve,lo_curve,bin;
    ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS);

    /* low frequency curves are measured with greater resolution than
       the MDCT/FFT will actually give us; we want the curve applied
       to the tone data to be pessimistic and thus apply the minimum
       masking possible for a given bin.  That means that a single bin
       could span more than one octave and that the curve will be a
       composite of multiple octaves.  It also may mean that a single
       bin may span > an eighth of an octave and that the eighth
       octave values may also be composited. */

    /* which octave curves will we be compositing? */
    bin      = (int)ogg_floor(fromOC(i * .5f)/binHz);
    lo_curve = (int)ogg_ceil(toOC(bin*binHz+1)*2);
    hi_curve = (int)ogg_floor(toOC((bin+1)*binHz)*2);
    if(lo_curve>i)lo_curve=i;
    if(lo_curve<0)lo_curve=0;
    if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1;

    for(m=0;m<P_LEVELS;m++){
      ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2));

      for(j=0;j<n;j++)brute_buffer[j]=999.0f;

      /* render the curve into bins, then pull values back into curve.
     The point is that any inherent subsampling aliasing results in
     a safe minimum */
      for(k=lo_curve;k<=hi_curve;k++){
    int l=0;

    for(j=0;j<EHMER_MAX;j++){
      int lo_bin= (int)(fromOC(j*.125f + k* .5f - 2.0625f)/binHz);
      int hi_bin= (int)(fromOC(j*.125f + k* .5f - 1.9375f)/binHz+1);

      if(lo_bin<0)lo_bin=0;
      if(lo_bin>n)lo_bin=n;
      if(lo_bin<l)l=lo_bin;
      if(hi_bin<0)hi_bin=0;
      if(hi_bin>n)hi_bin=n;

      for(;l<hi_bin && l<n;l++)
        if(brute_buffer[l]>workc[k][m][j])
          brute_buffer[l]=workc[k][m][j];
    }

    for(;l<n;l++)
      if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
        brute_buffer[l]=workc[k][m][EHMER_MAX-1];

      }

      /* be equally paranoid about being valid up to next half ocatve */
      if(i+1<P_BANDS){
    int l=0;
    k=i+1;
    for(j=0;j<EHMER_MAX;j++){
      int lo_bin= (int)(fromOC(j*.125f + i * .5f - 2.0625f)/binHz);
      int hi_bin= (int)(fromOC(j*.125f + i * .5f - 1.9375f)/binHz+1);

      if(lo_bin<0)lo_bin=0;
      if(lo_bin>n)lo_bin=n;
      if(lo_bin<l)l=lo_bin;
      if(hi_bin<0)hi_bin=0;
      if(hi_bin>n)hi_bin=n;

      for(;l<hi_bin && l<n;l++)
        if(brute_buffer[l]>workc[k][m][j])
          brute_buffer[l]=workc[k][m][j];
    }

    for(;l<n;l++)
      if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
        brute_buffer[l]=workc[k][m][EHMER_MAX-1];

      }


      for(j=0;j<EHMER_MAX;j++){
    int bin=(int)(fromOC(j*.125f+i*.5f-2.0f)/binHz);
    if(bin<0){
      ret[i][m][j+2]=-999.0f;
    }else{
      if(bin>=n){
        ret[i][m][j+2]=-999.0f;
      }else{
        ret[i][m][j+2]=brute_buffer[bin];
      }
    }
      }

      /* add fenceposts */
      for(j=0;j<EHMER_OFFSET;j++)
    if(ret[i][m][j+2]>-200.f)break;  
      ret[i][m][0]=(float)j;

      for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
    if(ret[i][m][j+2]>-200.f)
      break;
      ret[i][m][1]=(float)j;
    
    }
  }

  return(ret);
}

void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi, 
          vorbis_info_psy_global *gi,int n,ogg_int32_t rate){
  ogg_int32_t i,j,lo=-99,hi=0;
  ogg_int32_t maxoc;
  memset(p,0,sizeof(*p));

  p->eighth_octave_lines=gi->eighth_octave_lines;
  p->shiftoc=(ogg_int32_t)ogg_rint(ogg_log(gi->eighth_octave_lines*8.f)/log(2.f))-1;

  p->firstoc=(ogg_int32_t)toOC(.25f*rate*.5f/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
  maxoc=(ogg_int32_t)(toOC((n+.25f)*rate*.5f/n)*(1<<(p->shiftoc+1))+.5f);
  p->total_octave_lines=maxoc-p->firstoc+1;
  p->ath=_ogg_malloc(n*sizeof(*p->ath));

  p->octave=_ogg_malloc(n*sizeof(*p->octave));
  p->bark=_ogg_malloc(n*sizeof(*p->bark));
  p->vi=vi;
  p->n=n;
  p->rate=rate;

    /* set up the lookups for a given blocksize and sample rate */

  for(i=0,j=0;i<MAX_ATH-1;i++){
    int endpos=(int)ogg_rint(fromOC((i+1)*.125f - 2.0f)*2*n/rate);
    float base=ATH[i];
    if(j<endpos){
      float delta=(ATH[i+1]-base)/(endpos-j);
      for(;j<endpos && j<n;j++){
        p->ath[j]=base+100.0f;
        base+=delta;
      }
    }
  }

  for(i=0;i<n;i++){
    float bark=toBARK(rate/(2*n)*i); 

    for (;lo+vi->noisewindowlomin<i && 
      toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++){}
          
    for (;hi<n && (hi<i+vi->noisewindowhimin || 
      toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++){}
      
    p->bark[i]=((lo-1)<<16)+(hi-1);

  }

  for(i=0;i<n;i++)
    p->octave[i]=(ogg_int32_t)(toOC((i+.25f)*.5f*rate/n)*(1<<(p->shiftoc+1))+.5f);

  p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5f/n,n,
                  vi->tone_centerboost,vi->tone_decay);

  /* set up rolling noise median */
  p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
  for(i=0;i<P_NOISECURVES;i++)
    p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));

  for(i=0;i<n;i++){
    float halfoc=toOC((i+.5f)*rate/(2.0f*n))*2.0f;
    int inthalfoc;
    float del;

    if(halfoc<0)halfoc=0;
    if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
    inthalfoc=(int)halfoc;
    del=halfoc-inthalfoc;

    for(j=0;j<P_NOISECURVES;j++) 
      p->noiseoffset[j][i]= 
    p->vi->noiseoff[j][inthalfoc]*(1.0f-del) + 
    p->vi->noiseoff[j][inthalfoc+1]*del;

  }
}

void _vp_psy_clear(vorbis_look_psy *p){
  int i,j;
  if(p){
    if(p->ath)_ogg_free(p->ath);
    if(p->octave)_ogg_free(p->octave);
    if(p->bark)_ogg_free(p->bark);
    if(p->tonecurves){
      for(i=0;i<P_BANDS;i++){
	for(j=0;j<P_LEVELS;j++){
	  _ogg_free(p->tonecurves[i][j]);
	}
	_ogg_free(p->tonecurves[i]);
      }
      _ogg_free(p->tonecurves);
    }
    if(p->noiseoffset){
      for(i=0;i<P_NOISECURVES;i++){
        _ogg_free(p->noiseoffset[i]);
      }
      _ogg_free(p->noiseoffset);
    }
    memset(p,0,sizeof(*p));
  }
}

/* octave/(8*eighth_octave_lines) x scale and dB y scale */
static void seed_curve(float *seed,
		       const float **curves,
		       float amp,
		       int oc, int n,
		       int linesper,float dBoffset){
  int i,post1;