# vorbis makefile configured for use with gcc on any platform
-# $Id: Makefile.in,v 1.11 1999/10/10 20:32:17 xiphmont Exp $
+# $Id: Makefile.in,v 1.12 1999/10/18 09:41:07 xiphmont Exp $
###############################################################################
# #
RANLIB=@RANLIB@
LIBS=@LIBS@ -lm
-HFILES = mdct.h codec.h bitwise.h envelope.h lpc.h lsp.h modes.h\
+HFILES = mdct.h codec.h bitwise.h envelope.h lpc.h lsp.h \
psy.h smallft.h window.h xlogmap.h
LFILES = framing.o mdct.o smallft.o block.o envelope.o window.o\
debug:
$(MAKE) target CFLAGS="$(DEBUG)"
+analysis:
+ $(MAKE) target CFLAGS="$(DEBUG) -DANALYSIS"
+
profile:
$(MAKE) target CFLAGS="$(PROFILE)"
$(LFILES): $(HFILES)
+info.o: modes.h
+
.c.o:
$(CC) $(CFLAGS) -c $<
function: single-block PCM analysis
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 7 1999
+ last modification date: Oct 15 1999
********************************************************************/
we go... --Monty 19991004 */
int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
- static int frameno=0;
-
int i;
double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+ psy_lookup *vp=&vb->vd->vp[vb->W];
lpc_lookup *vl=&vb->vd->vl[vb->W];
vorbis_dsp_state *vd=vb->vd;
vorbis_info *vi=vd->vi;
int n=vb->pcmend;
int spectral_order=vi->floororder[vb->W];
+ vb->gluebits=0;
+ vb->time_envelope_bits=0;
+ vb->spectral_envelope_bits=0;
+ vb->spectral_residue_bits=0;
+
/*lpc_lookup *vbal=&vb->vd->vbal[vb->W];
double balance_v[vbal->m];
double balance_amp;*/
/* just do by channel. No coupling yet */
{
for(i=0;i<vi->channels;i++){
+ static int frameno=0;
+ int j;
double floor[n/2];
double curve[n/2];
double *lpc=vb->lpc[i];
memset(floor,0,sizeof(double)*n/2);
- _vp_noise_floor(vb->pcm[i],floor,n/2);
- _vp_mask_floor(vb->pcm[i],floor,n/2);
+ _vp_noise_floor(vp,vb->pcm[i],floor);
+
+#ifdef ANALYSIS
+ {
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"spectrum.m");
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",vb->pcm[i][j]);
+ fclose(out);
+
+ sprintf(buffer,"noise.m");
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",floor[j]);
+ fclose(out);
+ }
+#endif
+
+ _vp_mask_floor(vp,vb->pcm[i],floor);
+
+#ifdef ANALYSIS
+ {
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"premask.m");
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",floor[j]);
+ fclose(out);
+ }
+#endif
/* Convert our floor to a set of lpc coefficients */
vb->amp[i]=sqrt(vorbis_curve_to_lpc(floor,lpc,vl));
/* this may do various interesting massaging too...*/
_vs_residue_quantize(vb->pcm[i],curve,vi,n/2);
+#ifdef ANALYSIS
+ {
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"mask.m");
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",curve[j]);
+ fclose(out);
+
+ sprintf(buffer,"res.m");
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",vb->pcm[i][j]);
+ fclose(out);
+ }
+#endif
+
/* encode the residue */
_vs_residue_encode(vb,vb->pcm[i]);
function: PCM data vector blocking, windowing and dis/reassembly
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 2 1999
+ last modification date: Oct 15 1999
Handle windowing, overlap-add, etc of the PCM vectors. This is made
more amusing by Vorbis' current two allowed block sizes.
#include "envelope.h"
#include "mdct.h"
#include "lpc.h"
+#include "bitwise.h"
+#include "psy.h"
/* pcm accumulator examples (not exhaustive):
/* the coder init is different for read/write */
v->analysisp=1;
+ _vp_psy_init(&v->vp[0],vi,vi->blocksize[0]/2);
+ _vp_psy_init(&v->vp[1],vi,vi->blocksize[1]/2);
/* Yes, wasteful to have four lookups. This will get collapsed once
things crystallize */
- lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/8,
+
+ lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
vi->floororder[0],vi->flooroctaves[0],1);
- lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/8,
+ lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/2,
vi->floororder[0],vi->flooroctaves[0],1);
/*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
mdct_clear(&v->vm[1]);
lpc_clear(&v->vl[0]);
lpc_clear(&v->vl[1]);
+ _vp_psy_clear(&v->vp[0]);
+ _vp_psy_clear(&v->vp[1]);
/*lpc_clear(&v->vbal[0]);
lpc_clear(&v->vbal[1]);*/
memset(v,0,sizeof(vorbis_dsp_state));
memcpy(vb->mult[i],v->multipliers[i]+beginM,vi->blocksize[v->W]/
vi->envelopesa*sizeof(double));
- vb->frameno=v->frame;
+ vb->sequence=v->sequence;
+ vb->frameno=v->frameno;
/* handle eof detection: eof==0 means that we've not yet received EOF
eof>0 marks the last 'real' sample in pcm[]
v->W=v->nW;
v->centerW=new_centerNext;
- v->frame++;
- v->samples+=movementW;
+ v->sequence++;
+ v->frameno+=movementW;
if(v->eofflag)
v->eofflag-=movementW;
/* Yes, wasteful to have four lookups. This will get collapsed once
things crystallize */
- lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/8,
+ lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
vi->floororder[0],vi->flooroctaves[0],0);
- lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/8,
+ lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[0]/2,
vi->floororder[1],vi->flooroctaves[1],0);
/*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
vi->balanceorder,vi->balanceoctaves,0);
/* Shift out any PCM that we returned previously */
+ v->sequence++;
if(v->pcm_returned && v->centerW>vi->blocksize[1]/2){
/* don't shift too much; we need to have a minimum PCM buffer of
v->lW=v->W;
v->W=vb->W;
+ v->gluebits+=vb->gluebits;
+ v->time_envelope_bits+=vb->time_envelope_bits;
+ v->spectral_envelope_bits+=vb->spectral_envelope_bits;
+ v->spectral_residue_bits+=vb->spectral_residue_bits;
+
{
int sizeW=vi->blocksize[v->W];
int centerW=v->centerW+vi->blocksize[v->lW]/4+sizeW/4;
typedef struct {
int n;
+ struct vorbis_info *vi;
+
+ double *noisethresh;
+ double *maskthresh;
+
+} psy_lookup;
+
+typedef struct {
+ int n;
int log2n;
double *trig;
compression/decompression mode in progress (eg, psychoacoustic settings,
channel setup, options, codebook etc) *********************************/
+#define THRESH_POINTS 20
+
typedef struct vorbis_info{
int channels;
int rate;
int balanceoctaves;
this would likely need to be by channel and blocksize anyway*/
+ /* Encode-side settings for analysis */
+
double preecho_thresh;
double preecho_clamp;
+ double *threshhold_points;
+ double noisethresh[THRESH_POINTS];
+ double lnoise;
+ double hnoise;
+ double noisebias;
+ double maskthresh[THRESH_POINTS];
+ double lroll;
+ double hroll;
+ double maskbias;
+
/* local storage, only used on the encoding size. This way the
application does not need to worry about freeing some packets'
memory and not others'. Packet storage is always tracked */
mdct_lookup vm[2];
lpc_lookup vl[2];
lpc_lookup vbal[2];
+ psy_lookup vp[2];
double **pcm;
double **pcmret;
long nW;
long centerW;
- long frame;
- long samples;
+ long frameno;
+ long sequence;
size64 gluebits;
size64 time_envelope_bits;
int eofflag;
int frameno;
+ int sequence;
vorbis_dsp_state *vd; /* For read-only access of configuration */
long gluebits;
--- /dev/null
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
+ * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
+ * PLEASE READ THESE TERMS DISTRIBUTING. *
+ * *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
+ * http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: coder backend; handles encoding into and decoding from
+ codewords, codebook generation, and codebook transmission
+ author: Monty <xiphmont@mit.edu>
+ modifications by: Monty
+ last modification date: Oct 13 1999
+
+ ********************************************************************/
+
+#include <stdlib.h>
+
+#include "codec.h"
+
+/* several basic functions:
+
+ create codebooks from representation
+ create codebooks from data
+ encode codebooks
+ decode codebooks
+
+ encode using codebook
+ decode using codebook
+ map an entry to a sparse codebook's (eg, VQ) codeword index
+*/
+
+typedef struct vorbis_codebook{
+
+ int entries;
+
+ /*** codebook side ****************************************************/
+
+ /* encode side tree structure (indice->codeword) */
+ int *codewords; /* null if indice==codeword (ie fixed len codeword) */
+ int *codelengths; /* null if indice==codeword */
+ int len; /* -1 if varlength */
+
+ /* decode side tree structure (codeword->indice) */
+
+ int *tree; /* null if codeword==indice */
+
+ /*** mapping side *****************************************************/
+
+ int wordsper;
+
+ /* decode is easy */
+ int remaining;
+ int *mapping;
+
+ /* encode in the VQ case is what's hard. We need to find the
+ 'closest match' out of a large n-dimentional space */
+
+
+} vorbis_codebook;
while((samples=vorbis_synthesis_pcmout(&vd,&pcm))>0){
int j;
+ int clipflag=0;
int out=(samples<convsize?samples:convsize);
/* convert doubles to 16 bit signed ints (host order) and
for(j=0;j<out;j++){
int val=rint(mono[j]*32767.);
/* might as well guard clipping */
- if(val>32767)val=32767;
- if(val<-32768)val=-32768;
+ if(val>32767){
+ val=32767;
+ clipflag=1;
+ }
+ if(val<-32768){
+ val=-32768;
+ clipflag=1;
+ }
*ptr=val;
ptr+=2;
}
}
+
+ if(clipflag)
+ fprintf(stderr,"Clipping in frame %ld\n",vd.sequence);
+
+
fwrite(convbuffer,2*vi.channels,out,stdout);
vorbis_synthesis_read(&vd,out); /* tell libvorbis how
function: PCM data envelope analysis and manipulation
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 06 1999
+ last modification date: Oct 17 1999
Vorbis manipulates the dynamic range of the incoming PCM data
envelope to minimise time-domain energy leakage from percussive and
#include <math.h>
#include "codec.h"
+#include "mdct.h"
#include "envelope.h"
#include "bitwise.h"
return(1);
}
-/* right now, we do things simple and dirty (read: our current preecho
- is a joke). Should this prove inadequate, then we'll think of
- something different. The details of the encoding format do not
- depend on the exact behavior, only the format of the bits that come
- out.
-
- Mark Taylor probably has much witter ways of doing this... Let's
- see if simple delta analysis gives us acceptible results for now. */
+/* use MDCT for spectral power estimation */
static void _ve_deltas(double *deltas,double *pcm,int n,double *win,
int winsize){
- int i,j,p=winsize/2;
+ int i,j;
+ mdct_lookup m;
+ double out[winsize/2];
+
+ mdct_init(&m,winsize);
+
for(j=0;j<n;j++){
- p-=winsize/2;
- for(i=0;i<winsize-1;i++,p++){
- double temp=fabs(win[i]*pcm[p]-win[i+1]*pcm[p+1]);
- if(deltas[j]<temp)deltas[j]=temp;
- }
- p++;
+ double acc=0.;
+
+ mdct_forward(&m,pcm+j*winsize/2,out,win);
+ for(i=0;i<winsize/2;i++)
+ acc+=fabs(out[i]);
+ if(deltas[j]<acc)deltas[j]=acc;
}
+
+ mdct_clear(&m);
+
+#ifdef ANALYSIS
+ {
+ static int frameno=0;
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"delta%d.m",frameno++);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n;j++)
+ fprintf(out,"%g\n",deltas[j]);
+ fclose(out);
+ }
+#endif
}
void _ve_envelope_multipliers(vorbis_dsp_state *v){
double clamp;
int i;
+#ifdef ANALYSIS
+ {
+ static int frameno=0.;
+ FILE *out;
+ int j;
+ char buffer[80];
+
+ sprintf(buffer,"del%d.m",frameno++);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n;j++)
+ fprintf(out,"%g\n",mult[j]);
+ fclose(out);
+ }
+#endif
+
/* are we going to multiply anything? */
for(i=1;i<n;i++){
for(i=0;i<begin;i++)mult[i]=0;
- last=1;
for(;i<n;i++){
- if(mult[i]/last>clamp*vi->preecho_thresh){
- last=mult[i]/vi->preecho_clamp;
+ if(mult[i]>=last*vi->preecho_thresh){
+ last=mult[i];
- mult[i]=floor(log(mult[i]/clamp/vi->preecho_clamp)/log(2))-1;
+ mult[i]=floor(log(mult[i]/clamp)/log(2));
if(mult[i]>15)mult[i]=15;
if(mult[i]<1)mult[i]=1;
}
}else
memset(mult,0,sizeof(double)*n);
+
+#ifdef ANALYSIS
+ {
+ static int frameno=0.;
+ FILE *out;
+ int j;
+ char buffer[80];
+
+ sprintf(buffer,"sparse%d.m",frameno++);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n;j++)
+ fprintf(out,"%g\n",mult[j]);
+ fclose(out);
+ }
+#endif
+
+
}
}
void _ve_envelope_apply(vorbis_block *vb,int multp){
+ static int frameno=0;
vorbis_info *vi=vb->vd->vi;
double env[vb->multend*vi->envelopesa];
envelope_lookup *look=&vb->vd->ve;
- int i,j;
+ int i,j,k;
for(i=0;i<vi->envelopech;i++){
double *mult=vb->mult[i];
/* compute the envelope curve */
if(_ve_envelope_generate(mult,env,look->window,vb->multend,
vi->envelopesa)){
- if(multp)
+#ifdef ANALYSIS
+ {
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"env%d.m",frameno);
+ out=fopen(buffer,"w+");
for(j=0;j<vb->multend*vi->envelopesa;j++)
- vb->pcm[i][j]*=env[j];
- else
+ fprintf(out,"%g\n",env[j]);
+ fclose(out);
+ sprintf(buffer,"prepcm%d.m",frameno);
+ out=fopen(buffer,"w+");
for(j=0;j<vb->multend*vi->envelopesa;j++)
- vb->pcm[i][j]/=env[j];
+ fprintf(out,"%g\n",vb->pcm[i][j]);
+ fclose(out);
+ }
+#endif
+ for(k=0;k<vi->channels;k++){
+
+ if(multp)
+ for(j=0;j<vb->multend*vi->envelopesa;j++)
+ vb->pcm[k][j]*=env[j];
+ else
+ for(j=0;j<vb->multend*vi->envelopesa;j++)
+ vb->pcm[k][j]/=env[j];
+
+ }
+
+#ifdef ANALYSIS
+ {
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"pcm%d.m",frameno);
+ out=fopen(buffer,"w+");
+ for(j=0;j<vb->multend*vi->envelopesa;j++)
+ fprintf(out,"%g\n",vb->pcm[i][j]);
+ fclose(out);
+ }
+#endif
}
+ frameno++;
}
}
/* handle the flat settings first */
memcpy(vi,&(predef_modes[mode]),sizeof(vorbis_info));
+ vi->threshhold_points=threshhold_points;
vi->user_comments=calloc(1,sizeof(char *));
- vi->vendor=strdup("Xiphophorus libVorbis I 19991012");
+ vi->vendor=strdup("Xiphophorus libVorbis I 19991018");
return(0);
}
int n=l->ln;
int m=l->m;
double aut[m+1],work[n+n],error;
- double fscale=1./n;
+ double fscale=.5/n;
int i,j;
/* input is a real curve. make it complex-real */
for(i=0;i<n;i++){
/* how much 'real estate' in the log domain does the bin in the
linear domain represent? */
- double logA=LOG_X(i-.5,bias);
- double logB=LOG_X(i+.5,bias);
+ double logA=LOG_X(i,bias);
+ double logB=LOG_X(i+1.,bias);
l->norm[i]=logB-logA; /* this much */
}
/* decode; encode may use this too */
drft_init(&l->fft,mapped*2);
- {
- double w=1./oct*M_PI;
- for(i=0;i<n;i++){
- l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
- if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
- }
+ for(i=0;i<n;i++){
+ l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
+ if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
}
}
optimization; it could both be more precise as well as not compute
quite a few unused values */
-static void _vlpc_de_helper(double *curve,double *lpc,double amp,
+void _vlpc_de_helper(double *curve,double *lpc,double amp,
lpc_lookup *l){
int i;
memset(curve,0,sizeof(double)*l->ln*2);
-
+ if(amp==0)return;
+
for(i=0;i<l->m;i++){
curve[i*2+1]=lpc[i]/4/amp;
curve[i*2+2]=-lpc[i]/4/amp;
{
int l2=l->ln*2;
double unit=1./amp;
- curve[0]=(1./(curve[0]+unit));
+ curve[0]=(1./(curve[0]*2+unit));
for(i=1;i<l->ln;i++){
double real=(curve[i]+curve[l2-i]);
double imag=(curve[i]-curve[l2-i]);
int i;
_vlpc_de_helper(lcurve,lpc,amp,l);
+ if(amp==0)return;
for(i=0;i<l->n;i++)
curve[i]=lcurve[l->iscale[i]]*l->norm[i];
double lcurve[l->ln*2];
int i;
- _vlpc_de_helper(lcurve,lpc,amp,l);
+ if(amp==0){
+ memset(residue,0,l->n*sizeof(double));
+ }else{
+
+ _vlpc_de_helper(lcurve,lpc,amp,l);
- for(i=0;i<l->n;i++)
- residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
+ for(i=0;i<l->n;i++)
+ residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
+ }
}
+
function: predefined encoding modes
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 2 1999
+ last modification date: Oct 15 1999
********************************************************************/
#include <stdio.h>
#include "codec.h"
+double threshhold_points[THRESH_POINTS]=
+/* 0Hz 24kHz
+ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 */
+{0.,.01,.02,.03,.04,.06,.08,.1,.15,.2,.25,.3,.34,.4,.45,.5,.6,.7,.8,1.};
+
vorbis_info predef_modes[]={
/* CD quality stereo, no channel coupling */
- { 2, 44100, 0, NULL, 0, NULL, {512, 2048}, {30,30}, {5,5}, 64,2,2, 10, .5,
- NULL,NULL,NULL},
+ /* channels, sample rate, dummy, dummy, dummy, dummy */
+ { 2, 44100, 0, NULL, 0, NULL,
+ /* smallblock, largeblock, LPC order (small, large) */
+ {512, 2048}, {16,16},
+ /* spectral octaves (small, large), envelope segment size */
+ {5,5}, 64,
+ /* envelope channels, spectrum channels */
+ 1,2,
+ /* preecho clamp trigger threshhold, clamp range, dummy */
+ 4, 2, NULL,
+ /* noise masking curve dB attenuation levels [20] */
+ {-12,-12,-18,-18,-18,-18,-18,-10,-5,-2,
+ 0,0,0,0,1,2,3,3,4,5},
+ /* noise masking scale biases */
+ .95,1.01,.001,
+ /* tone masking curve dB attenuation levels [20] */
+ {-24,-24,-24,-24,-24,-24,-24,-24,-24,-24,
+ -24,-24,-24,-24,-24,-24,-24,-24,-24,-24},
+ /* tone masking rolloff settings (dB per octave), octave bias */
+ 90,60,.001,
+ NULL,NULL,NULL},
+
};
#define predef_mode_max 0
#endif
-
function: random psychoacoustics (not including preecho)
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Aug 26 1999
+ last modification date: Oct 15 1999
********************************************************************/
+#include <stdlib.h>
#include <math.h>
#include <string.h>
-#include "stdio.h"
+#include <stdio.h>
#include "codec.h"
#include "psy.h"
#include "lpc.h"
#include "smallft.h"
#include "xlogmap.h"
-#define NOISEdB -6
+/* Set up decibel threshhold 'curves'. Actually, just set a level at
+ log frequency intervals, interpolate, and call it a curve. */
+
+static void set_curve(double *points,
+ double *ref,int rn,double rrate,
+ double *c,int n, double crate){
+ int i,j=0;
+
+ for(i=0;i<rn-1;i++){
+ int endpos=points[i+1]*n*rrate/crate;
+ double base=ref[i];
+ double delta=(ref[i+1]-base)/(endpos-j);
+ for(;j<endpos && j<n;j++){
+ c[j]=base;
+ base+=delta;
+ }
+ }
+}
-#define MASKdB 20
-#define HROLL 60
-#define LROLL 90
-#define MASKBIAS 10
+void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n){
+ memset(p,0,sizeof(psy_lookup));
+ p->noisethresh=malloc(n*sizeof(double));
+ p->maskthresh=malloc(n*sizeof(double));
+ p->vi=vi;
+ p->n=n;
-#define LNOISE .95
-#define HNOISE 1.01
-#define NOISEBIAS 20
+ /* set up the curves for a given blocksize and sample rate */
+
+ set_curve(vi->threshhold_points,
+ vi->noisethresh,THRESH_POINTS,48000,p->noisethresh,n,vi->rate);
+ set_curve(vi->threshhold_points,
+ vi->maskthresh, THRESH_POINTS,48000,p->maskthresh, n,vi->rate);
+
+#ifdef ANALYSIS
+ {
+ int j;
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"noise_threshhold_%d.m",n);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n;j++)
+ fprintf(out,"%g\n",p->noisethresh[j]);
+ fclose(out);
+ sprintf(buffer,"mask_threshhold_%d.m",n);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n;j++)
+ fprintf(out,"%g\n",p->maskthresh[j]);
+ fclose(out);
+ }
+#endif
+
+}
+
+void _vp_psy_clear(psy_lookup *p){
+ if(p){
+ if(p->noisethresh)free(p->noisethresh);
+ if(p->maskthresh)free(p->maskthresh);
+ memset(p,0,sizeof(psy_lookup));
+ }
+}
/* Find the mean log energy of a given 'band'; used to evaluate tones
against background noise */
scale we seek, and gives OK results. So, that means it's a good
hack */
-/* To add: f scale noise attenuation curve */
-
-void _vp_noise_floor(double *f, double *m,int n){
+void _vp_noise_floor(psy_lookup *p, double *f, double *m){
+ int n=p->n;
+ vorbis_info *vi=p->vi;
+
long lo=0,hi=0;
- double acc=0,div=0;
+ double acc=0.;
+ double div=0.;
int i,j;
+
+ double bias=n*vi->noisebias;
- for(i=100;i<n;i++){
- long newlo=i*LNOISE-NOISEBIAS;
- long newhi=i*HNOISE+NOISEBIAS;
+ for(i=0;i<n;i++){
+ long newlo=i*vi->lnoise-bias;
+ long newhi=i*vi->hnoise+bias;
double temp;
if(newhi>n)newhi=n;
if(newlo<0)newlo=0;
-
+
for(j=hi;j<newhi;j++){
acc+=todB(f[j]);
- div++;
+ div+=1.;
}
for(j=lo;j<newlo;j++){
acc-=todB(f[j]);
- div--;
+ div-=1.;
}
hi=newhi;
lo=newlo;
- temp=fromdB(acc/div+NOISEdB); /* The NOISEdB constant should be an
- attenuation curve */
+ /* attenuate by the noise threshhold curve */
+ temp=fromdB(acc/div+p->noisethresh[i]);
if(m[i]<temp)m[i]=temp;
}
}
-/* figure the masking curve. linear rolloff on a dB scale, adjusted
- by octave */
-void _vp_mask_floor(double *f, double *m,int n){
+/* Masking curve: linear rolloff on a dB scale, adjusted by octave,
+ attenuated by maskthresh */
+
+void _vp_mask_floor(psy_lookup *p,double *f, double *m){
+ int n=p->n;
+ double hroll=p->vi->hroll;
+ double lroll=p->vi->lroll;
double ocSCALE=1./log(2);
double curmask=-9.e40;
- double curoc=log(MASKBIAS)*ocSCALE;
+ double maskbias=n*p->vi->maskbias;
+ double curoc=log(maskbias)*ocSCALE;
long i;
/* run mask forward then backward */
for(i=0;i<n;i++){
- double newmask=todB(f[i])-MASKdB;
- double newoc=log(i+MASKBIAS)*ocSCALE;
- double roll=curmask-(newoc-curoc)*HROLL;
- double lroll;
+ double newmask=todB(f[i])+p->maskthresh[i];
+ double newoc=log(i+maskbias)*ocSCALE;
+ double roll=curmask-(newoc-curoc)*hroll;
+ double troll;
if(newmask>roll){
roll=curmask=newmask;
curoc=newoc;
}
- lroll=fromdB(roll);
- if(m[i]<lroll)m[i]=lroll;
+ troll=fromdB(roll);
+ if(m[i]<troll)m[i]=troll;
}
curmask=-9.e40;
- curoc=log(n+MASKBIAS)*ocSCALE;
+ curoc=log(n+maskbias)*ocSCALE;
for(i=n-1;i>=0;i--){
- double newmask=todB(f[i])-MASKdB;
- double newoc=log(i+MASKBIAS)*ocSCALE;
- double roll=curmask-(curoc-newoc)*LROLL;
- double lroll;
+ double newmask=todB(f[i])+p->maskthresh[i];
+ double newoc=log(i+maskbias)*ocSCALE;
+ double roll=curmask-(curoc-newoc)*lroll;
+ double troll;
if(newmask>roll){
roll=curmask=newmask;
curoc=newoc;
}
- lroll=fromdB(roll);
- if(m[i]<lroll)m[i]=lroll;
+ troll=fromdB(roll);
+ if(m[i]<troll)m[i]=troll;
}
}
#ifndef _V_PSY_H_
#define _V_PSY_H_
+extern void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n);
+extern void _vp_psy_clear(psy_lookup *p);
+
+extern void _vp_noise_floor(psy_lookup *p, double *f, double *m);
+extern void _vp_mask_floor(psy_lookup *p,double *f, double *m);
+
+
extern double _vp_balance_compute(double *A, double *B, double *lpc,
lpc_lookup *vb);
extern void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
function: spectrum envelope and residue code/decode
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 10 1999
+ last modification date: Oct 17 1999
********************************************************************/
/* this code is still seriously abbreviated. I'm filling in pieces as
we go... --Monty 19991004 */
+/* unlike other LPC-based coders, we never apply the filter, only
+ inspect the frequency response, thus we don't need to guard against
+ instability. However, two coefficients quantising to the same
+ value will cause the response to explode. */
+
int _vs_spectrum_encode(vorbis_block *vb,double amp,double *lsp){
/* no real coding yet. Just write out full sized words for now
because people need bitstreams to work with */
int scale=vb->W;
+ int m=vb->vd->vi->floororder[scale];
int n=vb->pcmend/2;
int last=0;
+ double dlast=0.;
+ double min=M_PI/n/2.;
+
int bits=rint(log(n)/log(2));
int i;
-
_oggpack_write(&vb->opb,amp*327680,18);
- for(i=0;i<vb->vd->vi->floororder[scale];i++){
+ for(i=0;i<m;i++){
int val=rint(lsp[i]/M_PI*n-last);
- lsp[i]=(last+=val)*M_PI/n;
_oggpack_write(&vb->opb,val,bits);
+ lsp[i]=(last+=val)*M_PI/n;
+
+ /* Underpowered but sufficient */
+ if(lsp[i]<dlast+min)lsp[i]=dlast+min;
+ dlast=lsp[i];
}
return(0);
}
int _vs_spectrum_decode(vorbis_block *vb,double *amp,double *lsp){
int scale=vb->W;
+ int m=vb->vd->vi->floororder[scale];
int n=vb->pcmend/2;
int last=0;
+ double dlast=0.;
int bits=rint(log(n)/log(2));
int i;
+ double min=M_PI/n/2.;
*amp=_oggpack_read(&vb->opb,18)/327680.;
- for(i=0;i<vb->vd->vi->floororder[scale];i++){
+ for(i=0;i<m;i++){
int val=_oggpack_read(&vb->opb,bits);
lsp[i]=(last+=val)*M_PI/n;
+ /* Underpowered but sufficient */
+ if(lsp[i]<dlast+min)lsp[i]=dlast+min;
+ dlast=lsp[i];
}
return(0);
}
if(val>16)val=16;
if(val<-16)val=-16;
- /*if(val==0){
+ if(val==0 || val==2 || val==-2){
if(data[i]<0){
val=-1;
}else{
val=1;
}
- }*/
+ }
data[i]=val;
/*if(val<0){
extern void _vs_residue_quantize(double *data,double *curve,
vorbis_info *vi,int n);
extern int _vs_residue_encode(vorbis_block *vb,double *data);
+extern int _vs_residue_decode(vorbis_block *vb,double *data);
#endif
#include "spectrum.h"
int vorbis_synthesis(vorbis_block *vb,ogg_packet *op){
- static int frameno=0;
vorbis_dsp_state *vd=vb->vd;
vorbis_info *vi=vd->vi;
oggpack_buffer *opb=&vb->opb;