Commit e74b2b02 authored by Monty's avatar Monty
Browse files

All new LSP->freq envelope curve computation code.

It turns out that LSP->LPC using the impulse response algorithm is
*very* sensitive to noise, and doubles really are necessary.
Unfortunate, that.

Reimplmented the code with a direct LSP->curve computation, skipping
the LPC intermediary step.  This also eliminates any need for the LPC
or iFFT code in decode/synthesis.

Monty

svn path=/trunk/vorbis/; revision=597
parent 45e1cc5e
......@@ -12,7 +12,7 @@
********************************************************************
function: floor backend 0 implementation
last mod: $Id: floor0.c,v 1.20 2000/08/15 09:09:42 xiphmont Exp $
last mod: $Id: floor0.c,v 1.21 2000/08/19 11:46:28 xiphmont Exp $
********************************************************************/
......@@ -41,6 +41,7 @@ typedef struct {
vorbis_info_floor0 *vi;
lpc_lookup lpclook;
double *lsp_look;
} vorbis_look_floor0;
......@@ -82,6 +83,7 @@ static void free_look(vorbis_look_floor *i){
vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
if(i){
if(look->linearmap)free(look->linearmap);
if(look->lsp_look)free(look->lsp_look);
lpc_clear(&look->lpclook);
memset(look,0,sizeof(vorbis_look_floor0));
free(look);
......@@ -145,7 +147,9 @@ static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
look->n=vi->blocksizes[mi->blockflag]/2;
look->ln=info->barkmap;
look->vi=info;
lpc_init(&look->lpclook,look->ln,look->m);
if(vd->analysisp)
lpc_init(&look->lpclook,look->ln,look->m);
/* we choose a scaling constant so that:
floor(bark(rate/2-1)*C)=mapped-1
......@@ -166,6 +170,10 @@ static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
look->linearmap[j]=val;
}
look->lsp_look=malloc(look->ln*sizeof(double));
for(j=0;j<look->ln;j++)
look->lsp_look[j]=2*cos(M_PI/look->ln*j);
return look;
}
......@@ -235,31 +243,17 @@ double _curve_to_lpc(double *curve,double *lpc,
/* generate the whole freq response curve of an LPC IIR filter */
void _lpc_to_curve(double *curve,double *lpc,double amp,
void _lsp_to_curve(double *curve,double *lsp,double amp,
vorbis_look_floor0 *l,char *name,long frameno){
/* l->m+1 must be less than l->ln, but guard in case we get a bad stream */
double *lcurve=alloca(sizeof(double)*max(l->ln*2,l->m*2+2));
double *lcurve=alloca(sizeof(double)*l->ln);
int i;
if(amp==0){
memset(curve,0,sizeof(double)*l->n);
return;
}
vorbis_lpc_to_curve(lcurve,lpc,amp,&(l->lpclook));
#if 0
{ /******************/
FILE *of;
char buffer[80];
int i;
sprintf(buffer,"%s_%d.m",name,frameno);
of=fopen(buffer,"w");
for(i=0;i<l->ln;i++)
fprintf(of,"%g\n",lcurve[i]);
fclose(of);
}
#endif
vorbis_lsp_to_curve(lcurve,l->ln,lsp,l->m,amp,l->lsp_look);
for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
......@@ -335,7 +329,7 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
#ifdef ANALYSIS
if(vb->W==0){fprintf(stderr,"%d ",seq);}
vorbis_lsp_to_lpc(out,work,look->m);
_lpc_to_curve(work,work,amp,look,"Ffloor",seq);
_lsp_to_curve(work,work,amp,look,"Ffloor",seq);
for(j=0;j<look->n;j++)work[j]-=info->ampdB;
_analysis_output("rawfloor",seq,work,look->n,0,0);
{
......@@ -394,8 +388,7 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
#endif
/* take the coefficients back to a spectral envelope curve */
vorbis_lsp_to_lpc(work,out,look->m);
_lpc_to_curve(out,out,amp,look,"Ffloor",seq++);
_lsp_to_curve(out,work,amp,look,"Ffloor",seq++);
for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
return(1);
}
......@@ -430,8 +423,7 @@ static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
}
/* take the coefficients back to a spectral envelope curve */
vorbis_lsp_to_lpc(out,out,look->m);
_lpc_to_curve(out,out,amp,look,"",0);
_lsp_to_curve(out,out,amp,look,"",0);
for(j=0;j<look->n;j++)out[j]=fromdB(out[j]-info->ampdB);
return(1);
......
......@@ -12,7 +12,7 @@
********************************************************************
function: LPC low level routines
last mod: $Id: lpc.c,v 1.23 2000/08/15 09:09:43 xiphmont Exp $
last mod: $Id: lpc.c,v 1.24 2000/08/19 11:46:28 xiphmont Exp $
********************************************************************/
......@@ -165,127 +165,3 @@ void lpc_clear(lpc_lookup *l){
}
}
/* One can do this the long way by generating the transfer function in
the time domain and taking the forward FFT of the result. The
results from direct calculation are cleaner and faster.
This version does a linear curve generation and then later
interpolates the log curve from the linear curve. */
void vorbis_lpc_to_curve(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);
}
drft_backward(&l->fft,curve); /* reappropriated ;-) */
{
int l2=l->ln*2;
double unit=1./amp;
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]);
double a = real + unit;
curve[i] = 1.0 / FAST_HYPOT(a, imag);
}
}
}
/* subtract or add an lpc filter to data. */
void vorbis_lpc_filter(double *coeff,double *prime,int m,
double *data,long n,double amp){
/* in: coeff[0...m-1] LPC coefficients
prime[0...m-1] initial values
data[0...n-1] data samples
out: data[0...n-1] residuals from LPC prediction */
long i,j;
double *work=alloca(sizeof(double)*(m+n));
double y;
if(!prime)
for(i=0;i<m;i++)
work[i]=0;
else
for(i=0;i<m;i++)
work[i]=prime[i];
for(i=0;i<n;i++){
y=0;
for(j=0;j<m;j++)
y-=work[i+j]*coeff[m-j-1];
data[i]=work[i+m]=data[i]+y;
}
}
void vorbis_lpc_residue(double *coeff,double *prime,int m,
double *data,long n){
/* in: coeff[0...m-1] LPC coefficients
prime[0...m-1] initial values
data[0...n-1] data samples
out: data[0...n-1] residuals from LPC prediction */
long i,j;
double *work=alloca(sizeof(double)*(m+n));
double y;
if(!prime)
for(i=0;i<m;i++)
work[i]=0;
else
for(i=0;i<m;i++)
work[i]=prime[i];
for(i=0;i<n;i++){
y=0;
for(j=0;j<m;j++)
y-=work[i+j]*coeff[m-j-1];
work[i+m]=data[i];
data[i]-=y;
}
}
void vorbis_lpc_predict(double *coeff,double *prime,int m,
double *data,long n){
/* in: coeff[0...m-1] LPC coefficients
prime[0...m-1] initial values (allocated size of n+m-1)
data[0...n-1] residuals from LPC prediction
out: data[0...n-1] data samples */
long i,j,o,p;
double y;
double *work=alloca(sizeof(double)*(m+n));
if(!prime)
for(i=0;i<m;i++)
work[i]=0.;
else
for(i=0;i<m;i++)
work[i]=prime[i];
for(i=0;i<n;i++){
y=data[i];
o=i;
p=m;
for(j=0;j<m;j++)
y-=work[o++]*coeff[--p];
data[i]=work[o]=y;
}
}
......@@ -12,7 +12,7 @@
********************************************************************
function: LPC low level routines
last mod: $Id: lpc.h,v 1.11 2000/07/12 09:36:18 xiphmont Exp $
last mod: $Id: lpc.h,v 1.12 2000/08/19 11:46:28 xiphmont Exp $
********************************************************************/
......@@ -37,15 +37,5 @@ extern void lpc_clear(lpc_lookup *l);
/* simple linear scale LPC code */
extern double vorbis_lpc_from_data(double *data,double *lpc,int n,int m);
extern double vorbis_lpc_from_curve(double *curve,double *lpc,lpc_lookup *l);
extern void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,
lpc_lookup *l);
/* standard lpc stuff */
extern void vorbis_lpc_filter(double *coeff,double *prime,int m,
double *data,long n,double amp);
extern void vorbis_lpc_residue(double *coeff,double *prime,int m,
double *data,long n);
extern void vorbis_lpc_predict(double *coeff,double *prime,int m,
double *data,long n);
#endif
......@@ -12,7 +12,7 @@
********************************************************************
function: LSP (also called LSF) conversion routines
last mod: $Id: lsp.c,v 1.8 2000/05/08 20:49:49 xiphmont Exp $
last mod: $Id: lsp.c,v 1.9 2000/08/19 11:46:28 xiphmont Exp $
The LSP generation code is taken (with minimal modification) from
"On the Computation of the LSP Frequencies" by Joseph Rothweiler
......@@ -39,51 +39,21 @@
#include "os.h"
#include "misc.h"
void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
int i,j,m2=m/2;
double *O=alloca(sizeof(double)*m2);
double *E=alloca(sizeof(double)*m2);
double A;
double *Ae=alloca(sizeof(double)*(m2+1));
double *Ao=alloca(sizeof(double)*(m2+1));
double B;
double *Be=alloca(sizeof(double)*(m2));
double *Bo=alloca(sizeof(double)*(m2));
double temp;
/* even/odd roots setup */
for(i=0;i<m2;i++){
O[i]=-2.*cos(lsp[i*2]);
E[i]=-2.*cos(lsp[i*2+1]);
}
/* set up impulse response */
for(j=0;j<m2;j++){
Ae[j]=0.;
Ao[j]=1.;
Be[j]=0.;
Bo[j]=1.;
}
Ao[j]=1.;
Ae[j]=1.;
/* run impulse response */
for(i=1;i<m+1;i++){
A=B=0.;
for(j=0;j<m2;j++){
temp=O[j]*Ao[j]+Ae[j];
Ae[j]=Ao[j];
Ao[j]=A;
A+=temp;
temp=E[j]*Bo[j]+Be[j];
Be[j]=Bo[j];
Bo[j]=B;
B+=temp;
void vorbis_lsp_to_curve(double *curve,int n,double *lsp,int m,double amp,
double *w){
int i,j,k;
double *coslsp=alloca(m*sizeof(double));
for(i=0;i<m;i++)coslsp[i]=2*cos(lsp[i]);
for(k=0;k<n;k++){
double p=.70710678118654752440;
double q=.70710678118654752440;
for(j=0;j<m;){
p*= *w-coslsp[j++];
q*= *w-coslsp[j++];
}
lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
Ao[j]=A;
Ae[j]=B;
curve[k]=amp/sqrt(p*p*(1.+ *w*.5)+q*q*(1.- *w*.5));
w++;
}
}
......
......@@ -12,7 +12,7 @@
********************************************************************
function: LSP (also called LSF) conversion routines
last mod: $Id: lsp.h,v 1.3 1999/12/30 07:26:43 xiphmont Exp $
last mod: $Id: lsp.h,v 1.4 2000/08/19 11:46:28 xiphmont Exp $
********************************************************************/
......@@ -20,7 +20,9 @@
#ifndef _V_LSP_H_
#define _V_LSP_H_
extern void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m);
extern void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m);
extern void vorbis_lsp_to_curve(double *curve,int n,
double *lsp,int m,double amp,
double *w);
#endif
......@@ -13,7 +13,7 @@
function: simple utility that runs audio through the psychoacoustics
without encoding
last mod: $Id: psytune.c,v 1.5 2000/08/15 09:09:43 xiphmont Exp $
last mod: $Id: psytune.c,v 1.6 2000/08/19 11:46:28 xiphmont Exp $
********************************************************************/
......@@ -151,7 +151,7 @@ typedef struct {
extern double _curve_to_lpc(double *curve,double *lpc,vorbis_look_floor0 *l,
long frameno);
extern void _lpc_to_curve(double *curve,double *lpc,double amp,
extern void _lsp_to_curve(double *curve,double *lpc,double amp,
vorbis_look_floor0 *l,char *name,long frameno);
long frameno=0;
......@@ -299,7 +299,8 @@ int main(int argc,char *argv[]){
for(j=0;j<framesize/2;j++)floor[j]=todB(floor[j])+150;
amp=_curve_to_lpc(floor,lpc,&floorlook,frameno);
_lpc_to_curve(floor,lpc,sqrt(amp),&floorlook,"Ffloor",frameno);
vorbis_lpc_to_lsp(lpc,lpc,order);
_lsp_to_curve(floor,lpc,sqrt(amp),&floorlook,"Ffloor",frameno);
for(j=0;j<framesize/2;j++)floor[j]=fromdB(floor[j]-150);
analysis("floor",frameno,floor,framesize/2,1,1);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment