lpc.c 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342
  1. /*
  2. * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
  3. * Universitaet Berlin. See the accompanying file "COPYRIGHT" for
  4. * details. THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
  5. */
  6. /* $Header: /tmp_amd/presto/export/kbs/jutta/src/gsm/RCS/lpc.c,v 1.5 1994/12/30 23:14:54 jutta Exp $ */
  7. #include "config.h"
  8. #include <stdio.h>
  9. #include <assert.h>
  10. #include "private.h"
  11. #include "gsm.h"
  12. #include "proto.h"
  13. #undef P
  14. /*
  15. * 4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
  16. */
  17. /* 4.2.4 */
  18. static void Autocorrelation P2((s, L_ACF),
  19. word * s, /* [0..159] IN/OUT */
  20. longword * L_ACF) /* [0..8] OUT */
  21. /*
  22. * The goal is to compute the array L_ACF[k]. The signal s[i] must
  23. * be scaled in order to avoid an overflow situation.
  24. */
  25. {
  26. register int k, i;
  27. word temp, smax, scalauto;
  28. #ifdef USE_FLOAT_MUL
  29. float float_s[160];
  30. #endif
  31. /* Dynamic scaling of the array s[0..159]
  32. */
  33. /* Search for the maximum.
  34. */
  35. smax = 0;
  36. for (k = 0; k <= 159; k++) {
  37. temp = GSM_ABS( s[k] );
  38. if (temp > smax) smax = temp;
  39. }
  40. /* Computation of the scaling factor.
  41. */
  42. if (smax == 0) scalauto = 0;
  43. else {
  44. assert(smax > 0);
  45. scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
  46. }
  47. /* Scaling of the array s[0...159]
  48. */
  49. if (scalauto > 0) {
  50. # ifdef USE_FLOAT_MUL
  51. # define SCALE(n) \
  52. case n: for (k = 0; k <= 159; k++) \
  53. float_s[k] = (float) \
  54. (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
  55. break;
  56. # else
  57. # define SCALE(n) \
  58. case n: for (k = 0; k <= 159; k++) \
  59. s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
  60. break;
  61. # endif /* USE_FLOAT_MUL */
  62. switch (scalauto) {
  63. SCALE(1)
  64. SCALE(2)
  65. SCALE(3)
  66. SCALE(4)
  67. }
  68. # undef SCALE
  69. }
  70. # ifdef USE_FLOAT_MUL
  71. else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
  72. # endif
  73. /* Compute the L_ACF[..].
  74. */
  75. {
  76. # ifdef USE_FLOAT_MUL
  77. register float * sp = float_s;
  78. register float sl = *sp;
  79. # define STEP(k) L_ACF[k] += (longword)(sl * sp[ -(k) ]);
  80. # else
  81. word * sp = s;
  82. word sl = *sp;
  83. # define STEP(k) L_ACF[k] += ((longword)sl * sp[ -(k) ]);
  84. # endif
  85. # define NEXTI sl = *++sp
  86. for (k = 9; k--; L_ACF[k] = 0) ;
  87. STEP (0);
  88. NEXTI;
  89. STEP(0); STEP(1);
  90. NEXTI;
  91. STEP(0); STEP(1); STEP(2);
  92. NEXTI;
  93. STEP(0); STEP(1); STEP(2); STEP(3);
  94. NEXTI;
  95. STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
  96. NEXTI;
  97. STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
  98. NEXTI;
  99. STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
  100. NEXTI;
  101. STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
  102. for (i = 8; i <= 159; i++) {
  103. NEXTI;
  104. STEP(0);
  105. STEP(1); STEP(2); STEP(3); STEP(4);
  106. STEP(5); STEP(6); STEP(7); STEP(8);
  107. }
  108. for (k = 9; k--; L_ACF[k] <<= 1) ;
  109. }
  110. /* Rescaling of the array s[0..159]
  111. */
  112. if (scalauto > 0) {
  113. assert(scalauto <= 4);
  114. for (k = 160; k--; *s++ <<= scalauto) ;
  115. }
  116. }
  117. #if defined(USE_FLOAT_MUL) && defined(FAST)
  118. static void Fast_Autocorrelation P2((s, L_ACF),
  119. word * s, /* [0..159] IN/OUT */
  120. longword * L_ACF) /* [0..8] OUT */
  121. {
  122. register int k, i;
  123. float f_L_ACF[9];
  124. float scale;
  125. float s_f[160];
  126. register float *sf = s_f;
  127. for (i = 0; i < 160; ++i) sf[i] = s[i];
  128. for (k = 0; k <= 8; k++) {
  129. register float L_temp2 = 0;
  130. register float *sfl = sf - k;
  131. for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
  132. f_L_ACF[k] = L_temp2;
  133. }
  134. scale = MAX_LONGWORD / f_L_ACF[0];
  135. for (k = 0; k <= 8; k++) {
  136. L_ACF[k] = f_L_ACF[k] * scale;
  137. }
  138. }
  139. #endif /* defined (USE_FLOAT_MUL) && defined (FAST) */
  140. /* 4.2.5 */
  141. static void Reflection_coefficients P2( (L_ACF, r),
  142. longword * L_ACF, /* 0...8 IN */
  143. register word * r /* 0...7 OUT */
  144. )
  145. {
  146. register int i, m, n;
  147. register word temp;
  148. register longword ltmp;
  149. word ACF[9]; /* 0..8 */
  150. word P[ 9]; /* 0..8 */
  151. word K[ 9]; /* 2..8 */
  152. /* Schur recursion with 16 bits arithmetic.
  153. */
  154. if (L_ACF[0] == 0) {
  155. for (i = 8; i--; *r++ = 0) ;
  156. return;
  157. }
  158. assert( L_ACF[0] != 0 );
  159. temp = gsm_norm( L_ACF[0] );
  160. assert(temp >= 0 && temp < 32);
  161. /* ? overflow ? */
  162. for (i = 0; i <= 8; i++) ACF[i] = SASR( L_ACF[i] << temp, 16 );
  163. /* Initialize array P[..] and K[..] for the recursion.
  164. */
  165. for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
  166. for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
  167. /* Compute reflection coefficients
  168. */
  169. for (n = 1; n <= 8; n++, r++) {
  170. temp = P[1];
  171. temp = GSM_ABS(temp);
  172. if (P[0] < temp) {
  173. for (i = n; i <= 8; i++) *r++ = 0;
  174. return;
  175. }
  176. *r = gsm_div( temp, P[0] );
  177. assert(*r >= 0);
  178. if (P[1] > 0) *r = -*r; /* r[n] = sub(0, r[n]) */
  179. assert (*r != MIN_WORD);
  180. if (n == 8) return;
  181. /* Schur recursion
  182. */
  183. temp = GSM_MULT_R( P[1], *r );
  184. P[0] = GSM_ADD( P[0], temp );
  185. for (m = 1; m <= 8 - n; m++) {
  186. temp = GSM_MULT_R( K[ m ], *r );
  187. P[m] = GSM_ADD( P[ m+1 ], temp );
  188. temp = GSM_MULT_R( P[ m+1 ], *r );
  189. K[m] = GSM_ADD( K[ m ], temp );
  190. }
  191. }
  192. }
  193. /* 4.2.6 */
  194. static void Transformation_to_Log_Area_Ratios P1((r),
  195. register word * r /* 0..7 IN/OUT */
  196. )
  197. /*
  198. * The following scaling for r[..] and LAR[..] has been used:
  199. *
  200. * r[..] = integer( real_r[..]*32768. ); -1 <= real_r < 1.
  201. * LAR[..] = integer( real_LAR[..] * 16384 );
  202. * with -1.625 <= real_LAR <= 1.625
  203. */
  204. {
  205. register word temp;
  206. register int i;
  207. /* Computation of the LAR[0..7] from the r[0..7]
  208. */
  209. for (i = 1; i <= 8; i++, r++) {
  210. temp = *r;
  211. temp = GSM_ABS(temp);
  212. assert(temp >= 0);
  213. if (temp < 22118) {
  214. temp >>= 1;
  215. } else if (temp < 31130) {
  216. assert( temp >= 11059 );
  217. temp -= 11059;
  218. } else {
  219. assert( temp >= 26112 );
  220. temp -= 26112;
  221. temp <<= 2;
  222. }
  223. *r = *r < 0 ? -temp : temp;
  224. assert( *r != MIN_WORD );
  225. }
  226. }
  227. /* 4.2.7 */
  228. static void Quantization_and_coding P1((LAR),
  229. register word * LAR /* [0..7] IN/OUT */
  230. )
  231. {
  232. register word temp;
  233. longword ltmp;
  234. /* This procedure needs four tables; the following equations
  235. * give the optimum scaling for the constants:
  236. *
  237. * A[0..7] = integer( real_A[0..7] * 1024 )
  238. * B[0..7] = integer( real_B[0..7] * 512 )
  239. * MAC[0..7] = maximum of the LARc[0..7]
  240. * MIC[0..7] = minimum of the LARc[0..7]
  241. */
  242. # undef STEP
  243. # define STEP( A, B, MAC, MIC ) \
  244. temp = GSM_MULT( A, *LAR ); \
  245. temp = GSM_ADD( temp, B ); \
  246. temp = GSM_ADD( temp, 256 ); \
  247. temp = SASR( temp, 9 ); \
  248. *LAR = temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
  249. LAR++;
  250. STEP( 20480, 0, 31, -32 );
  251. STEP( 20480, 0, 31, -32 );
  252. STEP( 20480, 2048, 15, -16 );
  253. STEP( 20480, -2560, 15, -16 );
  254. STEP( 13964, 94, 7, -8 );
  255. STEP( 15360, -1792, 7, -8 );
  256. STEP( 8534, -341, 3, -4 );
  257. STEP( 9036, -1144, 3, -4 );
  258. # undef STEP
  259. }
  260. void Gsm_LPC_Analysis P3((S, s,LARc),
  261. struct gsm_state *S,
  262. word * s, /* 0..159 signals IN/OUT */
  263. word * LARc) /* 0..7 LARc's OUT */
  264. {
  265. longword L_ACF[9];
  266. #if defined(USE_FLOAT_MUL) && defined(FAST)
  267. if (S->fast) Fast_Autocorrelation (s, L_ACF );
  268. else
  269. #endif
  270. Autocorrelation (s, L_ACF );
  271. Reflection_coefficients (L_ACF, LARc );
  272. Transformation_to_Log_Area_Ratios (LARc);
  273. Quantization_and_coding (LARc);
  274. }