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);
}
/* 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);
}
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));
}
}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
_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;
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;
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;
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;
int mtemp[]={0,1};
int mtemp2[]={0,1};
int frame=0;
+ int i;
signed char buffer[READ*4+44];
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;
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 */
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];
}
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. */
work[i++]=work[j];
work[j++]=temp;
}
-
+
/* autocorrelation, p+1 lag coefficients */
j=m+1;
/* we need the error value to know how big an impulse to hit the
filter with later */
-
+
return error;
}
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;
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));
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 */
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);
}
/* 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;
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 */
#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 */
+extern double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l);
+extern double vorbis_lpc_magnitude(double w,double *lpc, int m);
+
+/* log scale layer */
extern double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l);
extern void vorbis_lpc_to_curve(double *curve,double *lpc, double amp,
lpc_lookup *l);
lsp[i+1] = acos(g2r[i/2]*.5);
}
}
-
extern void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m);
extern void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m);
-
+
#endif
#ifndef _OGG_mdct_H_
#define _OGG_mdct_H_
-typedef struct {
- int n;
- int log2n;
-
- double *trig;
- int *bitrev;
-
-} mdct_lookup;
+#include "codec.h"
extern void mdct_init(mdct_lookup *lookup,int n);
extern void mdct_clear(mdct_lookup *l);
function: random psychoacoustics (not including preecho)
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Aug 08 1999
+ last modification date: Aug 26 1999
********************************************************************/
#include <math.h>
+#include <string.h>
+#include "stdio.h"
#include "codec.h"
#include "psy.h"
+#include "lpc.h"
+#include "smallft.h"
+#include "xlogmap.h"
-#define NOISEdB 6
+#define NOISEdB 10
-#define MASKdB 12
+#define MASKdB 18
#define HROLL 60
#define LROLL 90
-#define MASKBIAS 40
+#define MASKBIAS 10
-#define LNOISE .8
+#define LNOISE .95
#define HNOISE 1.01
#define NOISEBIAS 20
double acc=0,div=0;
int i,j;
- for(i=0;i<n;i++){
+ for(i=100;i<n;i++){
long newlo=i*LNOISE-NOISEBIAS;
long newhi=i*HNOISE+NOISEBIAS;
double temp;
int val=rint(f[i]/m[i]);
if(val>16)val=16;
if(val<-16)val=-16;
- f[i]=val*m[i];
+ if(val==0 || val==2){
+ if(f[i]<0){
+ f[i]=-1;
+ }else{
+ f[i]=1;
+ }
+ }
+
+ f[i]=val;
}
}
+void _vp_psy_unquantize(double *f, double *m,int n){
+ int i;
+ for(i=0;i<n;i++)
+ f[i]*=m[i];
+}
void _vp_psy_sparsify(double *f, double *m,int n){
int i;
for(i=0;i<n;i++)
- if(fabs(f[i])<m[i])f[i]=0;
+ if(fabs(f[i])<m[i]*.5)f[i]=0;
+}
+
+/* s must be padded at the end with m-1 zeroes */
+static void time_convolve(double *s,double *r,int n,int m){
+ int i;
+
+ for(i=0;i<n;i++){
+ int j;
+ double acc=0;
+
+ for(j=0;j<m;j++)
+ acc+=s[i+j]*r[m-j-1];
+
+ s[i]=acc;
+ }
}
-void _vp_psy_make_lsp(vorbis_block *vb){
+/*************************************************************************/
+/* Continuous balance analysis/synthesis *********************************/
+
+/* Compute the average continual spectral balance of the given vectors
+ (in radians); encode into LPC coefficients */
+
+double _vp_balance_compute(double *A, double *B, double *lpc,lpc_lookup *vb){
+ /* correlate in time (the response function is small). Log
+ frequency scale, small mapping */
+
+ int n=vb->n;
+ int mapped=vb->ln;
+
+ /* 256/15 are arbitrary but unimportant to decoding */
+ int resp=15;
+ int i;
+ /* This is encode side. Don't think too hard about it */
+
+ double workA[mapped+resp];
+ double workB[mapped+resp];
+ double p[mapped];
+ double workC[resp];
+ memset(workA,0,sizeof(workA));
+ memset(workB,0,sizeof(workB));
+ memset(workC,0,sizeof(workC));
+
+ for(i=0;i<n;i++){
+ int j=vb->bscale[i]+(resp>>1);
+ double mag_sq=A[i]*A[i]+B[i]*B[i];
+ double phi;
+
+ if(B[i]==0)
+ phi=M_PI/2.;
+ else{
+ phi=atan(A[i]/B[i]);
+ if((A[i]<0 && B[i]>0)||
+ (A[i]>0 && B[i]<0)){
+ /* rotate II and IV into the first quadrant. III is already there */
+ phi+=M_PI/2;
+ }
+ }
+
+ workA[j]+=mag_sq*sin(phi);
+ workB[j]+=mag_sq*cos(phi);
+ }
+
+ /* prepare convolution vector. Play with a few different shapes */
+
+ for(i=0;i<resp;i++){
+ workC[i]=sin(M_PI*(i+1)/(float)(resp+1));
+ workC[i]*=workC[i];
+ }
+
+ time_convolve(workA,workC,mapped,resp);
+ time_convolve(workB,workC,mapped,resp);
+
+ {
+ double amp1;
+
+ for(i=0;i<mapped;i++){
+ p[i]=atan2(workA[i],workB[i]);
+ }
+
+ amp1=sqrt(vorbis_gen_lpc(p,lpc,vb));
+
+ return(amp1);
+ }
}
+void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
+ lpc_lookup *vb,int divp){
+ int i;
+ for(i=0;i<vb->n;i++){
+ double mag=sqrt(A[i]*A[i]+B[i]*B[i]);
+ double del=vorbis_lpc_magnitude(vb->dscale[i],lpc,vb->m)*amp;
+ double phi=atan2(A[i],B[i]);
+
+ if(divp)
+ phi-=del;
+ else
+ phi+=del;
+
+ A[i]=mag*sin(phi);
+ B[i]=mag*cos(phi);
+ }
+}
#ifndef _V_PSY_H_
#define _V_PSY_H_
+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,
+ lpc_lookup *vb,int divp);
#endif
}
}
+static void dradfg(int ido,int ip,int l1,int idl1,double *cc,double *c1,
+ double *c2,double *ch,double *ch2,double *wa){
+
+ static double tpi=6.28318530717958647692528676655900577;
+ int idij,ipph,i,j,k,l,ic,ik,is;
+ int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+ double dc2,ai1,ai2,ar1,ar2,ds2;
+ int nbd;
+ double dcp,arg,dsp,ar1h,ar2h;
+ int idp2,ipp2;
+
+ arg=tpi/(double)ip;
+ dcp=cos(arg);
+ dsp=sin(arg);
+ ipph=(ip+1)>>1;
+ ipp2=ip;
+ idp2=ido;
+ nbd=(ido-1)>>1;
+ t0=l1*ido;
+ t10=ip*ido;
+
+ if(ido==1)goto L119;
+ for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
+
+ t1=0;
+ for(j=1;j<ip;j++){
+ t1+=t0;
+ t2=t1;
+ for(k=0;k<l1;k++){
+ ch[t2]=c1[t2];
+ t2+=ido;
+ }
+ }
+
+ is=-ido;
+ t1=0;
+ if(nbd>l1){
+ for(j=1;j<ip;j++){
+ t1+=t0;
+ is+=ido;
+ t2= -ido+t1;
+ for(k=0;k<l1;k++){
+ idij=is-1;
+ t2+=ido;
+ t3=t2;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t3+=2;
+ ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+ ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+ }
+ }
+ }
+ }else{
+
+ for(j=1;j<ip;j++){
+ is+=ido;
+ idij=is-1;
+ t1+=t0;
+ t2=t1;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t2+=2;
+ t3=t2;
+ for(k=0;k<l1;k++){
+ ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+ ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+ t3+=ido;
+ }
+ }
+ }
+ }
+
+ t1=0;
+ t2=ipp2*t0;
+ if(nbd<l1){
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5=t3-ido;
+ t6=t4-ido;
+ for(k=0;k<l1;k++){
+ t5+=ido;
+ t6+=ido;
+ c1[t5-1]=ch[t5-1]+ch[t6-1];
+ c1[t6-1]=ch[t5]-ch[t6];
+ c1[t5]=ch[t5]+ch[t6];
+ c1[t6]=ch[t6-1]-ch[t5-1];
+ }
+ }
+ }
+ }else{
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ c1[t5-1]=ch[t5-1]+ch[t6-1];
+ c1[t6-1]=ch[t5]-ch[t6];
+ c1[t5]=ch[t5]+ch[t6];
+ c1[t6]=ch[t6-1]-ch[t5-1];
+ }
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+ }
+
+L119:
+ for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+ t1=0;
+ t2=ipp2*idl1;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1-ido;
+ t4=t2-ido;
+ for(k=0;k<l1;k++){
+ t3+=ido;
+ t4+=ido;
+ c1[t3]=ch[t3]+ch[t4];
+ c1[t4]=ch[t4]-ch[t3];
+ }
+ }
+
+ ar1=1.;
+ ai1=0.;
+ t1=0;
+ t2=ipp2*idl1;
+ t3=(ip-1)*idl1;
+ for(l=1;l<ipph;l++){
+ t1+=idl1;
+ t2-=idl1;
+ ar1h=dcp*ar1-dsp*ai1;
+ ai1=dcp*ai1+dsp*ar1;
+ ar1=ar1h;
+ t4=t1;
+ t5=t2;
+ t6=t3;
+ t7=idl1;
+
+ for(ik=0;ik<idl1;ik++){
+ ch2[t4++]=c2[ik]+ar1*c2[t7++];
+ ch2[t5++]=ai1*c2[t6++];
+ }
+
+ dc2=ar1;
+ ds2=ai1;
+ ar2=ar1;
+ ai2=ai1;
+
+ t4=idl1;
+ t5=(ipp2-1)*idl1;
+ for(j=2;j<ipph;j++){
+ t4+=idl1;
+ t5-=idl1;
+
+ ar2h=dc2*ar2-ds2*ai2;
+ ai2=dc2*ai2+ds2*ar2;
+ ar2=ar2h;
+
+ t6=t1;
+ t7=t2;
+ t8=t4;
+ t9=t5;
+ for(ik=0;ik<idl1;ik++){
+ ch2[t6++]+=ar2*c2[t8++];
+ ch2[t7++]+=ai2*c2[t9++];
+ }
+ }
+ }
+
+ t1=0;
+ for(j=1;j<ipph;j++){
+ t1+=idl1;
+ t2=t1;
+ for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
+ }
+
+ if(ido<l1)goto L132;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t4=t2;
+ for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
+ t1+=ido;
+ t2+=t10;
+ }
+
+ goto L135;
+
+ L132:
+ for(i=0;i<ido;i++){
+ t1=i;
+ t2=i;
+ for(k=0;k<l1;k++){
+ cc[t2]=ch[t1];
+ t1+=ido;
+ t2+=t10;
+ }
+ }
+
+ L135:
+ t1=0;
+ t2=ido<<1;
+ t3=0;
+ t4=ipp2*t0;
+ for(j=1;j<ipph;j++){
+
+ t1+=t2;
+ t3+=t0;
+ t4-=t0;
+
+ t5=t1;
+ t6=t3;
+ t7=t4;
+
+ for(k=0;k<l1;k++){
+ cc[t5-1]=ch[t6];
+ cc[t5]=ch[t7];
+ t5+=t10;
+ t6+=ido;
+ t7+=ido;
+ }
+ }
+
+ if(ido==1)return;
+ if(nbd<l1)goto L141;
+
+ t1=-ido;
+ t3=0;
+ t4=0;
+ t5=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t2;
+ t3+=t2;
+ t4+=t0;
+ t5-=t0;
+ t6=t1;
+ t7=t3;
+ t8=t4;
+ t9=t5;
+ for(k=0;k<l1;k++){
+ for(i=2;i<ido;i+=2){
+ ic=idp2-i;
+ cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
+ cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
+ cc[i+t7]=ch[i+t8]+ch[i+t9];
+ cc[ic+t6]=ch[i+t9]-ch[i+t8];
+ }
+ t6+=t10;
+ t7+=t10;
+ t8+=ido;
+ t9+=ido;
+ }
+ }
+ return;
+
+ L141:
+
+ t1=-ido;
+ t3=0;
+ t4=0;
+ t5=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t2;
+ t3+=t2;
+ t4+=t0;
+ t5-=t0;
+ for(i=2;i<ido;i+=2){
+ t6=idp2+t1-i;
+ t7=i+t3;
+ t8=i+t4;
+ t9=i+t5;
+ for(k=0;k<l1;k++){
+ cc[t7-1]=ch[t8-1]+ch[t9-1];
+ cc[t6-1]=ch[t8-1]-ch[t9-1];
+ cc[t7]=ch[t8]+ch[t9];
+ cc[t6]=ch[t9]-ch[t8];
+ t6+=t10;
+ t7+=t10;
+ t8+=ido;
+ t9+=ido;
+ }
+ }
+ }
+}
+
static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
int i,k1,l1,l2;
int na,kh,nf;
goto L110;
L104:
- return; /* We're restricted to powers of two. just fail */
+ if(ido==1)na=1-na;
+ if(na!=0)goto L109;
+
+ dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+ na=1;
+ goto L110;
+
+ L109:
+ dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+ na=0;
L110:
l2=l1;
}
}
+static void dradb3(int ido,int l1,double *cc,double *ch,double *wa1,
+ double *wa2){
+ static double taur = -.5;
+ static double taui = .86602540378443864676372317075293618;
+ int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+ double ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
+ t0=l1*ido;
+
+ t1=0;
+ t2=t0<<1;
+ t3=ido<<1;
+ t4=ido+(ido<<1);
+ t5=0;
+ for(k=0;k<l1;k++){
+ tr2=cc[t3-1]+cc[t3-1];
+ cr2=cc[t5]+(taur*tr2);
+ ch[t1]=cc[t5]+tr2;
+ ci3=taui*(cc[t3]+cc[t3]);
+ ch[t1+t0]=cr2-ci3;
+ ch[t1+t2]=cr2+ci3;
+ t1+=ido;
+ t3+=t4;
+ t5+=t4;
+ }
+
+ if(ido==1)return;
+
+ t1=0;
+ t3=ido<<1;
+ for(k=0;k<l1;k++){
+ t7=t1+(t1<<1);
+ t6=(t5=t7+t3);
+ t8=t1;
+ t10=(t9=t1+t0)+t0;
+
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6-=2;
+ t7+=2;
+ t8+=2;
+ t9+=2;
+ t10+=2;
+ tr2=cc[t5-1]+cc[t6-1];
+ cr2=cc[t7-1]+(taur*tr2);
+ ch[t8-1]=cc[t7-1]+tr2;
+ ti2=cc[t5]-cc[t6];
+ ci2=cc[t7]+(taur*ti2);
+ ch[t8]=cc[t7]+ti2;
+ cr3=taui*(cc[t5-1]-cc[t6-1]);
+ ci3=taui*(cc[t5]+cc[t6]);
+ dr2=cr2-ci3;
+ dr3=cr2+ci3;
+ di2=ci2+cr3;
+ di3=ci2-cr3;
+ ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
+ ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
+ ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
+ ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
+ }
+ t1+=ido;
+ }
+}
+
static void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
double *wa2,double *wa3){
static double sqrt2=1.4142135623730950488016887242097;
}
}
+static void dradbg(int ido,int ip,int l1,int idl1,double *cc,double *c1,
+ double *c2,double *ch,double *ch2,double *wa){
+ static double tpi=6.28318530717958647692528676655900577;
+ int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
+ t11,t12;
+ double dc2,ai1,ai2,ar1,ar2,ds2;
+ int nbd;
+ double dcp,arg,dsp,ar1h,ar2h;
+ int ipp2;
+
+ t10=ip*ido;
+ t0=l1*ido;
+ arg=tpi/(double)ip;
+ dcp=cos(arg);
+ dsp=sin(arg);
+ nbd=(ido-1)>>1;
+ ipp2=ip;
+ ipph=(ip+1)>>1;
+ if(ido<l1)goto L103;
+
+ t1=0;
+ t2=0;
+ for(k=0;k<l1;k++){
+ t3=t1;
+ t4=t2;
+ for(i=0;i<ido;i++){
+ ch[t3]=cc[t4];
+ t3++;
+ t4++;
+ }
+ t1+=ido;
+ t2+=t10;
+ }
+ goto L106;
+
+ L103:
+ t1=0;
+ for(i=0;i<ido;i++){
+ t2=t1;
+ t3=t1;
+ for(k=0;k<l1;k++){
+ ch[t2]=cc[t3];
+ t2+=ido;
+ t3+=t10;
+ }
+ t1++;
+ }
+
+ L106:
+ t1=0;
+ t2=ipp2*t0;
+ t7=(t5=ido<<1);
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ t6=t5;
+ for(k=0;k<l1;k++){
+ ch[t3]=cc[t6-1]+cc[t6-1];
+ ch[t4]=cc[t6]+cc[t6];
+ t3+=ido;
+ t4+=ido;
+ t6+=t10;
+ }
+ t5+=t7;
+ }
+
+ if (ido == 1)goto L116;
+ if(nbd<l1)goto L112;
+
+ t1=0;
+ t2=ipp2*t0;
+ t7=0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+
+ t7+=(ido<<1);
+ t8=t7;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ t9=t8;
+ t11=t8;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ t9+=2;
+ t11-=2;
+ ch[t5-1]=cc[t9-1]+cc[t11-1];
+ ch[t6-1]=cc[t9-1]-cc[t11-1];
+ ch[t5]=cc[t9]-cc[t11];
+ ch[t6]=cc[t9]+cc[t11];
+ }
+ t3+=ido;
+ t4+=ido;
+ t8+=t10;
+ }
+ }
+ goto L116;
+
+ L112:
+ t1=0;
+ t2=ipp2*t0;
+ t7=0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ t7+=(ido<<1);
+ t8=t7;
+ t9=t7;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t8+=2;
+ t9-=2;
+ t5=t3;
+ t6=t4;
+ t11=t8;
+ t12=t9;
+ for(k=0;k<l1;k++){
+ ch[t5-1]=cc[t11-1]+cc[t12-1];
+ ch[t6-1]=cc[t11-1]-cc[t12-1];
+ ch[t5]=cc[t11]-cc[t12];
+ ch[t6]=cc[t11]+cc[t12];
+ t5+=ido;
+ t6+=ido;
+ t11+=t10;
+ t12+=t10;
+ }
+ }
+ }
+
+L116:
+ ar1=1.;
+ ai1=0.;
+ t1=0;
+ t9=(t2=ipp2*idl1);
+ t3=(ip-1)*idl1;
+ for(l=1;l<ipph;l++){
+ t1+=idl1;
+ t2-=idl1;
+
+ ar1h=dcp*ar1-dsp*ai1;
+ ai1=dcp*ai1+dsp*ar1;
+ ar1=ar1h;
+ t4=t1;
+ t5=t2;
+ t6=0;
+ t7=idl1;
+ t8=t3;
+ for(ik=0;ik<idl1;ik++){
+ c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
+ c2[t5++]=ai1*ch2[t8++];
+ }
+ dc2=ar1;
+ ds2=ai1;
+ ar2=ar1;
+ ai2=ai1;
+
+ t6=idl1;
+ t7=t9-idl1;
+ for(j=2;j<ipph;j++){
+ t6+=idl1;
+ t7-=idl1;
+ ar2h=dc2*ar2-ds2*ai2;
+ ai2=dc2*ai2+ds2*ar2;
+ ar2=ar2h;
+ t4=t1;
+ t5=t2;
+ t11=t6;
+ t12=t7;
+ for(ik=0;ik<idl1;ik++){
+ c2[t4++]+=ar2*ch2[t11++];
+ c2[t5++]+=ai2*ch2[t12++];
+ }
+ }
+ }
+
+ t1=0;
+ for(j=1;j<ipph;j++){
+ t1+=idl1;
+ t2=t1;
+ for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
+ }
+
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ ch[t3]=c1[t3]-c1[t4];
+ ch[t4]=c1[t3]+c1[t4];
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+
+ if(ido==1)goto L132;
+ if(nbd<l1)goto L128;
+
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(k=0;k<l1;k++){
+ t5=t3;
+ t6=t4;
+ for(i=2;i<ido;i+=2){
+ t5+=2;
+ t6+=2;
+ ch[t5-1]=c1[t5-1]-c1[t6];
+ ch[t6-1]=c1[t5-1]+c1[t6];
+ ch[t5]=c1[t5]+c1[t6-1];
+ ch[t6]=c1[t5]-c1[t6-1];
+ }
+ t3+=ido;
+ t4+=ido;
+ }
+ }
+ goto L132;
+
+ L128:
+ t1=0;
+ t2=ipp2*t0;
+ for(j=1;j<ipph;j++){
+ t1+=t0;
+ t2-=t0;
+ t3=t1;
+ t4=t2;
+ for(i=2;i<ido;i+=2){
+ t3+=2;
+ t4+=2;
+ t5=t3;
+ t6=t4;
+ for(k=0;k<l1;k++){
+ ch[t5-1]=c1[t5-1]-c1[t6];
+ ch[t6-1]=c1[t5-1]+c1[t6];
+ ch[t5]=c1[t5]+c1[t6-1];
+ ch[t6]=c1[t5]-c1[t6-1];
+ t5+=ido;
+ t6+=ido;
+ }
+ }
+ }
+
+L132:
+ if(ido==1)return;
+
+ for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+ t1=0;
+ for(j=1;j<ip;j++){
+ t2=(t1+=t0);
+ for(k=0;k<l1;k++){
+ c1[t2]=ch[t2];
+ t2+=ido;
+ }
+ }
+
+ if(nbd>l1)goto L139;
+
+ is= -ido-1;
+ t1=0;
+ for(j=1;j<ip;j++){
+ is+=ido;
+ t1+=t0;
+ idij=is;
+ t2=t1;
+ for(i=2;i<ido;i+=2){
+ t2+=2;
+ idij+=2;
+ t3=t2;
+ for(k=0;k<l1;k++){
+ c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+ c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+ t3+=ido;
+ }
+ }
+ }
+ return;
+
+ L139:
+ is= -ido-1;
+ t1=0;
+ for(j=1;j<ip;j++){
+ is+=ido;
+ t1+=t0;
+ t2=t1;
+ for(k=0;k<l1;k++){
+ idij=is;
+ t3=t2;
+ for(i=2;i<ido;i+=2){
+ idij+=2;
+ t3+=2;
+ c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+ c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+ }
+ t2+=ido;
+ }
+ }
+}
+
static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
int i,k1,l1,l2;
int na;
goto L115;
L106:
- return; /* silently fail. we only do powers of two in this version */
+ if(ip!=3)goto L109;
+
+ ix2=iw+ido;
+ if(na!=0)
+ dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
+ else
+ dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
+ na=1-na;
+ goto L115;
+
+ L109:
+/* The radix five case can be translated later..... */
+/* if(ip!=5)goto L112;
+
+ ix2=iw+ido;
+ ix3=ix2+ido;
+ ix4=ix3+ido;
+ if(na!=0)
+ dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+ else
+ dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+ na=1-na;
+ goto L115;
+
+ L112:*/
+ if(na!=0)
+ dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+ else
+ dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+ if(ido==1)na=1-na;
L115:
l1=l2;
#ifndef _V_SMFT_H_
#define _V_SMFT_H_
-typedef struct {
- int n;
- double *trigcache;
- int *splitcache;
-} drft_lookup;
+#include "codec.h"
extern void drft_forward(drft_lookup *l,double *data);
extern void drft_backward(drft_lookup *l,double *data);
int i;
vorbis_info *vi=&vb->vd->vi;
double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+ lpc_lookup *vl=&vb->vd->vl[vb->W];
/* get the LSP params back to LPC params. This will be for each
spectral floor curve */
vorbis_lpc_apply(vb->pcm[i],vb->pcmend,vb->lpc[vi->floormap[i]],
vb->amp[vi->floormap[i]],vi->floororder);*/
- /* XXXX Mix back from C+D stereo now */
- {
- int i;
- int n=vb->pcmend;
- 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++){
- L[i]=C[i]+D[i];
- R[i]=C[i]-D[i];
- }
- }
-
for(i=0;i<vb->pcm_channels;i++)
mdct_backward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
_ve_envelope_apply(vb,1);
projects / platform / upstream / libvorbis.git / commitdiff