Commit ff375f51 authored by Monty's avatar Monty
Browse files

Finished LPC and LSP code.  Adding cut down version of full FFT from
old Ogg (the encoder-side needs it).

Monty 19990803

svn path=/trunk/vorbis/; revision=18
parent 36c0bd8c
......@@ -14,7 +14,7 @@
function: PCM data vector blocking, windowing and dis/reassembly
author: Monty <xiphmont@mit.edu>
modifications by: Monty
last modification date: Jul 28 1999
last modification date: Jul 29 1999
Handle windowing, overlap-add, etc of the PCM vectors. This is made
more amusing by Vorbis' current two allowed block sizes (512 and 2048
......@@ -33,10 +33,9 @@
#include "envelope.h"
#include "mdct.h"
/* pcm accumulator and multipliers
examples (not exhaustive):
/* pcm accumulator examples (not exhaustive):
<-------------- lW----------------->
<-------------- lW ---------------->
<--------------- W ---------------->
: .....|..... _______________ |
: .''' | '''_--- | |\ |
......@@ -45,7 +44,6 @@
|<------ Sl ------>| > Sr < |endW
|beginSl |endSl | |endSr
|beginW |endlW |beginSr
mult[0] mult[n]
|< lW >|
......@@ -59,20 +57,18 @@
|beginW | |endlW
mult[0] |beginSl mult[n]
<-------------- lW----------------->
|<-W-->|
: .............. __ | |
: .''' |`/ \ | |
:.....''' |/`...\|...|
:.........................|__||__|___|
|Sl||Sr|endW
| || |endSr
| ||beginSr
| |endSl
<-------------- lW ----------------->
|<--W-->|
: .............. ___ | |
: .''' |`/ \ | |
:.....''' |/`....\|...|
:.........................|___|___|___|
|Sl |Sr |endW
| | |endSr
| |beginSr
| |endSl
|beginSl
|beginW
mult[0]
mult[n]
*/
static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
......
/********************************************************************
* *
* 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: LPC low level routines
author: Monty <monty@xiph.org>
modifications by: Monty
last modification date: Aug 02 1999
********************************************************************/
/* Some of these routines (autocorrelator, LPC coefficient estimator)
are directly derived from and/or modified from code written by
Jutta Degener and Carsten Bormann; thus we include their copyright
below. The entirety of this file is freely redistributable on the
condition that both of these copyright notices are preserved
without modification. */
/* Preserved Copyright: *********************************************/
/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
Technische Universita"t Berlin
Any use of this software is permitted provided that this notice is not
removed and that neither the authors nor the Technische Universita"t
Berlin are deemed to have made any representations as to the
suitability of this software for any purpose nor are held responsible
for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
THIS SOFTWARE.
As a matter of courtesy, the authors request to be informed about uses
this software has found, about bugs in this software, and about any
improvements that may be of general interest.
Berlin, 28.11.1994
Jutta Degener
Carsten Bormann
*********************************************************************/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "smallft.h"
#include "lpc.h"
/* 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
invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
/* Input : n element envelope curve
Output: m lpc coefficients, excitation energy */
double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
double aut[m+1],work[n+n],error;
int i,j;
/* input is a real curve. make it complex-real */
for(i=0;i<n;i++){
work[i*2]=curve[i];
work[i*2+1]=0;
}
n*=2;
fft_backward(n,work,NULL,NULL);
/* The autocorrelation will not be circular. Shift, else we lose
most of the power in the edges. */
for(i=0,j=n/2;i<n/2;){
double temp=work[i];
work[i++]=work[j];
work[j++]=temp;
}
/* autocorrelation, p+1 lag coefficients */
j=m+1;
while(j--){
double d=0;
for(i=j;i<n;i++)d+=work[i]*work[i-j];
aut[j]=d;
}
/* Generate lpc coefficients from autocorr values */
error=aut[0];
if(error==0){
memset(lpc,0,m*sizeof(double));
return 0;
}
for(i=0;i<m;i++){
double r=-aut[i+1];
/* Sum up this iteration's reflection coefficient; note that in
Vorbis we don't save it. If anyone wants to recycle this code
and needs reflection coefficients, save the results of 'r' from
each iteration. */
for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
r/=error;
/* Update LPC coefficients and total error */
lpc[i]=r;
for(j=0;j<i/2;j++){
double tmp=lpc[j];
lpc[j]+=r*lpc[i-1-j];
lpc[i-1-j]+=r*tmp;
}
if(i%2)lpc[j]+=lpc[j]*r;
error*=1.0-r*r;
}
/* we need the error value to know how big an impulse to hit the
filter with later */
return error;
}
/* 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, however */
static double vorbis_lpc_magnitude(double w,double *lpc, int m){
int k;
double real=1,imag=0;
for(k=0;k<m;k++){
real+=lpc[k]*cos((k+1)*w);
imag+=lpc[k]*sin((k+1)*w);
}
return(1/sqrt(real*real+imag*imag));
}
/* generate the whole freq response curve on an LPC IIR filter */
void vorbis_lpc_to_curve(double *curve,int n,double *lpc, double amp,int m){
int i;
double w=1./n*M_PI;
for(i=0;i<n;i++)
curve[i]=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
}
/* find frequency response of LPC filter only at nonsero residue
points and apply the envelope to the residue */
void vorbis_lpc_apply(double *residue,int n,double *lpc, double amp,int m){
int i;
double w=1./n*M_PI;
for(i=0;i<n;i++)
if(residue[i])
residue[i]*=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
}
/********************************************************************
* *
* 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: LPC low level routines
********************************************************************/
#ifndef _V_LPC_H_
#define _V_LPC_H_
extern double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m);
extern void vorbis_lpc_to_curve(double *curve,int n,double *lpc,
double amp,int m);
extern void vorbis_lpc_apply(double *residue,int n,double *lpc,
double amp,int m);
#endif
/********************************************************************
* *
* 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: LSP (also called LSF) conversion routines
author: Monty <monty@xiph.org>
modifications by: Monty
last modification date: Aug 03 1999
The LSP generation code is taken (with minimal modification) from
"On the Computation of the LSP Frequencies" by Joseph Rothweiler
<rothwlr@altavista.net>, available at:
http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html
********************************************************************/
#include <math.h>
#include <string.h>
#include <stdlib.h>
void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
int i,j,m2=m/2;
double O[m2],E[m2];
double A,Ae[m2+1],Ao[m2+1];
double B,Be[m2],Bo[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;
}
lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
Ao[j]=A;
Ae[j]=B;
}
}
static void kw(double *r,int n) {
double s[n/2+1];
double c[n+1];
int i, j, k;
s[0] = 1.0;
s[1] = -2.0;
s[2] = 2.0;
for(i=3;i<=n/2;i++) s[i] = s[i-2];
for(k=0;k<=n;k++) {
c[k] = r[k];
j = 1;
for(i=k+2;i<=n;i+=2) {
c[k] += s[j]*r[i];
s[j] -= s[j-1];
j++;
}
}
for(k=0;k<=n;k++) r[k] = c[k];
}
static int comp(const void *a,const void *b){
return(*(double *)a<*(double *)b);
}
/* CACM algorithm 283. */
static void cacm283(double *a,int ord,double *r){
int i, k;
double val, p, delta, error;
double rooti;
for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
for(error=1 ; error > 1.e-12; ) {
error = 0;
for( i=0; i<ord; i++) { /* Update each point. */
rooti = r[i];
val = a[ord];
p = a[ord];
for(k=ord-1; k>= 0; k--) {
val = val * rooti + a[k];
if (k != i) p *= rooti - r[k];
}
delta = val/p;
r[i] -= delta;
error += delta*delta;
}
}
/* Replaced the original bubble sort with a real sort. With your
help, we can eliminate the bubble sort in our lifetime. --Monty */
qsort(r,ord,sizeof(double),comp);
}
/* Convert lpc coefficients to lsp coefficients */
void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
int order2=m/2;
double g1[order2+1], g2[order2+1];
double g1r[order2+1], g2r[order2+1];
int i;
/* Compute the lengths of the x polynomials. */
/* Compute the first half of K & R F1 & F2 polynomials. */
/* Compute half of the symmetric and antisymmetric polynomials. */
/* Remove the roots at +1 and -1. */
g1[order2] = 1.0;
for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
g2[order2] = 1.0;
for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
for(i=0; i<order2;i++) g1[order2-i-1] -= g1[order2-i];
for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
/* Convert into polynomials in cos(alpha) */
kw(g1,order2);
kw(g2,order2);
/* Find the roots of the 2 even polynomials.*/
cacm283(g1,order2,g1r);
cacm283(g2,order2,g2r);
for(i=0;i<m;i+=2){
lsp[i] = acos(g1r[i/2]*.5);
lsp[i+1] = acos(g2r[i/2]*.5);
}
}
/********************************************************************
* *
* 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/ *
* *
********************************************************************/
#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);
#endif
/******************************************************************
* CopyPolicy: GNU Public License 2 applies
* Copyright (C) 1998 Monty xiphmont@mit.edu, monty@xiph.org
*
* FFT implementation from OggSquish, minus cosine transforms,
* minus all but radix 2/4 case. In Vorbis we only need this
* cut-down version (used by the encoder, not decoder).
*
* To do more than just power-of-two sized vectors, see the full
* version I wrote for NetLib.
*
******************************************************************/
#include <stdlib.h>
#include <math.h>
#include "smallft.h"
static void drfti1(int n, double *wa, int *ifac){
static int ntryh[4] = { 4,2,3,5 };
static double tpi = 6.28318530717958647692528676655900577;
double arg,argh,argld,fi;
int ntry=0,i,j=-1;
int k1, l1, l2, ib;
int ld, ii, ip, is, nq, nr;
int ido, ipm, nfm1;
int nl=n;
int nf=0;
L101:
j++;
if (j < 4)
ntry=ntryh[j];
else
ntry+=2;
L104:
nq=nl/ntry;
nr=nl-ntry*nq;
if (nr!=0) goto L101;
nf++;
ifac[nf+1]=ntry;
nl=nq;
if(ntry!=2)goto L107;
if(nf==1)goto L107;
for (i=1;i<nf;i++){
ib=nf-i+1;
ifac[ib+1]=ifac[ib];
}
ifac[2] = 2;
L107:
if(nl!=1)goto L104;
ifac[0]=n;
ifac[1]=nf;
argh=tpi/n;
is=0;
nfm1=nf-1;
l1=1;
if(nfm1==0)return;
for (k1=0;k1<nfm1;k1++){
ip=ifac[k1+2];
ld=0;
l2=l1*ip;
ido=n/l2;
ipm=ip-1;
for (j=0;j<ipm;j++){
ld+=l1;
i=is;
argld=(double)ld*argh;
fi=0.;
for (ii=2;ii<ido;ii+=2){
fi+=1.;
arg=fi*argld;
wa[i++]=cos(arg);
wa[i++]=sin(arg);
}
is+=ido;
}
l1=l2;
}
}
static void fdrffti(int n, double *wsave, int *ifac){
if (n == 1) return;
drfti1(n, wsave+n, ifac);
}
static void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){
int i,k;
double ti2,tr2;
int t0,t1,t2,t3,t4,t5,t6;
t1=0;
t0=(t2=l1*ido);
t3=ido<<1;
for(k=0;k<l1;k++){
ch[t1<<1]=cc[t1]+cc[t2];
ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
t1+=ido;
t2+=ido;
}
if(ido<2)return;
if(ido==2)goto L105;
t1=0;
t2=t0;
for(k=0;k<l1;k++){
t3=t2;
t4=(t1<<1)+(ido<<1);
t5=t1;
t6=t1+t1;
for(i=2;i<ido;i+=2){
t3+=2;
t4-=2;
t5+=2;
t6+=2;
tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
ch[t6]=cc[t5]+ti2;
ch[t4]=ti2-cc[t5];
ch[t6-1]=cc[t5-1]+tr2;
ch[t4-1]=cc[t5-1]-tr2;
}
t1+=ido;
t2+=ido;
}
if(ido%2==1)return;
L105:
t3=(t2=(t1=ido)-1);
t2+=t0;
for(k=0;k<l1;k++){
ch[t1]=-cc[t2];
ch[t1-1]=cc[t3];
t1+=ido<<1;
t2+=ido;
t3+=ido;
}
}
static void dradf4(int ido,int l1,double *cc,double *ch,double *wa1,
double *wa2,double *wa3){
static double hsqt2 = .70710678118654752440084436210485;
int i,k,t0,t1,t2,t3,t4,t5,t6;
double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
t0=l1*ido;
t1=t0;
t4=t1<<1;
t2=t1+(t1<<1);
t3=0;
for(k=0;k<l1;k++){
tr1=cc[t1]+cc[t2];
tr2=cc[t3]+cc[t4];
ch[t5=t3<<2]=tr1+tr2;
ch[(ido<<2)+t5-1]=tr2-tr1;
ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
ch[t5]=cc[t2]-cc[t1];
t1+=ido;
t2+=ido;
t3+=ido;
t4+=ido;
}
if(ido<2)return;
if(ido==2)goto L105;
t1=0;