vrshoot

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

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