1cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com/******************************************************************** 2cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * * 3cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * 4cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * 5cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * 6cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * 7cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * * 8cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 * 9cf2cfa174ca878c144e17e9fc60ca8e9070d7dededisonn@google.com * by the Xiph.Org Foundation http://www.xiph.org/ * 10ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com * * 11e91260c3e75d9d26e391be6823ea37373b27d1ffedisonn@google.com ******************************************************************** 12ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 13ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com function: LPC low level routines 142af2ad9cc0b2c7d911aed2e8d2ac77c0b7d3b5dfedisonn@google.com last mod: $Id: lpc.c 16227 2009-07-08 06:58:46Z xiphmont $ 152af2ad9cc0b2c7d911aed2e8d2ac77c0b7d3b5dfedisonn@google.com 162af2ad9cc0b2c7d911aed2e8d2ac77c0b7d3b5dfedisonn@google.com ********************************************************************/ 172af2ad9cc0b2c7d911aed2e8d2ac77c0b7d3b5dfedisonn@google.com 182af2ad9cc0b2c7d911aed2e8d2ac77c0b7d3b5dfedisonn@google.com/* Some of these routines (autocorrelator, LPC coefficient estimator) 19d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.com are derived from code written by Jutta Degener and Carsten Bormann; 20d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.com thus we include their copyright below. The entirety of this file 21d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.com is freely redistributable on the condition that both of these 22d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.com copyright notices are preserved without modification. */ 23d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.com 242af2ad9cc0b2c7d911aed2e8d2ac77c0b7d3b5dfedisonn@google.com/* Preserved Copyright: *********************************************/ 251f2f338e23789f3eef168dcbd8171a28820ba6c1robertphillips@google.com 26ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann, 27ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comTechnische Universita"t Berlin 28ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 291f2f338e23789f3eef168dcbd8171a28820ba6c1robertphillips@google.comAny use of this software is permitted provided that this notice is not 30ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comremoved and that neither the authors nor the Technische Universita"t 31ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comBerlin are deemed to have made any representations as to the 32ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comsuitability of this software for any purpose nor are held responsible 331f2f338e23789f3eef168dcbd8171a28820ba6c1robertphillips@google.comfor any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR 34ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comTHIS SOFTWARE. 35ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 36ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comAs a matter of courtesy, the authors request to be informed about uses 37ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comthis software has found, about bugs in this software, and about any 38d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.comimprovements that may be of general interest. 39d03c2c732e2e914f8429cfc8b077a7b9b853dd8eedisonn@google.com 40ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comBerlin, 28.11.1994 41ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comJutta Degener 42ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comCarsten Bormann 43ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 44ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com*********************************************************************/ 45ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 46ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include <stdlib.h> 47ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include <string.h> 48ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include <math.h> 49ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include "os.h" 50ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include "smallft.h" 51ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include "lpc.h" 52ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include "scales.h" 53ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com#include "misc.h" 54ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 55ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com/* Autocorrelation LPC coeff generation algorithm invented by 56ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com N. Levinson in 1947, modified by J. Durbin in 1959. */ 57ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 58ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com/* Input : n elements of time doamin data 59ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com Output: m lpc coefficients, excitation energy */ 60ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 61ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comfloat vorbis_lpc_from_data(float *data,float *lpci,int n,int m){ 62ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double *aut=alloca(sizeof(*aut)*(m+1)); 63ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double *lpc=alloca(sizeof(*lpc)*(m)); 64ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double error; 65ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double epsilon; 66ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com int i,j; 67ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 68ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* autocorrelation, p+1 lag coefficients */ 69ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com j=m+1; 70ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com while(j--){ 71ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double d=0; /* double needed for accumulator depth */ 72ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(i=j;i<n;i++)d+=(double)data[i]*data[i-j]; 73ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com aut[j]=d; 74ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 75ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 76ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* Generate lpc coefficients from autocorr values */ 77ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 78ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* set our noise floor to about -100dB */ 79ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com error=aut[0] * (1. + 1e-10); 80ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com epsilon=1e-9*aut[0]+1e-10; 81ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 82ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(i=0;i<m;i++){ 83ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double r= -aut[i+1]; 84ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 85ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com if(error<epsilon){ 86ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com memset(lpc+i,0,(m-i)*sizeof(*lpc)); 87ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com goto done; 88ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 89ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 90ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* Sum up this iteration's reflection coefficient; note that in 91ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com Vorbis we don't save it. If anyone wants to recycle this code 92ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com and needs reflection coefficients, save the results of 'r' from 93ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com each iteration. */ 94ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 95ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(j=0;j<i;j++)r-=lpc[j]*aut[i-j]; 96ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com r/=error; 97ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 98ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* Update LPC coefficients and total error */ 9950bbdb4f3edc5005a71c438651732f53c00a5331edisonn@google.com 100ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com lpc[i]=r; 101ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(j=0;j<i/2;j++){ 102ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double tmp=lpc[j]; 103ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 104ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com lpc[j]+=r*lpc[i-1-j]; 105ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com lpc[i-1-j]+=r*tmp; 106ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 107ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com if(i&1)lpc[j]+=lpc[j]*r; 108ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 109ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com error*=1.-r*r; 110ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 111ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 1122c8177767a4b3b6c27c6ac071c8619b557472521edisonn@google.com 1132c8177767a4b3b6c27c6ac071c8619b557472521edisonn@google.com done: 114ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 1152c8177767a4b3b6c27c6ac071c8619b557472521edisonn@google.com /* slightly damp the filter */ 116ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com { 117ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double g = .99; 118ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com double damp = g; 119ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(j=0;j<m;j++){ 120ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com lpc[j]*=damp; 121ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com damp*=g; 122ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 123ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 124ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 125ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(j=0;j<m;j++)lpci[j]=(float)lpc[j]; 126ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 127ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* we need the error value to know how big an impulse to hit the 128ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com filter with later */ 129ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 130ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com return error; 131ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com} 132ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 133ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.comvoid vorbis_lpc_predict(float *coeff,float *prime,int m, 134ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com float *data,long n){ 135ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 136ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com /* in: coeff[0...m-1] LPC coefficients 137ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com prime[0...m-1] initial values (allocated size of n+m-1) 138ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com out: data[0...n-1] data samples */ 139ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 140ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com long i,j,o,p; 141ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com float y; 142ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com float *work=alloca(sizeof(*work)*(m+n)); 143ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 144ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com if(!prime) 145ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(i=0;i<m;i++) 146ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com work[i]=0.f; 147ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com else 148e50d9a1fcd9c4298079ff54f9a40c9708d30f8c6edisonn@google.com for(i=0;i<m;i++) 149e50d9a1fcd9c4298079ff54f9a40c9708d30f8c6edisonn@google.com work[i]=prime[i]; 150ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 151ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(i=0;i<n;i++){ 152ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com y=0; 1531f2f338e23789f3eef168dcbd8171a28820ba6c1robertphillips@google.com o=i; 154ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com p=m; 155ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com for(j=0;j<m;j++) 156ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com y-=work[o++]*coeff[--p]; 157ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com 158ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com data[i]=work[o]=y; 159ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com } 160ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com} 161ac03d91ee03599eab946a8ad25e33f9fc5f3166eedisonn@google.com