dbf-halloween2015

annotate libs/vorbis/lpc.c @ 1:c3f5c32cb210

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