Version 0.1.12 from FTP
[asterisk/asterisk.git] / codecs / lpc10 / pitsyn.c
1 /*
2
3 $Log$
4 Revision 1.11  2000/01/05 08:20:39  markster
5 Version 0.1.12 from FTP
6
7 Revision 1.2  2000/01/05 08:20:39  markster
8 Some OSS fixes and a few lpc changes to make it actually work
9
10  * Revision 1.2  1996/08/20  20:40:12  jaf
11  * Removed all static local variables that were SAVE'd in the Fortran
12  * code, and put them in struct lpc10_decoder_state that is passed as an
13  * argument.
14  *
15  * Removed init function, since all initialization is now done in
16  * init_lpc10_decoder_state().
17  *
18  * Revision 1.1  1996/08/19  22:31:12  jaf
19  * Initial revision
20  *
21
22 */
23
24 #ifdef P_R_O_T_O_T_Y_P_E_S
25 extern int pitsyn_(integer *order, integer *voice, integer *pitch, real *rms, real *rc, integer *lframe, integer *ivuv, integer *ipiti, real *rmsi, real *rci, integer *nout, real *ratio, struct lpc10_decoder_state *st);
26 #endif
27
28 /*  -- translated by f2c (version 19951025).
29    You must link the resulting object file with the libraries:
30         -lf2c -lm   (in that order)
31 */
32
33 #include "f2c.h"
34
35 /* ***************************************************************** */
36
37 /*      PITSYN Version 53 */
38
39 /* $Log$
40  * Revision 1.11  2000/01/05 08:20:39  markster
41  * Version 0.1.12 from FTP
42  *
43 /* Revision 1.2  2000/01/05 08:20:39  markster
44 /* Some OSS fixes and a few lpc changes to make it actually work
45 /*
46  * Revision 1.2  1996/08/20  20:40:12  jaf
47  * Removed all static local variables that were SAVE'd in the Fortran
48  * code, and put them in struct lpc10_decoder_state that is passed as an
49  * argument.
50  *
51  * Removed init function, since all initialization is now done in
52  * init_lpc10_decoder_state().
53  *
54  * Revision 1.1  1996/08/19  22:31:12  jaf
55  * Initial revision
56  * */
57 /* Revision 1.2  1996/03/25  18:49:07  jaf */
58 /* Added commments about which indices of array arguments are read or */
59 /* written. */
60
61 /* Rearranged local variable declarations to indicate which need to be */
62 /* saved from one invocation to the next.  Added entry INITPITSYN to */
63 /* reinitialize local state variables, if desired. */
64
65 /* Added lots of comments about proving that the maximum number of pitch */
66 /* periods (NOUT) that can be returned is 16.  The call to STOP that */
67 /* could happen if NOUT got too large was removed as a result. */
68
69 /* Also proved that the total number of samples returned from N calls, */
70 /* each with identical values of LFRAME, will always be in the range */
71 /* N*LFRAME-MAXPIT+1 to N*LFRAME. */
72
73 /* Revision 1.1  1996/02/07 14:48:18  jaf */
74 /* Initial revision */
75
76
77 /* ***************************************************************** */
78
79 /*   Synthesize a single pitch epoch */
80
81 /* Input: */
82 /*  ORDER  - Synthesis order (number of RC's) */
83 /*  VOICE  - Half frame voicing decisions */
84 /*           Indices 1 through 2 read. */
85 /*  LFRAME - Length of speech buffer */
86 /* Input/Output: */
87 /*  PITCH  - Pitch */
88 /*           This value should be in the range MINPIT (20) to MAXPIT */
89 /*           (156), inclusive. */
90 /*           PITCH can be modified under some conditions. */
91 /*  RMS    - Energy  (can be modified) */
92 /*           RMS is changed to 1 if the value passed in is less than 1. */
93 /*  RC     - Reflection coefficients */
94 /*           Indices 1 through ORDER can be temporarily overwritten with */
95 /*           RCO, and then replaced with original values, under some */
96 /*           conditions. */
97 /* Output: */
98 /*  IVUV   - Pitch epoch voicing decisions */
99 /*           Indices (I) of IVUV, IPITI, and RMSI are written, */
100 /*           and indices (J,I) of RCI are written, */
101 /*           where I ranges from 1 to NOUT, and J ranges from 1 to ORDER. */
102 /*  IPITI  - Pitch epoch length */
103 /*  RMSI   - Pitch epoch energy */
104 /*  RCI    - Pitch epoch RC's */
105 /*  NOUT   - Number of pitch periods in this frame */
106 /*           This is at least 0, at least 1 if MAXPIT .LT. LFRAME (this */
107 /*           is currently true on every call), and can never be more than */
108 /*           (LFRAME+MAXPIT-1)/PITCH, which is currently 16 with */
109 /*           LFRAME=180, MAXPIT=156, and PITCH .GE. 20, as SYNTHS */
110 /*           guarantees when it calls this subroutine. */
111 /*  RATIO  - Previous to present energy ratio */
112 /*           Always assigned a value. */
113
114 /* Subroutine */ int pitsyn_(integer *order, integer *voice, 
115         integer *pitch, real *rms, real *rc, integer *lframe, integer *ivuv, 
116         integer *ipiti, real *rmsi, real *rci, integer *nout, real *ratio,
117                                struct lpc10_decoder_state *st)
118 {
119     /* Initialized data */
120
121     real *rmso;
122     logical *first;
123
124     /* System generated locals */
125     integer rci_dim1, rci_offset, i__1, i__2;
126     real r__1;
127
128     /* Builtin functions */
129     double log(doublereal), exp(doublereal);
130
131     /* Local variables */
132     real alrn, alro, yarc[10], prop;
133     integer i__, j, vflag, jused, lsamp;
134     integer *jsamp;
135     real slope;
136     integer *ipito;
137     real uvpit;
138     integer ip, nl, ivoice;
139     integer *ivoico;
140     integer istart;
141     real *rco;
142     real xxy;
143
144 /*       Arguments */
145 /* $Log$
146  * Revision 1.11  2000/01/05 08:20:39  markster
147  * Version 0.1.12 from FTP
148  *
149 /* Revision 1.2  2000/01/05 08:20:39  markster
150 /* Some OSS fixes and a few lpc changes to make it actually work
151 /*
152  * Revision 1.2  1996/08/20  20:40:12  jaf
153  * Removed all static local variables that were SAVE'd in the Fortran
154  * code, and put them in struct lpc10_decoder_state that is passed as an
155  * argument.
156  *
157  * Removed init function, since all initialization is now done in
158  * init_lpc10_decoder_state().
159  *
160  * Revision 1.1  1996/08/19  22:31:12  jaf
161  * Initial revision
162  * */
163 /* Revision 1.3  1996/03/29  22:03:47  jaf */
164 /* Removed definitions for any constants that were no longer used. */
165
166 /* Revision 1.2  1996/03/26  19:34:33  jaf */
167 /* Added comments indicating which constants are not needed in an */
168 /* application that uses the LPC-10 coder. */
169
170 /* Revision 1.1  1996/02/07  14:43:51  jaf */
171 /* Initial revision */
172
173 /*   LPC Configuration parameters: */
174 /* Frame size, Prediction order, Pitch period */
175 /*       Local variables that need not be saved */
176 /*       LSAMP is initialized in the IF (FIRST) THEN clause, but it is */
177 /*       not used the first time through, and it is given a value before 
178 */
179 /*       use whenever FIRST is .FALSE., so it appears unnecessary to */
180 /*       assign it a value when FIRST is .TRUE. */
181 /*       Local state */
182 /* FIRST  - .TRUE. only on first call to PITSYN. */
183 /* IVOICO - Previous VOICE(2) value. */
184 /* IPITO  - Previous PITCH value. */
185 /* RMSO   - Previous RMS value. */
186 /* RCO    - Previous RC values. */
187
188 /* JSAMP  - If this routine is called N times with identical values of */
189 /*          LFRAME, then the total length of all pitch periods returned */
190 /*          is always N*LFRAME-JSAMP, and JSAMP is always in the range 0 
191 */
192 /*          to MAXPIT-1 (see below for why this is so).  Thus JSAMP is */
193 /*          the number of samples "left over" from the previous call to */
194 /*          PITSYN, that haven't been "used" in a pitch period returned */
195 /*          from this subroutine.  Every time this subroutine is called, 
196 */
197 /*          it returns pitch periods with a total length of at most */
198 /*          LFRAME+JSAMP. */
199
200 /* IVOICO, IPITO, RCO, and JSAMP need not be assigned an initial value */
201 /* with a DATA statement, because they are always initialized on the */
202 /* first call to PITSYN. */
203
204 /* FIRST and RMSO should be initialized with DATA statements, because */
205 /* even on the first call, they are used before being initialized. */
206     /* Parameter adjustments */
207     if (rc) {
208         --rc;
209         }
210     if (rci) {
211         rci_dim1 = *order;
212         rci_offset = rci_dim1 + 1;
213         rci -= rci_offset;
214         }
215     if (voice) {
216         --voice;
217         }
218     if (ivuv) {
219         --ivuv;
220         }
221     if (ipiti) {
222         --ipiti;
223         }
224     if (rmsi) {
225         --rmsi;
226         }
227
228     /* Function Body */
229     ivoico = &(st->ivoico);
230     ipito = &(st->ipito);
231     rmso = &(st->rmso);
232     rco = &(st->rco[0]);
233     jsamp = &(st->jsamp);
234     first = &(st->first_pitsyn);
235
236     if (*rms < 1.f) {
237         *rms = 1.f;
238     }
239     if (*rmso < 1.f) {
240         *rmso = 1.f;
241     }
242     uvpit = 0.f;
243     *ratio = *rms / (*rmso + 8.f);
244     if (*first) {
245         lsamp = 0;
246         ivoice = voice[2];
247         if (ivoice == 0) {
248             *pitch = *lframe / 4;
249         }
250         *nout = *lframe / *pitch;
251         *jsamp = *lframe - *nout * *pitch;
252
253 /*          SYNTHS only calls this subroutine with PITCH in the range 
254 20 */
255 /*          to 156.  LFRAME = MAXFRM = 180, so NOUT is somewhere in th
256 e */
257 /*          range 1 to 9. */
258
259 /*          JSAMP is "LFRAME mod PITCH", so it is in the range 0 to */
260 /*          (PITCH-1), or 0 to MAXPIT-1=155, after the first call. */
261
262         i__1 = *nout;
263         for (i__ = 1; i__ <= i__1; ++i__) {
264             i__2 = *order;
265             for (j = 1; j <= i__2; ++j) {
266                 rci[j + i__ * rci_dim1] = rc[j];
267             }
268             ivuv[i__] = ivoice;
269             ipiti[i__] = *pitch;
270             rmsi[i__] = *rms;
271         }
272         *first = FALSE_;
273     } else {
274         vflag = 0;
275         lsamp = *lframe + *jsamp;
276         slope = (*pitch - *ipito) / (real) lsamp;
277         *nout = 0;
278         jused = 0;
279         istart = 1;
280         if (voice[1] == *ivoico && voice[2] == voice[1]) {
281             if (voice[2] == 0) {
282 /* SSUV - -   0  ,  0  ,  0 */
283                 *pitch = *lframe / 4;
284                 *ipito = *pitch;
285                 if (*ratio > 8.f) {
286                     *rmso = *rms;
287                 }
288             }
289 /* SSVC - -   1  ,  1  ,  1 */
290             slope = (*pitch - *ipito) / (real) lsamp;
291             ivoice = voice[2];
292         } else {
293             if (*ivoico != 1) {
294                 if (*ivoico == voice[1]) {
295 /* UV2VC2 - -  0  ,  0  ,  1 */
296                     nl = lsamp - *lframe / 4;
297                 } else {
298 /* UV2VC1 - -  0  ,  1  ,  1 */
299                     nl = lsamp - *lframe * 3 / 4;
300                 }
301                 ipiti[1] = nl / 2;
302                 ipiti[2] = nl - ipiti[1];
303                 ivuv[1] = 0;
304                 ivuv[2] = 0;
305                 rmsi[1] = *rmso;
306                 rmsi[2] = *rmso;
307                 i__1 = *order;
308                 for (i__ = 1; i__ <= i__1; ++i__) {
309                     rci[i__ + rci_dim1] = rco[i__ - 1];
310                     rci[i__ + (rci_dim1 << 1)] = rco[i__ - 1];
311                     rco[i__ - 1] = rc[i__];
312                 }
313                 slope = 0.f;
314                 *nout = 2;
315                 *ipito = *pitch;
316                 jused = nl;
317                 istart = nl + 1;
318                 ivoice = 1;
319             } else {
320                 if (*ivoico != voice[1]) {
321 /* VC2UV1 - -   1  ,  0  ,  0 */
322                     lsamp = *lframe / 4 + *jsamp;
323                 } else {
324 /* VC2UV2 - -   1  ,  1  ,  0 */
325                     lsamp = *lframe * 3 / 4 + *jsamp;
326                 }
327                 i__1 = *order;
328                 for (i__ = 1; i__ <= i__1; ++i__) {
329                     yarc[i__ - 1] = rc[i__];
330                     rc[i__] = rco[i__ - 1];
331                 }
332                 ivoice = 1;
333                 slope = 0.f;
334                 vflag = 1;
335             }
336         }
337 /* Here is the value of most variables that are used below, depending 
338 on */
339 /* the values of IVOICO, VOICE(1), and VOICE(2).  VOICE(1) and VOICE(2
340 ) */
341 /* are input arguments, and IVOICO is the value of VOICE(2) on the */
342 /* previous call (see notes for the IF (NOUT .NE. 0) statement near th
343 e */
344 /* end).  Each of these three values is either 0 or 1.  These three */
345 /* values below are given as 3-bit long strings, in the order IVOICO, 
346 */
347 /* VOICE(1), and VOICE(2).  It appears that the code above assumes tha
348 t */
349 /* the bit sequences 010 and 101 never occur, but I wonder whether a 
350 */
351 /* large enough number of bit errors in the channel could cause such a
352  */
353 /* thing to happen, and if so, could that cause NOUT to ever go over 1
354 1? */
355
356 /* Note that all of the 180 values in the table are really LFRAME, but
357  */
358 /* 180 has fewer characters, and it makes the table a little more */
359 /* concrete.  If LFRAME is ever changed, keep this in mind.  Similarly
360 , */
361 /* 135's are 3*LFRAME/4, and 45's are LFRAME/4.  If LFRAME is not a */
362 /* multiple of 4, then the 135 for NL-JSAMP is actually LFRAME-LFRAME/
363 4, */
364 /* and the 45 for NL-JSAMP is actually LFRAME-3*LFRAME/4. */
365
366 /* Note that LSAMP-JSAMP is given as the variable.  This was just for 
367 */
368 /* brevity, to avoid adding "+JSAMP" to all of the column entries. */
369 /* Similarly for NL-JSAMP. */
370
371 /* Variable    | 000  001    011,010  111       110       100,101 */
372 /* ------------+-------------------------------------------------- */
373 /* ISTART      | 1    NL+1   NL+1     1         1         1 */
374 /* LSAMP-JSAMP | 180  180    180      180       135       45 */
375 /* IPITO       | 45   PITCH  PITCH    oldPITCH  oldPITCH  oldPITCH */
376 /* SLOPE       | 0    0      0        seebelow  0         0 */
377 /* JUSED       | 0    NL     NL       0         0         0 */
378 /* PITCH       | 45   PITCH  PITCH    PITCH     PITCH     PITCH */
379 /* NL-JSAMP    | --   135    45       --        --        -- */
380 /* VFLAG       | 0    0      0        0         1         1 */
381 /* NOUT        | 0    2      2        0         0         0 */
382 /* IVOICE      | 0    1      1        1         1         1 */
383
384 /* while_loop  | once once   once     once      twice     twice */
385
386 /* ISTART      | --   --     --       --        JUSED+1   JUSED+1 */
387 /* LSAMP-JSAMP | --   --     --       --        180       180 */
388 /* IPITO       | --   --     --       --        oldPITCH  oldPITCH */
389 /* SLOPE       | --   --     --       --        0         0 */
390 /* JUSED       | --   --     --       --        ??        ?? */
391 /* PITCH       | --   --     --       --        PITCH     PITCH */
392 /* NL-JSAMP    | --   --     --       --        --        -- */
393 /* VFLAG       | --   --     --       --        0         0 */
394 /* NOUT        | --   --     --       --        ??        ?? */
395 /* IVOICE      | --   --     --       --        0         0 */
396
397
398 /* UVPIT is always 0.0 on the first pass through the DO WHILE (.TRUE.)
399  */
400 /* loop below. */
401
402 /* The only possible non-0 value of SLOPE (in column 111) is */
403 /* (PITCH-IPITO)/FLOAT(LSAMP) */
404
405 /* Column 101 is identical to 100.  Any good properties we can prove 
406 */
407 /* for 100 will also hold for 101.  Similarly for 010 and 011. */
408
409 /* SYNTHS calls this subroutine with PITCH restricted to the range 20 
410 to */
411 /* 156.  IPITO is similarly restricted to this range, after the first 
412 */
413 /* call.  IP below is also restricted to this range, given the */
414 /* definitions of IPITO, SLOPE, UVPIT, and that I is in the range ISTA
415 RT */
416 /* to LSAMP. */
417
418         while(TRUE_) {
419
420 /*             JUSED is the total length of all pitch periods curr
421 ently */
422 /*             in the output arrays, in samples. */
423
424 /*             An invariant of the DO I = ISTART,LSAMP loop below,
425  under */
426 /*             the condition that IP is always in the range 1 thro
427 ugh */
428 /*             MAXPIT, is: */
429
430 /*             (I - MAXPIT) .LE. JUSED .LE. (I-1) */
431
432 /*             Note that the final value of I is LSAMP+1, so that 
433 after */
434 /*             the DO loop is complete, we know: */
435
436 /*             (LSAMP - MAXPIT + 1) .LE. JUSED .LE. LSAMP */
437
438             i__1 = lsamp;
439             for (i__ = istart; i__ <= i__1; ++i__) {
440                 r__1 = *ipito + slope * i__;
441                 ip = r__1 + .5f;
442                 if (uvpit != 0.f) {
443                     ip = uvpit;
444                 }
445                 if (ip <= i__ - jused) {
446                     ++(*nout);
447
448 /*                   The following check is no longer nece
449 ssary, now that */
450 /*                   we can prove that NOUT will never go 
451 over 16. */
452
453 /*                  IF (NOUT .GT. 16) STOP 'PITSYN: too many epochs' 
454 */
455
456                     ipiti[*nout] = ip;
457                     *pitch = ip;
458                     ivuv[*nout] = ivoice;
459                     jused += ip;
460                     prop = (jused - ip / 2) / (real) lsamp;
461                     i__2 = *order;
462                     for (j = 1; j <= i__2; ++j) {
463                         alro = log((rco[j - 1] + 1) / (1 - rco[j - 1]));
464                         alrn = log((rc[j] + 1) / (1 - rc[j]));
465                         xxy = alro + prop * (alrn - alro);
466                         xxy = exp(xxy);
467                         rci[j + *nout * rci_dim1] = (xxy - 1) / (xxy + 1);
468                     }
469                     rmsi[*nout] = log(*rmso) + prop * (log(*rms) - log(*rmso));
470                     rmsi[*nout] = exp(rmsi[*nout]);
471                 }
472             }
473             if (vflag != 1) {
474                 goto L100;
475             }
476
477 /*             I want to prove what range UVPIT must lie in after 
478 the */
479 /*             assignments to it below.  To do this, I must determ
480 ine */
481 /*             what range (LSAMP-ISTART) must lie in, after the */
482 /*             assignments to ISTART and LSAMP below. */
483
484 /*             Let oldLSAMP be the value of LSAMP at this point in
485  the */
486 /*             execution.  This is 135+JSAMP in state 110, or 45+J
487 SAMP in */
488 /*             states 100 or 101. */
489
490 /*             Given the loop invariant on JUSED above, we know th
491 at: */
492
493 /*             (oldLSAMP - MAXPIT + 1) .LE. JUSED .LE. oldLSAMP */
494
495 /*             ISTART is one more than this. */
496
497 /*             Let newLSAMP be the value assigned to LSAMP below. 
498  This */
499 /*             is 180+JSAMP.  Thus (newLSAMP-oldLSAMP) is either 4
500 5 or */
501 /*             135, depending on the state. */
502
503 /*             Thus, the range of newLSAMP-ISTART is: */
504
505 /*             (newLSAMP-(oldLSAMP+1)) .LE. newLSAMP-ISTART */
506 /*             .LE. (newLSAMP-(oldLSAMP - MAXPIT + 2)) */
507
508 /*             or: */
509
510 /*             46 .LE. newLSAMP-ISTART .LE. 133+MAXPIT .EQ. 289 */
511
512 /*             Therefore, UVPIT is in the range 23 to 144 after th
513 e first */
514 /*             assignment to UVPIT below, and after the conditiona
515 l */
516 /*             assignment, it is in the range 23 to 90. */
517
518 /*             The important thing is that it is in the range 20 t
519 o 156, */
520 /*             so that in the loop above, IP is always in this ran
521 ge. */
522
523             vflag = 0;
524             istart = jused + 1;
525             lsamp = *lframe + *jsamp;
526             slope = 0.f;
527             ivoice = 0;
528             uvpit = (real) ((lsamp - istart) / 2);
529             if (uvpit > 90.f) {
530                 uvpit /= 2;
531             }
532             *rmso = *rms;
533             i__1 = *order;
534             for (i__ = 1; i__ <= i__1; ++i__) {
535                 rc[i__] = yarc[i__ - 1];
536                 rco[i__ - 1] = yarc[i__ - 1];
537             }
538         }
539 L100:
540         *jsamp = lsamp - jused;
541     }
542 /*       Given that the maximum pitch period MAXPIT .LT. LFRAME (this is 
543 */
544 /*       currently true on every call, since SYNTHS always sets */
545 /*       LFRAME=180), NOUT will always be .GE. 1 at this point. */
546     if (*nout != 0) {
547         *ivoico = voice[2];
548         *ipito = *pitch;
549         *rmso = *rms;
550         i__1 = *order;
551         for (i__ = 1; i__ <= i__1; ++i__) {
552             rco[i__ - 1] = rc[i__];
553         }
554     }
555     return 0;
556 } /* pitsyn_ */