Commit a8e19c3c authored by Monty's avatar Monty
Browse files

Keeping up to date.

Monty

svn path=/trunk/vorbis/; revision=24
parent 1736dd22
# vorbis makefile configured for use with gcc on any platform
# $Id: Makefile.in,v 1.7 1999/08/07 23:36:32 xiphmont Exp $
# $Id: Makefile.in,v 1.8 1999/08/09 21:11:38 xiphmont Exp $
###############################################################################
# #
......@@ -27,7 +27,8 @@ AR=@AR@
RANLIB=@RANLIB@
LIBS=@LIBS@ -lm
OFILES = framing.o mdct.o smallft.o block.o envelope.o window.o
OFILES = framing.o mdct.o smallft.o block.o envelope.o window.o\
lsp.o lpc.o analysis.o synthesis.o psy.o
TARGETFILES = libvorbis.a
all:
......@@ -48,9 +49,8 @@ selftest:
@./test_framing
@./test_bitwise
dsptest:
$(MAKE) target CFLAGS="$(DEBUG)"
$(CC) $(DEBUG) $(LDFLAGS) $(OFILES) dsptest.c -o dsptest -lm
dsptest: debug
$(CC) $(DEBUG) $(LDFLAGS) dsptest.c libvorbis.a -o dsptest -lm
target: $(TARGETFILES)
......@@ -64,4 +64,4 @@ $(OFILES): mdct.h
$(CC) $(CFLAGS) -c $<
clean:
-rm -f *.o *.a test* *~ *.out ogg config.*
-rm -f *.o *.a test* *~ *.out ogg config.* dsptest
......@@ -14,12 +14,126 @@
function: single-block PCM analysis
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Jul 28 1999
last modification date: Aug 09 1999
********************************************************************/
analysis_packetout(vorbis_dsp_state *v, vorbis_block *vb,
ogg_packet *op){
#include <stdio.h>
#include <math.h>
#include "lpc.h"
#include "codec.h"
#include "envelope.h"
#include "mdct.h"
#include "psy.h"
int vorbis_analysis(vorbis_block *vb){
int i;
double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
int n=vb->pcmend;
static int frameno=0;
_ve_envelope_sparsify(vb);
_ve_envelope_apply(vb,0);
for(i=0;i<vb->pcm_channels;i++)
mdct_forward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
/* handle vectors to be encoded one at a time. The PCM storage now
has twice the space we need in the MDCT domain, so mix down into
the second half of storage */
/* XXXX Mix down to C+D stereo now */
{
double *C=vb->pcm[0]+n/2;
double *D=vb->pcm[1]+n/2;
double *L=vb->pcm[0];
double *R=vb->pcm[1];
for(i=0;i<n/2;i++){
C[i]=(L[i]+R[i])*.5;
D[i]=(L[i]-R[i])*.5;
}
/* extract the spectral envelope and residue */
{
double floor1[n/2];
double floor2[n/2];
double curve[n/2];
double work[n/2];
double lpc1[40];
double lpc2[40];
double amp1,amp2;
memset(floor1,0,sizeof(floor1));
memset(floor2,0,sizeof(floor2));
memset(work,0,sizeof(work));
_vp_noise_floor(C,floor1,n/2);
/* _vp_noise_floor(D,floor1,n/2);*/
_vp_mask_floor(C,floor1,n/2);
/* _vp_mask_floor(D,floor1,n/2);*/
_vp_psy_sparsify(C,floor1,n/2);
_vp_psy_sparsify(D,floor1,n/2);
memmove(floor1+n/4,floor1,(n/4)*sizeof(double));
amp2=sqrt(vorbis_curve_to_lpc(floor1,n/2,lpc2,20));
vorbis_lpc_to_curve(work,n/2,lpc2,amp2,20);
/* amp2=sqrt(vorbis_curve_to_lpc(floor1,n/2/8,lpc2,12));
vorbis_lpc_to_curve(work,n/2/8,lpc2,amp2,12);
amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/8,n/2/8,lpc2,10));
vorbis_lpc_to_curve(work+n/2/8,n/2/8,lpc2,amp2,10);
amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/4,n/2/4,lpc2,6));
vorbis_lpc_to_curve(work+n/2/4,n/2/4,lpc2,amp2,6);
amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/2,n/2/2,lpc2,6));
vorbis_lpc_to_curve(work+n/2/2,n/2/2,lpc2,amp2,6);*/
_vp_psy_quantize(C,work,n/2);
_vp_psy_quantize(D,work,n/2);
{
FILE *out;
char buffer[80];
sprintf(buffer,"mdctC%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++)
fprintf(out,"%g\n",fabs(C[i]));
fclose(out);
sprintf(buffer,"mdctD%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++)
fprintf(out,"%g\n",fabs(D[i]));
fclose(out);
sprintf(buffer,"floor%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++)
fprintf(out,"%g\n",floor1[i]);
fclose(out);
sprintf(buffer,"lpc%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++)
fprintf(out,"%g\n",work[i]);
fclose(out);
frameno++;
}
}
}
return(0);
}
int vorbis_analysis_packetout(vorbis_block *vb, ogg_packet *op){
/* find block's envelope vector and apply it */
......
......@@ -227,10 +227,13 @@ int vorbis_analysis_wrote(vorbis_dsp_state *v, int vals){
int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb){
int i;
memset(vb,0,sizeof(vorbis_block));
vb->pcm_storage=v->block_size[1];
vb->pcm_channels=v->pcm_channels;
vb->mult_storage=v->block_size[1]/v->samples_per_envelope_step;
vb->mult_channels=v->envelope_channels;
vb->floor_channels=v->vi.floorch;
vb->floor_storage=v->vi.floororder;
vb->pcm=malloc(vb->pcm_channels*sizeof(double *));
for(i=0;i<vb->pcm_channels;i++)
......@@ -239,6 +242,15 @@ int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb){
vb->mult=malloc(vb->mult_channels*sizeof(double *));
for(i=0;i<vb->mult_channels;i++)
vb->mult[i]=malloc(vb->mult_storage*sizeof(double));
vb->lsp=malloc(vb->floor_channels*sizeof(double *));
vb->lpc=malloc(vb->floor_channels*sizeof(double *));
vb->amp=malloc(vb->floor_channels*sizeof(double));
for(i=0;i<vb->floor_channels;i++){
vb->lsp[i]=malloc(vb->floor_storage*sizeof(double));
vb->lpc[i]=malloc(vb->floor_storage*sizeof(double));
}
return(0);
}
......@@ -284,42 +296,46 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
/* overconserve for now; any block with a non-placeholder multiplier
should be minimal size. We can be greedy and only look at nW size */
if(v->W)
/* this is a long window; we start the search forward of centerW
because that's the fastest we could react anyway */
i=v->centerW+v->block_size[1]/4-v->block_size[0]/4;
else
/* short window. Search from centerW */
i=v->centerW;
i/=v->samples_per_envelope_step;
for(;i<v->envelope_current;i++){
for(j=0;j<v->envelope_channels;j++)
if(v->multipliers[j][i-1]*v->vi.preecho_thresh<
v->multipliers[j][i])break;
if(j<v->envelope_channels)break;
}
if(i<v->envelope_current){
/* Ooo, we hit a multiplier. Is it beyond the boundary to make the
upcoming block large ? */
long largebound;
if(v->vi.smallblock<v->vi.largeblock){
if(v->W)
largebound=v->centerW+v->block_size[1];
/* this is a long window; we start the search forward of centerW
because that's the fastest we could react anyway */
i=v->centerW+v->block_size[1]/4-v->block_size[0]/4;
else
largebound=v->centerW+v->block_size[0]/4+v->block_size[1]*3/4;
largebound/=v->samples_per_envelope_step;
if(i>=largebound)
/* short window. Search from centerW */
i=v->centerW;
i/=v->samples_per_envelope_step;
for(;i<v->envelope_current;i++){
for(j=0;j<v->envelope_channels;j++)
if(v->multipliers[j][i-1]*v->vi.preecho_thresh<
v->multipliers[j][i])break;
if(j<v->envelope_channels)break;
}
if(i<v->envelope_current){
/* Ooo, we hit a multiplier. Is it beyond the boundary to make the
upcoming block large ? */
long largebound;
if(v->W)
largebound=v->centerW+v->block_size[1];
else
largebound=v->centerW+v->block_size[0]/4+v->block_size[1]*3/4;
largebound/=v->samples_per_envelope_step;
if(i>=largebound)
v->nW=1;
else
v->nW=0;
}else{
/* Assume maximum; if the block is incomplete given current
buffered data, this will be detected below */
v->nW=1;
else
v->nW=0;
}else{
/* Assume maximum; if the block is incomplete given current
buffered data, this will be detected below */
}
}else
v->nW=1;
}
/* Do we actually have enough data *now* for the next block? The
reason to check is that if we had no multipliers, that could
......@@ -343,7 +359,9 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
vb->pcmend=v->block_size[v->W];
vb->multend=vb->pcmend / v->samples_per_envelope_step;
if(v->pcm_channels!=vb->pcm_channels ||
if(vb->floor_channels!=v->vi.floorch ||
vb->floor_storage!=v->vi.floororder ||
v->pcm_channels!=vb->pcm_channels ||
v->block_size[1]!=vb->pcm_storage ||
v->envelope_channels!=vb->mult_channels){
......
......@@ -53,7 +53,11 @@ typedef struct vorbis_info{
int envelopesa;
int envelopech;
int *envelopemap; /* which envelope applies to what pcm channel */
int **channelmap; /* which encoding channels produce what pcm */
int **Echannelmap; /* which encoding channels produce what pcm (decode) */
int **channelmap; /* which pcm channels produce what floors (encode) */
int floororder;
int floorch;
int *floormap;
double preecho_thresh;
double preecho_clamp;
......@@ -154,10 +158,16 @@ typedef struct vorbis_dsp_state{
typedef struct vorbis_block{
double **pcm;
double **mult;
double **lpc;
double **lsp;
double *amp;
int pcm_channels; /* allocated, not used */
int pcm_storage; /* allocated, not used */
int mult_channels; /* allocated, not used */
int mult_storage; /* allocated, not used */
int floor_channels;
int floor_storage;
long lW;
long W;
......@@ -225,8 +235,8 @@ extern double **vorbis_analysis_buffer(vorbis_dsp_state *vd,int vals);
extern int vorbis_analysis_wrote(vorbis_dsp_state *vd,int vals);
extern int vorbis_analysis_blockout(vorbis_dsp_state *vd,
vorbis_block *vb);
extern int vorbis_analysis_packetout(vorbis_dsp_state *vd,
vorbis_block *vb,
extern int vorbis_analysis(vorbis_block *vb);
extern int vorbis_analysis_packetout(vorbis_block *vb,
ogg_packet *op);
/* Vorbis PRIMITIVES: synthesis layer *******************************/
......@@ -234,6 +244,7 @@ extern int vorbis_analysis_packetout(vorbis_dsp_state *vd,
extern void vorbis_synthesis_free(vorbis_dsp_state *vd);
extern int vorbis_synthesis_init(vorbis_dsp_state *vd,vorbis_info *vi);
extern int vorbis_synthesis(vorbis_block *vb);
extern int vorbis_synthesis_blockin(vorbis_dsp_state *v,vorbis_block *vb);
extern int vorbis_synthesis_pcmout(vorbis_dsp_state *v,double ***pcm);
extern int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes);
......@@ -242,5 +253,10 @@ extern int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes);
extern int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb);
#define min(x,y) ((x)>(y)?(y):(x))
#define max(x,y) ((x)<(y)?(y):(x))
#define todB(x) (x==0?-9.e40:log(fabs(x))*8.6858896) /* 20log10(x) */
#define fromdB(x) (exp((x)*.11512925))
#endif
#include <math.h>
#include <stdio.h>
#include "codec.h"
#include "envelope.h"
#include "mdct.h"
#define READ 4096
#define READ 1024
int main(){
vorbis_dsp_state encode,decode;
......@@ -15,11 +13,11 @@ int main(){
int done=0;
char *temp[]={ "Test" ,"the Test band", "test records",NULL };
int mtemp[]={0,1};
int mtemp2[]={0,1};
int frame=0;
signed char buffer[READ*4+44];
vi.channels=2;
vi.rate=44100;
vi.version=0;
......@@ -31,9 +29,12 @@ int main(){
vi.envelopesa=64;
vi.envelopech=2;
vi.envelopemap=mtemp;
vi.floormap=mtemp2;
vi.floororder=20;
vi.floorch=2;
vi.channelmap=NULL;
vi.preecho_thresh=10.;
vi.preecho_clamp=4.;
vi.preecho_clamp=.5;
vorbis_analysis_init(&encode,&vi);
vorbis_synthesis_init(&decode,&vi);
......@@ -60,45 +61,14 @@ int main(){
while(vorbis_analysis_blockout(&encode,&vb)){
double **pcm;
int avail;
double *window=encode.window[vb.W][vb.lW][vb.nW];
/* analysis */
/* apply envelope */
_ve_envelope_sparsify(&vb);
_ve_envelope_apply(&vb,0);
for(i=0;i<vb.pcm_channels;i++)
mdct_forward(&vb.vd->vm[vb.W],vb.pcm[i],vb.pcm[i],window);
vorbis_analysis(&vb);
/* synthesis */
for(i=0;i<vb.pcm_channels;i++)
mdct_backward(&vb.vd->vm[vb.W],vb.pcm[i],vb.pcm[i],window);
/*{
FILE *out;
char path[80];
int i;
int avail=encode.block_size[vb.W];
int beginW=countermid-avail/2;
sprintf(path,"ana%d",vb.frameno);
out=fopen(path,"w");
for(i=0;i<avail;i++)
fprintf(out,"%d %g\n",i+beginW,vb.pcm[0][i]);
fprintf(out,"\n");
for(i=0;i<avail;i++)
fprintf(out,"%d %g\n",i+beginW,window[i]);
fclose(out);
}*/
_ve_envelope_apply(&vb,1);
vorbis_synthesis(&vb);
counterin+=bread/4;
vorbis_synthesis_blockin(&decode,&vb);
......
......@@ -53,7 +53,7 @@ Carsten Bormann
#include "smallft.h"
#include "lpc.h"
/* This is pared down for Vorbis, where we only use LPC to encode
/* This is pared down for Vorbis where we only use LPC to encode
spectral envelope curves. Thus we only are interested in
generating the coefficients and recovering the curve from the
coefficients. Autocorrelation LPC coeff generation algorithm
......@@ -64,6 +64,7 @@ Carsten Bormann
double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
double aut[m+1],work[n+n],error;
drft_lookup dl;
int i,j;
/* input is a real curve. make it complex-real */
......@@ -73,7 +74,9 @@ double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
}
n*=2;
fft_backward(n,work,NULL,NULL);
drft_init(&dl,n);
drft_backward(&dl,work);
drft_clear(&dl);
/* The autocorrelation will not be circular. Shift, else we lose
most of the power in the edges. */
......@@ -142,7 +145,7 @@ static double vorbis_lpc_magnitude(double w,double *lpc, int m){
real+=lpc[k]*cos((k+1)*w);
imag+=lpc[k]*sin((k+1)*w);
}
return(1/sqrt(real*real+imag*imag));
return(1./sqrt(real*real+imag*imag));
}
/* generate the whole freq response curve on an LPC IIR filter */
......
......@@ -145,13 +145,17 @@ static inline void _mdct_kernel(double *x,
for(r=0;r<(n4>>i);r+=4){
int w1=wbase;
int w2=wbase-(k0>>1);
double AEv= *AE,wA;
double AOv= *AO,wB;
wbase-=4;
for(s=0;s<(2<<i);s++){
x[w1+2]=w[w1+2]+w[w2+2];
x[w1] =w[w1]+w[w2];
x[w2+2]=(w[w1+2]-w[w2+2])* *AE-(w[w1]-w[w2])* *AO;
x[w2] =(w[w1]-w[w2])* *AE+(w[w1+2]-w[w2+2])* *AO;
wA =w[w1+2]-w[w2+2];
wB =w[w1]-w[w2];
x[w2+2]=wA*AEv - wB*AOv;
x[w2] =wB*AEv + wA*AOv;
w1-=k0;
w2-=k0;
}
......@@ -317,16 +321,16 @@ void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
double scale=n/4.;
for(i=0;i<n4;i++){
double BEv=BE[i];
double BOv=BO[i];
#ifdef VORBIS_SPECIFIC_MODIFICATIONS
double temp1= (*x * *BO - *(x+2) * *BE)/ scale;
double temp2= (*x * *BE + *(x+2) * *BO)/ -scale;
out[o1]=-temp1*window[o1];
out[o2]=temp1*window[o2];
out[o3]=out[o4]=temp2;
double temp= (*x * BOv - *(x+2) * BEv)/ scale;
out[o3]=out[o4]=(*x * BEv + *(x+2) * BOv)/ -scale;
out[o1]=-temp*window[o1];
out[o2]=temp*window[o2];
#else
double temp1= (*x * *BO - *(x+2) * *BE)* scale;
double temp2= (*x * *BE + *(x+2) * *BO)* -scale;
double temp1= (*x * BOv - *(x+2) * BEv)* scale;
double temp2= (*x * BEv + *(x+2) * BOv)* -scale;
out[o1]=-temp1*window[o1];
out[o2]=temp1*window[o2];
......@@ -339,7 +343,6 @@ void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
o3++;
o4--;
x+=4;
BE++;BO++;
}
}
}
......
/********************************************************************
* *
* 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: random psychoacoustics (not including preecho)
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Aug 08 1999
********************************************************************/
#include <math.h>
#include "codec.h"
#include "psy.h"
#define NOISEdB 6
#define MASKdB 12
#define HROLL 60
#define LROLL 90
#define MASKBIAS 40
#define LNOISE .8
#define HNOISE 1.01
#define NOISEBIAS 20
/* Find the mean log energy of a given 'band'; used to evaluate tones
against background noise */
/* This is faster than a real convolution, gives us roughly the log f
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){
long lo=0,hi=0;
double acc=0,div=0;
int i,j;
for(i=0;i<n;i++){
long newlo=i*LNOISE-NOISEBIAS;
long newhi=i*HNOISE+NOISEBIAS;
double temp;
if(newhi>n)newhi=n;
if(newlo<0)newlo=0;
for(j=hi;j<newhi;j++){
acc+=todB(f[j]);
div++;
}
for(j=lo;j<newlo;j++){
acc-=todB(f[j]);
div--;
}
hi=newhi;
lo=newlo;
temp=fromdB(acc/div+NOISEdB); /* The NOISEdB constant should be an
attenuation curve */
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){
double ocSCALE=1./log(2);
double curmask=-9.e40;
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;
if(newmask>roll){
roll=curmask=newmask;
curoc=newoc;
}
lroll=fromdB(roll);
if(m[i]<lroll)m[i]=lroll;
}
curmask=-9.e40;
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;
if(newmask>roll){
roll=curmask=newmask;