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