lpc.c 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260
  1. /********************************************************************
  2. * *
  3. * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
  4. * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
  5. * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
  6. * PLEASE READ THESE TERMS DISTRIBUTING. *
  7. * *
  8. * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
  9. * by Monty <monty@xiph.org> and The XIPHOPHORUS Company *
  10. * http://www.xiph.org/ *
  11. * *
  12. ********************************************************************
  13. function: LPC low level routines
  14. last mod: $Id: lpc.c,v 1.20.2.1 2000/05/24 21:17:01 xiphmont Exp $
  15. ********************************************************************/
  16. /* Some of these routines (autocorrelator, LPC coefficient estimator)
  17. are derived from code written by Jutta Degener and Carsten Bormann;
  18. thus we include their copyright below. The entirety of this file
  19. is freely redistributable on the condition that both of these
  20. copyright notices are preserved without modification. */
  21. /* Preserved Copyright: *********************************************/
  22. /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  23. Technische Universita"t Berlin
  24. Any use of this software is permitted provided that this notice is not
  25. removed and that neither the authors nor the Technische Universita"t
  26. Berlin are deemed to have made any representations as to the
  27. suitability of this software for any purpose nor are held responsible
  28. for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
  29. THIS SOFTWARE.
  30. As a matter of courtesy, the authors request to be informed about uses
  31. this software has found, about bugs in this software, and about any
  32. improvements that may be of general interest.
  33. Berlin, 28.11.1994
  34. Jutta Degener
  35. Carsten Bormann
  36. *********************************************************************/
  37. #include <stdlib.h>
  38. #include <string.h>
  39. #include <math.h>
  40. #include "os.h"
  41. #include "smallft.h"
  42. #include "lpc.h"
  43. #include "scales.h"
  44. #include "misc.h"
  45. /* Autocorrelation LPC coeff generation algorithm invented by
  46. N. Levinson in 1947, modified by J. Durbin in 1959. */
  47. /* Input : n elements of time doamin data
  48. Output: m lpc coefficients, excitation energy */
  49. double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
  50. double *aut=alloca(sizeof(double)*(m+1));
  51. double error;
  52. int i,j;
  53. /* autocorrelation, p+1 lag coefficients */
  54. j=m+1;
  55. while(j--){
  56. double d=0;
  57. for(i=j;i<n;i++)d+=data[i]*data[i-j];
  58. aut[j]=d;
  59. }
  60. /* Generate lpc coefficients from autocorr values */
  61. error=aut[0];
  62. if(error==0){
  63. memset(lpc,0,m*sizeof(double));
  64. return 0;
  65. }
  66. for(i=0;i<m;i++){
  67. double r=-aut[i+1];
  68. /* Sum up this iteration's reflection coefficient; note that in
  69. Vorbis we don't save it. If anyone wants to recycle this code
  70. and needs reflection coefficients, save the results of 'r' from
  71. each iteration. */
  72. for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
  73. r/=error;
  74. /* Update LPC coefficients and total error */
  75. lpc[i]=r;
  76. for(j=0;j<i/2;j++){
  77. double tmp=lpc[j];
  78. lpc[j]+=r*lpc[i-1-j];
  79. lpc[i-1-j]+=r*tmp;
  80. }
  81. if(i%2)lpc[j]+=lpc[j]*r;
  82. error*=1.0-r*r;
  83. }
  84. /* we need the error value to know how big an impulse to hit the
  85. filter with later */
  86. return error;
  87. }
  88. /* Input : n element envelope spectral curve
  89. Output: m lpc coefficients, excitation energy */
  90. double vorbis_lpc_from_curve(double *curve,double *lpc,lpc_lookup *l){
  91. int n=l->ln;
  92. int m=l->m;
  93. double *work=alloca(sizeof(double)*(n+n));
  94. double fscale=.5/n;
  95. int i,j;
  96. /* input is a real curve. make it complex-real */
  97. /* This mixes phase, but the LPC generation doesn't care. */
  98. for(i=0;i<n;i++){
  99. work[i*2]=curve[i]*fscale;
  100. work[i*2+1]=0;
  101. }
  102. work[n*2-1]=curve[n-1]*fscale;
  103. n*=2;
  104. drft_backward(&l->fft,work);
  105. /* The autocorrelation will not be circular. Shift, else we lose
  106. most of the power in the edges. */
  107. for(i=0,j=n/2;i<n/2;){
  108. double temp=work[i];
  109. work[i++]=work[j];
  110. work[j++]=temp;
  111. }
  112. return(vorbis_lpc_from_data(work,lpc,n,m));
  113. }
  114. void lpc_init(lpc_lookup *l,long mapped, int m){
  115. memset(l,0,sizeof(lpc_lookup));
  116. l->ln=mapped;
  117. l->m=m;
  118. /* we cheat decoding the LPC spectrum via FFTs */
  119. drft_init(&l->fft,mapped*2);
  120. }
  121. void lpc_clear(lpc_lookup *l){
  122. if(l){
  123. drft_clear(&l->fft);
  124. }
  125. }
  126. /* One can do this the long way by generating the transfer function in
  127. the time domain and taking the forward FFT of the result. The
  128. results from direct calculation are cleaner and faster.
  129. This version does a linear curve generation and then later
  130. interpolates the log curve from the linear curve. */
  131. void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,
  132. lpc_lookup *l){
  133. int i;
  134. memset(curve,0,sizeof(double)*l->ln*2);
  135. if(amp==0)return;
  136. for(i=0;i<l->m;i++){
  137. curve[i*2+1]=lpc[i]/(4*amp);
  138. curve[i*2+2]=-lpc[i]/(4*amp);
  139. }
  140. drft_backward(&l->fft,curve); /* reappropriated ;-) */
  141. {
  142. int l2=l->ln*2;
  143. double unit=1./amp;
  144. curve[0]=(1./(curve[0]*2+unit));
  145. for(i=1;i<l->ln;i++){
  146. double real=(curve[i]+curve[l2-i]);
  147. double imag=(curve[i]-curve[l2-i]);
  148. double a = real + unit;
  149. curve[i] = 1.0 / FAST_HYPOT(a, imag);
  150. }
  151. }
  152. }
  153. /* subtract or add an lpc filter to data. Vorbis doesn't actually use this. */
  154. void vorbis_lpc_residue(double *coeff,double *prime,int m,
  155. double *data,long n){
  156. /* in: coeff[0...m-1] LPC coefficients
  157. prime[0...m-1] initial values
  158. data[0...n-1] data samples
  159. out: data[0...n-1] residuals from LPC prediction */
  160. long i,j;
  161. double *work=alloca(sizeof(double)*(m+n));
  162. double y;
  163. if(!prime)
  164. for(i=0;i<m;i++)
  165. work[i]=0;
  166. else
  167. for(i=0;i<m;i++)
  168. work[i]=prime[i];
  169. for(i=0;i<n;i++){
  170. y=0;
  171. for(j=0;j<m;j++)
  172. y-=work[i+j]*coeff[m-j-1];
  173. work[i+m]=data[i];
  174. data[i]-=y;
  175. }
  176. }
  177. void vorbis_lpc_predict(double *coeff,double *prime,int m,
  178. double *data,long n){
  179. /* in: coeff[0...m-1] LPC coefficients
  180. prime[0...m-1] initial values (allocated size of n+m-1)
  181. data[0...n-1] residuals from LPC prediction
  182. out: data[0...n-1] data samples */
  183. long i,j,o,p;
  184. double y;
  185. double *work=alloca(sizeof(double)*(m+n));
  186. if(!prime)
  187. for(i=0;i<m;i++)
  188. work[i]=0.;
  189. else
  190. for(i=0;i<m;i++)
  191. work[i]=prime[i];
  192. for(i=0;i<n;i++){
  193. y=data[i];
  194. o=i;
  195. p=m;
  196. for(j=0;j<m;j++)
  197. y-=work[o++]*coeff[--p];
  198. data[i]=work[o]=y;
  199. }
  200. }