Commit f5d723f5 authored by Monty's avatar Monty
Browse files

Fixes to dsp

Monty 19990807

svn path=/trunk/vorbis/; revision=21
parent ea0d1fec
# vorbis makefile configured for use with gcc on any platform
# $Id: Makefile.in,v 1.6 1999/08/07 00:09:13 xiphmont Exp $
# $Id: Makefile.in,v 1.7 1999/08/07 23:36:32 xiphmont Exp $
###############################################################################
# #
......@@ -44,14 +44,13 @@ selftest:
$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST framing.c -o test_framing
$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST bitwise.c\
-o test_bitwise -lm
$(CC) $(DEBUG) -c envelope.c mdct.c window.c
$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST envelope.o mdct.o window.c\
block.c -o test_blocking -lm
@echo
@./test_framing
@./test_bitwise
@./test_blocking
dsptest:
$(MAKE) target CFLAGS="$(DEBUG)"
$(CC) $(DEBUG) $(LDFLAGS) $(OFILES) dsptest.c -o dsptest -lm
target: $(TARGETFILES)
......
......@@ -75,6 +75,8 @@ static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
memcpy(&v->vi,vi,sizeof(vorbis_info));
_ve_envelope_init(&v->ve,vi->envelopesa);
mdct_init(&v->vm[0],vi->smallblock);
mdct_init(&v->vm[1],vi->largeblock);
v->samples_per_envelope_step=vi->envelopesa;
v->block_size[0]=vi->smallblock;
......@@ -139,17 +141,9 @@ static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
/* arbitrary settings and spec-mandated numbers get filled in here */
int vorbis_analysis_init(vorbis_dsp_state *v,vorbis_info *vi){
vi->smallblock=512;
vi->largeblock=2048;
vi->envelopesa=64;
vi->envelopech=vi->channels;
vi->preecho_thresh=10.;
vi->preecho_thresh=4.;
vi->envelopemap=calloc(2,sizeof(int));
vi->envelopemap[0]=0;
vi->envelopemap[1]=1;
_vds_shared_init(v,vi);
drft_init(&v->vf[0],vi->smallblock); /* only used in analysis */
drft_init(&v->vf[1],vi->largeblock); /* only used in analysis */
return(0);
}
......@@ -173,6 +167,11 @@ void vorbis_analysis_clear(vorbis_dsp_state *v){
if(v->multipliers[i])free(v->multipliers[i]);
free(v->multipliers);
}
_ve_envelope_clear(&v->ve);
drft_clear(&v->vf[0]);
drft_clear(&v->vf[1]);
mdct_clear(&v->vm[0]);
mdct_clear(&v->vm[1]);
memset(v,0,sizeof(vorbis_dsp_state));
}
}
......@@ -524,138 +523,3 @@ int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes){
return(0);
}
#ifdef _V_SELFTEST
#include <stdio.h>
#include <math.h>
/* basic test of PCM blocking:
construct a PCM vector and block it using preset sizing in our fake
delta/multiplier generation. Immediately hand the block over to
'synthesis' and rebuild it. */
int main(){
int blocksize=1024;
int fini=100*1024;
vorbis_dsp_state encode,decode;
vorbis_info vi;
vorbis_block vb;
long counterin=0;
long countermid=0;
long counterout=0;
int done=0;
char *temp[]={ "Test" ,"the Test band", "test records",NULL };
int frame=0;
MDCT_lookup *ml[2];
vi.channels=2;
vi.rate=44100;
vi.version=0;
vi.mode=0;
vi.user_comments=temp;
vi.vendor="Xiphophorus";
vorbis_analysis_init(&encode,&vi);
vorbis_synthesis_init(&decode,&vi);
memset(&vb,0,sizeof(vorbis_block));
vorbis_block_init(&encode,&vb);
ml[0]=MDCT_init(encode.block_size[0]);
ml[1]=MDCT_init(encode.block_size[1]);
/* Submit 100K samples of data reading out blocks... */
while(!done){
int i;
double **buf=vorbis_analysis_buffer(&encode,blocksize);
for(i=0;i<blocksize;i++){
buf[0][i]=sin((counterin+i)%500/500.*M_PI*2)+2;
buf[1][i]=-1;
if((counterin+i)%15000>13000)buf[0][i]+=10;
}
i=(counterin+blocksize>fini?fini-counterin:blocksize);
vorbis_analysis_wrote(&encode,i);
counterin+=i;
while(vorbis_analysis_blockout(&encode,&vb)){
double **pcm;
int avail;
/* temp fixup */
double *window=encode.window[vb.W][vb.lW][vb.nW];
/* apply envelope */
_ve_envelope_sparsify(&vb);
_ve_envelope_apply(&vb,0);
for(i=0;i<vb.pcm_channels;i++)
MDCT(vb.pcm[i],vb.pcm[i],ml[vb.W],window);
for(i=0;i<vb.pcm_channels;i++)
iMDCT(vb.pcm[i],vb.pcm[i],ml[vb.W],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);
countermid+=encode.block_size[vb.W]/4+encode.block_size[vb.nW]/4;
}
_ve_envelope_apply(&vb,1);
vorbis_synthesis_blockin(&decode,&vb);
while((avail=vorbis_synthesis_pcmout(&decode,&pcm))){
FILE *out;
char path[80];
int i;
sprintf(path,"syn%d",frame);
out=fopen(path,"w");
for(i=0;i<avail;i++)
fprintf(out,"%ld %g\n",i+counterout,pcm[0][i]);
fprintf(out,"\n");
for(i=0;i<avail;i++)
fprintf(out,"%ld %g\n",i+counterout,pcm[1][i]);
fclose(out);
vorbis_synthesis_read(&decode,avail);
counterout+=avail;
frame++;
}
if(vb.eofflag){
done=1;
break;
}
}
}
return 0;
}
#endif
......@@ -21,6 +21,15 @@
#ifndef _vorbis_codec_h_
#define _vorbis_codec_h_
#include "mdct.h"
#include "smallft.h"
typedef struct {
int winlen;
double *window;
} envelope_lookup;
typedef struct {
long endbyte;
int endbit;
......@@ -107,17 +116,15 @@ typedef struct {
int bodybytes;
} ogg_sync_state;
typedef struct {
int winlen;
double *window;
} envelope_lookup;
typedef struct vorbis_dsp_state{
int samples_per_envelope_step;
int block_size[2];
double *window[2][2][2]; /* windowsize, leadin, leadout */
envelope_lookup ve;
drft_lookup vf[2];
mdct_lookup vm[2];
vorbis_info vi;
double **pcm;
......@@ -232,6 +239,8 @@ extern int vorbis_synthesis_pcmout(vorbis_dsp_state *v,double ***pcm);
extern int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes);
extern int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb);
#endif
......@@ -14,7 +14,7 @@
function: PCM data envelope analysis and manipulation
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Aug 05 1999
last modification date: Aug 07 1999
Vorbis manipulates the dynamic range of the incoming PCM data
envelope to minimise time-domain energy leakage from percussive and
......@@ -43,6 +43,11 @@ void _ve_envelope_init(envelope_lookup *e,int samples_per){
}
}
void _ve_envelope_clear(envelope_lookup *e){
if(e->window)free(e->window);
memset(e,0,sizeof(envelope_lookup));
}
/* initial and final blocks are special cases. Eg:
______
`--_
......@@ -64,46 +69,32 @@ void _ve_envelope_init(envelope_lookup *e,int samples_per){
span the threshhold (assuming the threshhold is active), we use an
abbreviated vector */
static void _ve_envelope_generate(double *mult,double *env,double *look,
int n,int step){
static int _ve_envelope_generate(double *mult,double *env,double *look,
int n,int step){
int i,j,p;
double m;
n*=step;
double m,mo;
for(j=0;j<n;j++)if(mult[j]!=1)break;
if(j==n)return(0);
n*=step;
/* first multiplier special case */
m=ldexp(2,mult[0]-1);
for(i=0;i<step/2;i++)env[i]=m;
p=i;
for(i=step;i<step*2;i++,p++)env[p]=m*look[i];
m=ldexp(1,mult[0]-1);
for(p=0;p<step/2;p++)env[p]=m;
/* mid multipliers normal case */
for(j=1;p<n-step/2;j++){
p-=step;
m=ldexp(2,mult[j]-1);
for(i=0;i<step;i++,p++)env[p]+=m*look[i];
for(;i<step*2;i++,p++)env[p]=m*look[i];
mo=m;
m=ldexp(1,mult[j]-1);
if(mo==m)
for(i=0;i<step;i++,p++)env[p]=m;
else
for(i=0;i<step;i++,p++)env[p]=m*look[i]+mo*look[i+step];
}
/* last multiplier special case */
p-=step;
m=ldexp(2,mult[j]-1);
for(i=0;i<step;i++,p++)env[p]+=m*look[i];
for(;p<n;p++)env[p]=m;
{
static int frameno=0;
FILE *out;
char path[80];
int i;
sprintf(path,"env%d",frameno);
out=fopen(path,"w");
for(i=0;i<n;i++)
fprintf(out,"%g\n",env[i]);
fclose(out);
frameno++;
}
return(1);
}
/* right now, we do things simple and dirty (read: our current preecho
......@@ -235,21 +226,22 @@ void _ve_envelope_apply(vorbis_block *vb,int multp){
}
/* compute the envelope curve */
_ve_envelope_generate(mult,env,look->window,vb->multend,vi->envelopesa);
/* apply the envelope curve */
for(j=0;j<vi->channels;j++){
/* check to see if the generated envelope applies to this channel */
if(vi->envelopemap[j]==i){
if(_ve_envelope_generate(mult,env,look->window,vb->multend,
vi->envelopesa)){
/* apply the envelope curve */
for(j=0;j<vi->channels;j++){
if(multp)
for(k=0;k<vb->multend*vi->envelopesa;k++)
vb->pcm[j][k]*=env[k];
else
for(k=0;k<vb->multend*vi->envelopesa;k++)
vb->pcm[j][k]/=env[k];
/* check to see if the generated envelope applies to this channel */
if(vi->envelopemap[j]==i){
if(multp)
for(k=0;k<vb->multend*vi->envelopesa;k++)
vb->pcm[j][k]*=env[k];
else
for(k=0;k<vb->multend*vi->envelopesa;k++)
vb->pcm[j][k]/=env[k];
}
}
}
}
......
......@@ -14,7 +14,7 @@
function: PCM data envelope analysis and manipulation
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Aug 06 1999
last modification date: Aug 07 1999
********************************************************************/
......@@ -25,6 +25,7 @@ extern void _ve_envelope_multipliers(vorbis_dsp_state *v);
extern void _ve_envelope_apply(vorbis_block *vb,int multp);
extern void _ve_envelope_sparsify(vorbis_block *vb);
extern void _ve_envelope_init(envelope_lookup *e,int samples_per);
extern void _ve_envelope_clear(envelope_lookup *e);
#endif
......@@ -40,17 +40,14 @@
#define VORBIS_SPECIFIC_MODIFICATIONS
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "mdct.h"
static double oPI = 3.14159265358979323846;
/* build lookups for trig functions; also pre-figure scaling and
some window function algebra. */
MDCT_lookup *MDCT_init(int n){
MDCT_lookup *lookup=malloc(sizeof(MDCT_lookup));
void mdct_init(mdct_lookup *lookup,int n){
int *bitrev=malloc(sizeof(int)*(n/4));
double *trig=malloc(sizeof(double)*(n+n/4));
double *AE=trig;
......@@ -72,14 +69,14 @@ MDCT_lookup *MDCT_init(int n){
/* trig lookups... */
for(i=0;i<n/4;i++){
AE[i]=cos((oPI/n)*(4*i));
AO[i]=-sin((oPI/n)*(4*i));
BE[i]=cos((oPI/(2*n))*(2*i+1));
BO[i]=sin((oPI/(2*n))*(2*i+1));
AE[i]=cos((M_PI/n)*(4*i));
AO[i]=-sin((M_PI/n)*(4*i));
BE[i]=cos((M_PI/(2*n))*(2*i+1));
BO[i]=sin((M_PI/(2*n))*(2*i+1));
}
for(i=0;i<n/8;i++){
CE[i]=cos((oPI/n)*(4*i+2));
CO[i]=-sin((oPI/n)*(4*i+2));
CE[i]=cos((M_PI/n)*(4*i+2));
CO[i]=-sin((M_PI/n)*(4*i+2));
}
/* bitreverse lookup... */
......@@ -95,21 +92,19 @@ MDCT_lookup *MDCT_init(int n){
bitB[i]=acc*2;
}
}
return(lookup);
}
void MDCT_free(MDCT_lookup *l){
void mdct_clear(mdct_lookup *l){
if(l){
if(l->trig)free(l->trig);
if(l->bitrev)free(l->bitrev);
free(l);
memset(l,0,sizeof(mdct_lookup));
}
}
static inline void _MDCT_kernel(double *x,
static inline void _mdct_kernel(double *x,
int n, int n2, int n4, int n8,
MDCT_lookup *init){
mdct_lookup *init){
double *w=x+1; /* interleaved access improves cache locality */
int i;
/* step 2 */
......@@ -205,7 +200,7 @@ static inline void _MDCT_kernel(double *x,
}
}
void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
void mdct_forward(mdct_lookup *init, double *in, double *out, double *window){
int n=init->n;
double *x=alloca(n*sizeof(double));
int n2=n>>1;
......@@ -255,7 +250,7 @@ void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
}
}
_MDCT_kernel(x,n,n2,n4,n8,init);
_mdct_kernel(x,n,n2,n4,n8,init);
/* step 8 */
......@@ -273,7 +268,7 @@ void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
}
}
void iMDCT(double *in, double *out, MDCT_lookup *init, double *window){
void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
int n=init->n;
double *x=alloca(n*sizeof(double));
int n2=n>>1;
......@@ -310,7 +305,7 @@ void iMDCT(double *in, double *out, MDCT_lookup *init, double *window){
}
_MDCT_kernel(x,n,n2,n4,n8,init);
_mdct_kernel(x,n,n2,n4,n8,init);
/* step 8 */
......
......@@ -15,8 +15,8 @@
********************************************************************/
#ifndef _OGG_MDCT_H_
#define _OGG_MDCT_H_
#ifndef _OGG_mdct_H_
#define _OGG_mdct_H_
typedef struct {
int n;
......@@ -25,12 +25,14 @@ typedef struct {
double *trig;
int *bitrev;
} MDCT_lookup;
} mdct_lookup;
extern MDCT_lookup *MDCT_init(int n);
extern void MDCT_free(MDCT_lookup *l);
extern void MDCT(double *in, double *out, MDCT_lookup *init, double *window);
extern void iMDCT(double *in, double *out, MDCT_lookup *init, double *window);
extern void mdct_init(mdct_lookup *lookup,int n);
extern void mdct_clear(mdct_lookup *l);
extern void mdct_forward(mdct_lookup *init, double *in,
double *out, double *window);
extern void mdct_backward(mdct_lookup *init, double *in,
double *out, double *window);
#endif
......
......@@ -12,6 +12,7 @@
******************************************************************/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "smallft.h"
......@@ -300,11 +301,6 @@ static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
for(i=0;i<n;i++)c[i]=ch[i];
}
static void fdrfftf(int n,double *r,double *wsave,int *ifac){
if(n==1)return;
drftf1(n,r,wsave,wsave+n,ifac);
}
static void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){
int i,k,t0,t1,t2,t3,t4,t5,t6;
double ti2,tr2;
......@@ -497,52 +493,29 @@ static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
for(i=0;i<n;i++)c[i]=ch[i];
}
static void fdrfftb(int n, double *r, double *wsave, int *ifac){
if (n == 1)return;
drftb1(n, r, wsave, wsave+n, ifac);
}
void fft_forward(int n, double *buf,double *trigcache,int *splitcache){
int flag=0;
if(!trigcache || !splitcache){
trigcache=calloc(3*n,sizeof(double));
splitcache=calloc(32,sizeof(int));
fdrffti(n, trigcache, splitcache);
flag=1;
}
fdrfftf(n, buf, trigcache, splitcache);
if(flag){
free(trigcache);
free(splitcache);
}
void drft_forward(drft_lookup *l,double *data){
if(l->n==1)return;
drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
}
void fft_backward(int n, double *buf, double *trigcache,int *splitcache){
void drft_backward(drft_lookup *l,double *data){
int i;
int flag=0;
if(!trigcache || !splitcache){
trigcache=calloc(3*n,sizeof(double));
splitcache=calloc(32,sizeof(int));
fdrffti(n, trigcache, splitcache);
flag=1;
}
fdrfftb(n, buf, trigcache, splitcache);
for(i=0;i<n;i++)buf[i]/=n;
if (l->n==1)return;
drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
for(i=0;i<l->n;i++)data[i]/=l->n;
}
if(flag){
free(trigcache);
free(splitcache);
}
void drft_init(drft_lookup *l,int n){
l->n=n;
l->trigcache=calloc(3*n,sizeof(double));
l->splitcache=calloc(32,sizeof(int));
fdrffti(n, l->trigcache, l->splitcache);
}
void fft_i(int n, double **trigcache, int **splitcache){
*trigcache=calloc(3*n,sizeof(double));
*splitcache=calloc(32,sizeof(int));
fdrffti(n, *trigcache, *splitcache);
void drft_clear(drft_lookup *l){
if(l){
if(l->trigcache)free(l->trigcache);
if(l->splitcache)free(l->splitcache);
memset(l,0,sizeof(drft_lookup));
}
}
/******************************************************************
* CopyPolicy: GNU Public License 2 applies
* Copyright (C) 1998 Monty xiphmont@mit.edu, monty@xiph.org
* Copyright (C) 1998-1999 Monty xiphmont@mit.edu, monty@xiph.org