Commit 117b96ee authored by Monty's avatar Monty
Browse files

First shot at the continuous balance code.  Analysis is not correct yet.

Monty

svn path=/trunk/vorbis/; revision=80
parent 9e5b2114
......@@ -14,122 +14,220 @@
function: single-block PCM analysis
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Aug 09 1999
last modification date: Aug 22 1999
********************************************************************/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "lpc.h"
#include "lsp.h"
#include "codec.h"
#include "envelope.h"
#include "mdct.h"
#include "psy.h"
extern void compute_balance(double *A, double *B, double *phi,int n);
int vorbis_analysis(vorbis_block *vb){
int i;
double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
lpc_lookup *vl=&vb->vd->vl[vb->W];
lpc_lookup *vbal=&vb->vd->vbal[vb->W];
int n=vb->pcmend;
static int frameno=0;
double balance_v[vbal->m];
double balance_amp;
_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 */
/* hardwired to stereo for 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;
}
double *C=vb->pcm[0];
double *D=vb->pcm[1];
/* Balance */
balance_amp=_vp_balance_compute(D,C,balance_v,vbal);
/*{
FILE *out;
char buffer[80];
sprintf(buffer,"com%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++){
fprintf(out," 0. 0.\n");
fprintf(out,"%g %g\n",C[i],D[i]);
fprintf(out,"\n");
}
fclose(out);
sprintf(buffer,"L%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++){
fprintf(out,"%g\n",C[i]);
}
fclose(out);
sprintf(buffer,"R%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++){
fprintf(out,"%g\n",D[i]);
}
fclose(out);
}*/
/* apply balance vectors, mix down to the vectors we actually encode */
_vp_balance_apply(D,C,balance_v,balance_amp,vbal,1);
/*{
FILE *out;
char buffer[80];
sprintf(buffer,"bal%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++){
fprintf(out," 0. 0.\n");
fprintf(out,"%g %g\n",C[i],D[i]);
fprintf(out,"\n");
}
fclose(out);
sprintf(buffer,"C%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++){
fprintf(out,"%g\n",C[i]);
}
fclose(out);
sprintf(buffer,"D%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++){
fprintf(out,"%g\n",D[i]);
}
fclose(out);
}*/
/* 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;
double lpc1[80];
double lsp1[80];
double lsp2[80];
double amp1;
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_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));
_vp_mask_floor(D,floor1,n/2);
memcpy(curve,floor1,sizeof(double)*n/2);
amp1=sqrt(vorbis_curve_to_lpc(curve,lpc1,vl));
/*vorbis_lpc_to_lsp(lpc1,lsp1,30);
{
int scale=1020;
int last=0;
for(i=0;i<30;i++){
double q=lsp1[i]/M_PI*scale;
int val=rint(q-last);
last+=val;
lsp1[i]=val;
}
}
{
int scale=1020;
double last=0;
for(i=0;i<30;i++){
last+=lsp1[i];
lsp2[i]=last*M_PI/scale;
}
}
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);*/
vorbis_lsp_to_lpc(lsp2,lpc1,30);*/
vorbis_lpc_to_curve(work,lpc1,amp1,vl);
_vp_psy_quantize(C,work,n/2);
_vp_psy_quantize(D,work,n/2);
_vp_psy_unquantize(C,work,n/2);
_vp_psy_unquantize(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++)
/*sprintf(buffer,"qC%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,"qD%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++)
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++)
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);
fclose(out);
sprintf(buffer,"curve%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<n/2;i++)
fprintf(out,"%g\n",curve[i]);
fclose(out);
sprintf(buffer,"lsp%d.m",frameno);
out=fopen(buffer,"w+");
for(i=0;i<30;i++){
fprintf(out,"%g 0.\n",lsp1[i]);
fprintf(out,"%g .1\n",lsp1[i]);
fprintf(out,"\n");
fprintf(out,"\n");
}
fclose(out);*/
frameno++;
}
}
}
/* unmix */
_vp_balance_apply(D,C,balance_v,balance_amp,vbal,0);
}
return(0);
}
......
......@@ -142,8 +142,18 @@ 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){
_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 */
/* Yes, wasteful to have four lookups. This will get collapsed once
things crystallize */
lpc_init(&v->vl[0],vi->smallblock/2,vi->smallblock/2,
vi->floororder,vi->flooroctaves,1);
lpc_init(&v->vl[1],vi->largeblock/2,vi->largeblock/2,
vi->floororder,vi->flooroctaves,1);
lpc_init(&v->vbal[0],vi->smallblock/2,256,
vi->balanceorder,vi->balanceoctaves,1);
lpc_init(&v->vbal[1],vi->largeblock/2,256,
vi->balanceorder,vi->balanceoctaves,1);
return(0);
}
......@@ -168,8 +178,6 @@ void vorbis_analysis_clear(vorbis_dsp_state *v){
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));
......@@ -336,6 +344,7 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
}
}else
v->nW=1;
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
......@@ -435,6 +444,18 @@ int vorbis_synthesis_init(vorbis_dsp_state *v,vorbis_info *vi){
_vds_shared_init(v,vi);
vi->envelopech=temp;
/* Yes, wasteful to have four lookups. This will get collapsed once
things crystallize */
lpc_init(&v->vl[0],vi->smallblock/2,vi->smallblock/2,
vi->floororder,vi->flooroctaves,0);
lpc_init(&v->vl[1],vi->largeblock/2,vi->largeblock/2,
vi->floororder,vi->flooroctaves,0);
lpc_init(&v->vbal[0],vi->smallblock/2,256,
vi->balanceorder,vi->balanceoctaves,0);
lpc_init(&v->vbal[1],vi->largeblock/2,256,
vi->balanceorder,vi->balanceoctaves,0);
/* Adjust centerW to allow an easier mechanism for determining output */
v->pcm_returned=v->centerW;
v->centerW-= v->block_size[v->W]/4+v->block_size[v->lW]/4;
......
......@@ -14,21 +14,48 @@
function: codec headers
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Jul 28 1999
last modification date: Aug 21 1999
********************************************************************/
#ifndef _vorbis_codec_h_
#define _vorbis_codec_h_
#include "mdct.h"
#include "smallft.h"
typedef struct {
int n;
int log2n;
double *trig;
int *bitrev;
} mdct_lookup;
typedef struct {
int n;
double *trigcache;
int *splitcache;
} drft_lookup;
typedef struct {
int winlen;
double *window;
} envelope_lookup;
typedef struct lpclook{
/* encode lookups */
int *bscale;
double *escale;
drft_lookup fft;
/* en/decode lookups */
double *dscale;
double *norm;
int n;
int ln;
int m;
} lpc_lookup;
typedef struct {
long endbyte;
......@@ -56,9 +83,13 @@ typedef struct vorbis_info{
int **Echannelmap; /* which encoding channels produce what pcm (decode) */
int **channelmap; /* which pcm channels produce what floors (encode) */
int floororder;
int flooroctaves;
int floorch;
int *floormap;
int balanceorder;
int balanceoctaves;
double preecho_thresh;
double preecho_clamp;
} vorbis_info;
......@@ -126,8 +157,9 @@ typedef struct vorbis_dsp_state{
double *window[2][2][2]; /* windowsize, leadin, leadout */
envelope_lookup ve;
drft_lookup vf[2];
mdct_lookup vm[2];
lpc_lookup vl[2];
lpc_lookup vbal[2];
vorbis_info vi;
......
......@@ -15,6 +15,7 @@ int main(){
int mtemp[]={0,1};
int mtemp2[]={0,1};
int frame=0;
int i;
signed char buffer[READ*4+44];
......@@ -30,8 +31,13 @@ int main(){
vi.envelopech=2;
vi.envelopemap=mtemp;
vi.floormap=mtemp2;
vi.floororder=20;
vi.floororder=30;
vi.flooroctaves=5;
vi.floorch=2;
vi.balanceorder=4;
vi.balanceoctaves=5;
vi.channelmap=NULL;
vi.preecho_thresh=10.;
vi.preecho_clamp=.5;
......@@ -43,11 +49,12 @@ int main(){
fread(buffer,1,44,stdin);
fwrite(buffer,1,44,stdout);
for(i=0;i<0;i++)
fread(buffer,1,READ*4,stdin);
while(!done){
long bread=fread(buffer,1,READ*4,stdin);
double **buf=vorbis_analysis_buffer(&encode,READ);
long i;
/* uninterleave samples */
......
......@@ -60,60 +60,15 @@ Carsten Bormann
coefficients. Autocorrelation LPC coeff generation algorithm
invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
/* Input : n element envelope curve
Output: m lpc coefficients, excitation energy */
/* Input : n element envelope curve
Output: m lpc coefficients, excitation energy */
double memcof(double *data, int n, int m, double *d){
int k,j,i;
double p=0.,wk1[n],wk2[n],wkm[m],xms;
memset(wk1,0,sizeof(wk1));
memset(wk2,0,sizeof(wk2));
memset(wkm,0,sizeof(wkm));
for (j=0;j<n;j++) p += data[j]*data[j];
xms=p/n;
wk1[0]=data[0];
wk2[n-2]=data[n-1];
for (j=2;j<=n-1;j++) {
wk1[j-1]=data[j-1];
wk2[j-2]=data[j-1];
}
for (k=1;k<=m;k++) {
double num=0.,denom=0.;
for (j=1;j<=(n-k);j++) {
num += wk1[j-1]*wk2[j-1];
denom += wk1[j-1]*wk1[j-1] + wk2[j-1]*wk2[j-1];
}
d[k-1]=2.0*num/denom;
xms *= (1.0-d[k-1]*d[k-1]);
for (i=1;i<=(k-1);i++)
d[i-1]=wkm[i-1]-d[k-1]*wkm[k-i-1];
if (k == m) return xms;
for (i=1;i<=k;i++) wkm[i-1]=d[i-1];
for (j=1;j<=(n-k-1);j++) {
wk1[j-1] -= wkm[k-1]*wk2[j-1];
wk2[j-1]=wk2[j]-wkm[k-1]*wk1[j];
}
}
}
static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l){
int n=l->ln;
int m=l->m;
double aut[m+1],work[n+n],error;
drft_lookup dl;
int i,j;
/* input is a real curve. make it complex-real */
for(i=0;i<n;i++){
work[i*2]=curve[i];
......@@ -121,9 +76,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
}
n*=2;
drft_init(&dl,n);
drft_backward(&dl,work);
drft_clear(&dl);
drft_backward(&l->fft,work);
/* The autocorrelation will not be circular. Shift, else we lose
most of the power in the edges. */
......@@ -133,7 +86,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
work[i++]=work[j];
work[j++]=temp;
}
/* autocorrelation, p+1 lag coefficients */
j=m+1;
......@@ -177,7 +130,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
/* we need the error value to know how big an impulse to hit the
filter with later */
return error;
}
......@@ -187,7 +140,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
looks at the below in the context of the calling function, there's
lots of redundant trig, but at least it's clear */
static double vorbis_lpc_magnitude(double w,double *lpc, int m){
double vorbis_lpc_magnitude(double w,double *lpc, int m){
int k;
double real=1,imag=0;
double wn=w;
......@@ -228,14 +181,16 @@ static double vorbis_lpc_magnitude(double w,double *lpc, int m){
oct == octaves in normalized scale, encode_p == encode (1) or
decode (0) */
double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
void lpc_init(lpc_lookup *l,int n, int mapped, int m, int oct, int encode_p){
double bias=LOG_BIAS(n,oct);
double scale=(float)n/(float)oct; /* where n==mapped */
double scale=(float)mapped/(float)oct; /* where n==mapped */
int i;
memset(l,0,sizeof(lpc_lookup));
l->n=n;
l->ln=mapped;
l->m=m;
l->escale=malloc(n*sizeof(double));
l->dscale=malloc(n*sizeof(double));
l->norm=malloc(n*sizeof(double));
......@@ -251,10 +206,15 @@ double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
if(encode_p){
/* encode */
l->bscale=malloc(n*sizeof(int));
l->escale=malloc(n*sizeof(double));
for(i=0;i<n;i++)
for(i=0;i<n;i++){
l->escale[i]=LINEAR_X(i/scale,bias);
l->bscale[i]=rint(LOG_X(i,bias)*scale);
}
drft_init(&l->fft,mapped*2);
}
/* decode; encode may use this too */
......@@ -267,7 +227,9 @@ double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
void lpc_clear(lpc_lookup *l){
if(l){
if(l->bscale)free(l->bscale);
if(l->escale)free(l->escale);
drft_clear(&l->fft);
free(l->dscale);
free(l->norm);
}
......@@ -281,10 +243,10 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
/* map the input curve to a log curve for encoding */
/* for clarity, mapped and n are both represented although setting
'em equal is a decent rule of thumb. */
'em equal is a decent rule of thumb. The below must be reworked
slightly if mapped != n */
int n=l->n;
int m=l->m;
int mapped=n;
double work[mapped];
int i;
......@@ -305,7 +267,7 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
memcpy(curve,work,sizeof(work));
return vorbis_gen_lpc(work,mapped,lpc,l->m);
return vorbis_gen_lpc(work,lpc,l);
}
/* generate the whole freq response curve on an LPC IIR filter */
......
......@@ -20,9 +20,15 @@
#include "codec.h"
extern double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p);
extern void lpc_init(lpc_lookup *l,int n, int mapped,
int m, int oct, int encode_p);
extern void lpc_clear(lpc_lookup *l);
/* simple linear scale LPC code */