18e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/******************************************************************** 28e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * * 38e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * 48e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * 58e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * 68e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * 78e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * * 88e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 * 98e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * by the Xiph.Org Foundation http://www.xiph.org/ * 108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * * 118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ******************************************************************** 128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels function: LPC low level routines 148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels last mod: $Id: lpc.c 16227 2009-07-08 06:58:46Z xiphmont $ 158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ********************************************************************/ 178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Some of these routines (autocorrelator, LPC coefficient estimator) 198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels are derived from code written by Jutta Degener and Carsten Bormann; 208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels thus we include their copyright below. The entirety of this file 218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels is freely redistributable on the condition that both of these 228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels copyright notices are preserved without modification. */ 238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Preserved Copyright: *********************************************/ 258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann, 278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsTechnische Universita"t Berlin 288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsAny use of this software is permitted provided that this notice is not 308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsremoved and that neither the authors nor the Technische Universita"t 318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsBerlin are deemed to have made any representations as to the 328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelssuitability of this software for any purpose nor are held responsible 338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsfor any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR 348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsTHIS SOFTWARE. 358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsAs a matter of courtesy, the authors request to be informed about uses 378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsthis software has found, about bugs in this software, and about any 388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsimprovements that may be of general interest. 398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsBerlin, 28.11.1994 418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsJutta Degener 428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas EckelsCarsten Bormann 438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels*********************************************************************/ 458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <stdlib.h> 478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <string.h> 488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <math.h> 498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "os.h" 508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "smallft.h" 518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "lpc.h" 528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "scales.h" 538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "misc.h" 548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Autocorrelation LPC coeff generation algorithm invented by 568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels N. Levinson in 1947, modified by J. Durbin in 1959. */ 578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Input : n elements of time doamin data 598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels Output: m lpc coefficients, excitation energy */ 608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsfloat vorbis_lpc_from_data(float *data,float *lpci,int n,int m){ 628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double *aut=alloca(sizeof(*aut)*(m+1)); 638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double *lpc=alloca(sizeof(*lpc)*(m)); 648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double error; 658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double epsilon; 668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int i,j; 678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* autocorrelation, p+1 lag coefficients */ 698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels j=m+1; 708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels while(j--){ 718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double d=0; /* double needed for accumulator depth */ 728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=j;i<n;i++)d+=(double)data[i]*data[i-j]; 738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels aut[j]=d; 748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* Generate lpc coefficients from autocorr values */ 778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* set our noise floor to about -100dB */ 798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels error=aut[0] * (1. + 1e-10); 808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels epsilon=1e-9*aut[0]+1e-10; 818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<m;i++){ 838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double r= -aut[i+1]; 848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(error<epsilon){ 868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memset(lpc+i,0,(m-i)*sizeof(*lpc)); 878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels goto done; 888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* Sum up this iteration's reflection coefficient; note that in 918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels Vorbis we don't save it. If anyone wants to recycle this code 928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels and needs reflection coefficients, save the results of 'r' from 938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels each iteration. */ 948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<i;j++)r-=lpc[j]*aut[i-j]; 968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels r/=error; 978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* Update LPC coefficients and total error */ 998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels lpc[i]=r; 1018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<i/2;j++){ 1028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double tmp=lpc[j]; 1038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels lpc[j]+=r*lpc[i-1-j]; 1058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels lpc[i-1-j]+=r*tmp; 1068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(i&1)lpc[j]+=lpc[j]*r; 1088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels error*=1.-r*r; 1108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels done: 1148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* slightly damp the filter */ 1168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels { 1178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double g = .99; 1188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels double damp = g; 1198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<m;j++){ 1208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels lpc[j]*=damp; 1218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels damp*=g; 1228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<m;j++)lpci[j]=(float)lpc[j]; 1268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* we need the error value to know how big an impulse to hit the 1288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels filter with later */ 1298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels return error; 1318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 1328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vorbis_lpc_predict(float *coeff,float *prime,int m, 1348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *data,long n){ 1358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* in: coeff[0...m-1] LPC coefficients 1378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels prime[0...m-1] initial values (allocated size of n+m-1) 1388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels out: data[0...n-1] data samples */ 1398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long i,j,o,p; 1418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float y; 1428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *work=alloca(sizeof(*work)*(m+n)); 1438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(!prime) 1458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<m;i++) 1468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels work[i]=0.f; 1478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels else 1488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<m;i++) 1498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels work[i]=prime[i]; 1508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<n;i++){ 1528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels y=0; 1538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels o=i; 1548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels p=m; 1558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<m;j++) 1568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels y-=work[o++]*coeff[--p]; 1578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels data[i]=work[o]=y; 1598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 161