Quantizing high-band excitation gains (SB-CELP).
[speexdsp.git] / libspeex / sb_celp.c
1 /* Copyright (C) 2002 Jean-Marc Valin 
2    File: speex.c
3
4    This library is free software; you can redistribute it and/or
5    modify it under the terms of the GNU Lesser General Public
6    License as published by the Free Software Foundation; either
7    version 2.1 of the License, or (at your option) any later version.
8    
9    This library is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12    Lesser General Public License for more details.
13    
14    You should have received a copy of the GNU Lesser General Public
15    License along with this library; if not, write to the Free Software
16    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17 */
18
19
20 #include "speex.h"
21 #include "sb_celp.h"
22 #include "stdlib.h"
23 #include "filters.h"
24 #include <math.h>
25 #include "lpc.h"
26 #include "lsp.h"
27 #include <stdio.h>
28 #include "stack_alloc.h"
29 #include "cb_search.h"
30 #include "quant_lsp.h"
31
32 #ifndef M_PI
33 #define M_PI           3.14159265358979323846  /* pi */
34 #endif
35
36 extern float stoc[];
37
38 #define sqr(x) ((x)*(x))
39
40 #if 0
41 #define QMF_ORDER 32
42 static float h0[32] = {
43    0.0006910579, -0.001403793,
44    -0.001268303, 0.004234195,
45    0.001414246, -0.009458318,
46    -0.0001303859, 0.01798145,
47    -0.004187483, -0.03123862,
48    0.01456844, 0.05294745,
49    -0.03934878, -0.09980243,
50    0.1285579, 0.4664053,
51    0.4664053, 0.1285579,
52    -0.09980243, -0.03934878,
53    0.05294745, 0.01456844,
54    -0.03123862, -0.004187483,
55    0.01798145, -0.0001303859,
56    -0.009458318, 0.001414246,
57    0.004234195, -0.001268303,
58    -0.001403793, 0.0006910579
59 };
60
61 static float h1[32] = {
62    0.0006910579, 0.001403793,
63    -0.001268303, -0.004234195,
64    0.001414246, 0.009458318,
65    -0.0001303859, -0.01798145,
66    -0.004187483, 0.03123862,
67    0.01456844, -0.05294745,
68    -0.03934878, 0.09980243,
69    0.1285579, -0.4664053,
70    0.4664053, -0.1285579,
71    -0.09980243, 0.03934878,
72    0.05294745, -0.01456844,
73    -0.03123862, 0.004187483,
74    0.01798145, 0.0001303859,
75    -0.009458318, -0.001414246,
76    0.004234195, 0.001268303,
77    -0.001403793, -0.0006910579
78 };
79 #else 
80 #define QMF_ORDER 64
81 static float h0[64] = {
82    3.596189e-05, -0.0001123515,
83    -0.0001104587, 0.0002790277,
84    0.0002298438, -0.0005953563,
85    -0.0003823631, 0.00113826,
86    0.0005308539, -0.001986177,
87    -0.0006243724, 0.003235877,
88    0.0005743159, -0.004989147,
89    -0.0002584767, 0.007367171,
90    -0.0004857935, -0.01050689,
91    0.001894714, 0.01459396,
92    -0.004313674, -0.01994365,
93    0.00828756, 0.02716055,
94    -0.01485397, -0.03764973,
95    0.026447, 0.05543245,
96    -0.05095487, -0.09779096,
97    0.1382363, 0.4600981,
98    0.4600981, 0.1382363,
99    -0.09779096, -0.05095487,
100    0.05543245, 0.026447,
101    -0.03764973, -0.01485397,
102    0.02716055, 0.00828756,
103    -0.01994365, -0.004313674,
104    0.01459396, 0.001894714,
105    -0.01050689, -0.0004857935,
106    0.007367171, -0.0002584767,
107    -0.004989147, 0.0005743159,
108    0.003235877, -0.0006243724,
109    -0.001986177, 0.0005308539,
110    0.00113826, -0.0003823631,
111    -0.0005953563, 0.0002298438,
112    0.0002790277, -0.0001104587,
113    -0.0001123515, 3.596189e-05
114 };
115
116 static float h1[64] = {
117    3.596189e-05, 0.0001123515,
118    -0.0001104587, -0.0002790277,
119    0.0002298438, 0.0005953563,
120    -0.0003823631, -0.00113826,
121    0.0005308539, 0.001986177,
122    -0.0006243724, -0.003235877,
123    0.0005743159, 0.004989147,
124    -0.0002584767, -0.007367171,
125    -0.0004857935, 0.01050689,
126    0.001894714, -0.01459396,
127    -0.004313674, 0.01994365,
128    0.00828756, -0.02716055,
129    -0.01485397, 0.03764973,
130    0.026447, -0.05543245,
131    -0.05095487, 0.09779096,
132    0.1382363, -0.4600981,
133    0.4600981, -0.1382363,
134    -0.09779096, 0.05095487,
135    0.05543245, -0.026447,
136    -0.03764973, 0.01485397,
137    0.02716055, -0.00828756,
138    -0.01994365, 0.004313674,
139    0.01459396, -0.001894714,
140    -0.01050689, 0.0004857935,
141    0.007367171, 0.0002584767,
142    -0.004989147, -0.0005743159,
143    0.003235877, 0.0006243724,
144    -0.001986177, -0.0005308539,
145    0.00113826, 0.0003823631,
146    -0.0005953563, -0.0002298438,
147    0.0002790277, 0.0001104587,
148    -0.0001123515, -3.596189e-05
149 };
150 #endif
151
152 void sb_encoder_init(SBEncState *st, SpeexMode *mode)
153 {
154    int i;
155    encoder_init(&st->st_low, mode);
156    st->full_frame_size = 2*st->st_low.frameSize;
157    st->frame_size = st->st_low.frameSize;
158    st->subframeSize = 40;
159    st->nbSubframes = 4;
160    st->windowSize = mode->windowSize;
161    st->lpcSize=8;
162
163    st->lag_factor = .002;
164    st->lpc_floor = 1.0001;
165    st->gamma1=.9;
166    st->gamma2=.6;
167    st->first=1;
168    st->stack = calloc(10000, sizeof(float));
169
170    st->x0=calloc(st->full_frame_size, sizeof(float));
171    st->x1=calloc(st->full_frame_size, sizeof(float));
172    st->x0d=calloc(st->frame_size, sizeof(float));
173    st->x1d=calloc(st->frame_size, sizeof(float));
174    st->high=calloc(st->full_frame_size, sizeof(float));
175    st->y0=calloc(st->full_frame_size, sizeof(float));
176    st->y1=calloc(st->full_frame_size, sizeof(float));
177
178    st->h0_mem=calloc(QMF_ORDER, sizeof(float));
179    st->h1_mem=calloc(QMF_ORDER, sizeof(float));
180    st->g0_mem=calloc(QMF_ORDER, sizeof(float));
181    st->g1_mem=calloc(QMF_ORDER, sizeof(float));
182
183    st->buf=calloc(st->windowSize, sizeof(float));
184    st->excBuf=calloc(2*st->frame_size, sizeof(float));
185    st->exc=st->excBuf+st->frame_size;
186    st->exc_alias=calloc(st->frame_size, sizeof(float));
187
188    st->res=calloc(st->frame_size, sizeof(float));
189    st->sw=calloc(st->frame_size, sizeof(float));
190    st->target=calloc(st->frame_size, sizeof(float));
191    st->window=calloc(st->windowSize, sizeof(float));
192    for (i=0;i<st->windowSize;i++)
193       st->window[i]=.5*(1-cos(2*M_PI*i/st->windowSize));
194
195    st->exc_window=calloc(st->frame_size, sizeof(float));
196    for (i=0;i<st->frame_size;i++)
197       st->exc_window[i]=.5*(1-cos(2*M_PI*i/st->frame_size));
198
199    st->lagWindow = malloc((st->lpcSize+1)*sizeof(float));
200    for (i=0;i<st->lpcSize+1;i++)
201       st->lagWindow[i]=exp(-.5*sqr(2*M_PI*st->lag_factor*i));
202
203    st->rc = malloc(st->lpcSize*sizeof(float));
204    st->autocorr = malloc((st->lpcSize+1)*sizeof(float));
205    st->lpc = malloc((st->lpcSize+1)*sizeof(float));
206    st->bw_lpc1 = malloc((st->lpcSize+1)*sizeof(float));
207    st->bw_lpc2 = malloc((st->lpcSize+1)*sizeof(float));
208    st->lsp = malloc(st->lpcSize*sizeof(float));
209    st->qlsp = malloc(st->lpcSize*sizeof(float));
210    st->old_lsp = malloc(st->lpcSize*sizeof(float));
211    st->old_qlsp = malloc(st->lpcSize*sizeof(float));
212    st->interp_lsp = malloc(st->lpcSize*sizeof(float));
213    st->interp_qlsp = malloc(st->lpcSize*sizeof(float));
214    st->interp_lpc = malloc((st->lpcSize+1)*sizeof(float));
215    st->interp_qlpc = malloc((st->lpcSize+1)*sizeof(float));
216
217    st->mem_sp = calloc(st->lpcSize, sizeof(float));
218    st->mem_sp2 = calloc(st->lpcSize, sizeof(float));
219    st->mem_sw = calloc(st->lpcSize, sizeof(float));
220    st->mem_exc = calloc(st->lpcSize, sizeof(float));
221
222 }
223
224 void sb_encoder_destroy(SBEncState *st)
225 {
226    encoder_destroy(&st->st_low);
227    free(st->x0);
228    free(st->x0d);
229    free(st->x1);
230    free(st->x1d);
231    free(st->high);
232    free(st->y0);
233    free(st->y1);
234    free(st->h0_mem);
235    free(st->h1_mem);
236    free(st->g0_mem);
237    free(st->g1_mem);
238    
239    free(st->buf);
240    free(st->window);
241    free(st->excBuf);
242    free(st->sw);
243    free(st->res);
244    free(st->target);
245    free(st->lagWindow);
246    free(st->rc);
247    free(st->autocorr);
248    free(st->lpc);
249    free(st->bw_lpc1);
250    free(st->bw_lpc2);
251    free(st->lsp);
252    free(st->qlsp);
253    free(st->old_lsp);
254    free(st->old_qlsp);
255    free(st->interp_lsp);
256    free(st->interp_qlsp);
257    free(st->interp_lpc);
258    free(st->interp_qlpc);
259
260    free(st->mem_sp);
261    free(st->mem_sw);
262
263    free(st->stack);
264    
265 }
266
267
268 void sb_encode(SBEncState *st, float *in, FrameBits *bits)
269 {
270    int i, roots, sub;
271    /* Compute the two sub-bands by filtering with h0 and h1*/
272    fir_mem(in, h0, st->x0, st->full_frame_size, QMF_ORDER, st->h0_mem);
273    fir_mem(in, h1, st->x1, st->full_frame_size, QMF_ORDER, st->h1_mem);
274    /* Down-sample x0 and x1 */
275    for (i=0;i<st->frame_size;i++)
276    {
277       st->x0d[i]=st->x0[i<<1];
278       st->x1d[i]=st->x1[i<<1];
279    }
280    /* Encode the narrowband part*/
281    encode(&st->st_low, st->x0d, bits);
282
283    /* High-band buffering / sync with low band */
284    for (i=0;i<st->frame_size;i++)
285    {
286       st->excBuf[i]=st->exc[i];
287       st->high[i]=st->high[st->frame_size+i];
288       st->high[st->frame_size+i]=st->x1d[i];
289    }
290    
291    /* Start encoding the high-band */
292
293    for (i=0;i<st->windowSize;i++)
294       st->buf[i] = st->high[i] * st->window[i];
295
296    /* Compute auto-correlation */
297    autocorr(st->buf, st->autocorr, st->lpcSize+1, st->windowSize);
298
299    st->autocorr[0] += 1;        /* prevents NANs */
300    st->autocorr[0] *= st->lpc_floor; /* Noise floor in auto-correlation domain */
301    /* Lag windowing: equivalent to filtering in the power-spectrum domain */
302    for (i=0;i<st->lpcSize+1;i++)
303       st->autocorr[i] *= st->lagWindow[i];
304
305    /* Levinson-Durbin */
306    wld(st->lpc+1, st->autocorr, st->rc, st->lpcSize);
307    st->lpc[0]=1;
308
309    /* LPC to LSPs (x-domain) transform */
310    roots=lpc_to_lsp (st->lpc, st->lpcSize, st->lsp, 6, 0.002, st->stack);
311    if (roots!=st->lpcSize)
312    {
313       fprintf (stderr, "roots!=st->lpcSize (found only %d roots)\n", roots);
314       exit(1);
315    }
316
317    {
318       for (i=0;i<st->frame_size;i++)
319          st->buf[i] = st->st_low.exc[i] * st->exc_window[i];
320       
321       /* Compute auto-correlation */
322       autocorr(st->buf, st->autocorr, st->lpcSize+1, st->frame_size);
323       
324       st->autocorr[0] += 1;        /* prevents NANs */
325       st->autocorr[0] *= st->lpc_floor; /* Noise floor in auto-correlation domain */
326       /* Lag windowing: equivalent to filtering in the power-spectrum domain */
327       for (i=0;i<st->lpcSize+1;i++)
328          st->autocorr[i] *= st->lagWindow[i];
329       
330       /* Levinson-Durbin */
331       wld(st->lpc+1, st->autocorr, st->rc, st->lpcSize);
332       st->lpc[0]=1;
333       printf ("exc_lpc: ");
334       for(i=0;i<=st->lpcSize;i++)
335          printf ("%f ", st->lpc[i]);
336       printf ("\n");
337       residue_mem(st->st_low.exc, st->lpc, st->exc_alias, st->frame_size, st->lpcSize, st->mem_exc);
338    }
339
340
341    /* x-domain to angle domain*/
342    for (i=0;i<st->lpcSize;i++)
343       st->lsp[i] = acos(st->lsp[i]);
344
345    /* LSP quantization */
346    lsp_quant_high(st->lsp, st->qlsp, st->lpcSize, bits);
347    
348    /*printf ("high_lsp:");
349    for (i=0;i<st->lpcSize;i++)
350       printf (" %f", st->lsp[i]);
351       printf ("\n");
352    for (i=0;i<st->lpcSize;i++)
353    st->qlsp[i]=st->lsp[i];
354    */
355
356    if (st->first)
357    {
358       for (i=0;i<st->lpcSize;i++)
359          st->old_lsp[i] = st->lsp[i];
360       for (i=0;i<st->lpcSize;i++)
361          st->old_qlsp[i] = st->qlsp[i];
362    }
363    
364    for (sub=0;sub<st->nbSubframes;sub++)
365    {
366       float *exc, *sp, *mem, *res, *target, *sw, tmp, filter_ratio;
367       int offset;
368       
369       offset = st->subframeSize*sub;
370       sp=st->high+offset;
371       exc=st->excBuf+offset;
372       res=st->res+offset;
373       target=st->target+offset;
374       sw=st->sw+offset;
375
376       mem=PUSH(st->stack, st->lpcSize);
377       
378       /* LSP interpolation (quantized and unquantized) */
379       tmp = (.5 + sub)/st->nbSubframes;
380       for (i=0;i<st->lpcSize;i++)
381          st->interp_lsp[i] = (1-tmp)*st->old_lsp[i] + tmp*st->lsp[i];
382       for (i=0;i<st->lpcSize;i++)
383          st->interp_qlsp[i] = (1-tmp)*st->old_qlsp[i] + tmp*st->qlsp[i];
384
385       /* Compute interpolated LPCs (quantized and unquantized) */
386       for (i=0;i<st->lpcSize;i++)
387          st->interp_lsp[i] = cos(st->interp_lsp[i]);
388       lsp_to_lpc(st->interp_lsp, st->interp_lpc, st->lpcSize,st->stack);
389
390       for (i=0;i<st->lpcSize;i++)
391          st->interp_qlsp[i] = cos(st->interp_qlsp[i]);
392       lsp_to_lpc(st->interp_qlsp, st->interp_qlpc, st->lpcSize, st->stack);
393
394       bw_lpc(st->gamma1, st->interp_lpc, st->bw_lpc1, st->lpcSize);
395       bw_lpc(st->gamma2, st->interp_lpc, st->bw_lpc2, st->lpcSize);
396
397       {
398          float tmp, rl=0, rh=0;
399          tmp=1;
400          for (i=0;i<=st->lpcSize;i++)
401          {
402             rh += tmp*st->interp_lpc[i];
403             tmp = -tmp;
404          }
405          rl = st->st_low.pi_gain[sub];
406          rl=1/(fabs(rl)+.001);
407          rh=1/(fabs(rh)+.001);
408          filter_ratio=fabs(.001+rh)/(.001+fabs(rl));
409          printf ("filter_ratio: %f\n", filter_ratio);
410       }
411 #if 0 /* 1 for spectral folding excitation, 0 for stochastic */
412       for (i=0;i<st->lpcSize;i++)
413          mem[i]=st->mem_sp[i];
414       residue_mem(sp, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, st->mem_sp);
415       {
416          float el=0,eh=0,g;
417          for (i=0;i<st->subframeSize;i++)
418             eh+=sqr(exc[i]);
419          /*for (i=0;i<st->subframeSize;i++)
420             el+=sqr(st->exc_alias[offset+i]);*/
421          for (i=0;i<st->subframeSize;i++)
422            el+=sqr(st->st_low.exc[offset+i]);
423          g=eh/(.01+el);
424          printf ("ratio: %f\n", g);
425          g=sqrt(g);
426          for (i=0;i<st->subframeSize;i++)
427            exc[i]=g*st->st_low.exc[offset+i];
428          /*for (i=0;i<st->subframeSize;i++)
429            exc[i]=g*st->exc_alias[offset+i];*/
430       }
431       syn_filt_mem(exc, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, mem);
432
433 #else
434       /* Reset excitation */
435       for (i=0;i<st->subframeSize;i++)
436          exc[i]=0;
437
438       /* Compute zero response of A(z/g1) / ( A(z/g2) * Aq(z) ) */
439       for (i=0;i<st->lpcSize;i++)
440          mem[i]=st->mem_sp[i];
441       syn_filt_mem(exc, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, mem);
442       for (i=0;i<st->lpcSize;i++)
443          mem[i]=st->mem_sp[i];
444       residue_mem(exc, st->bw_lpc1, res, st->subframeSize, st->lpcSize, mem);
445       for (i=0;i<st->lpcSize;i++)
446          mem[i]=st->mem_sw[i];
447       syn_filt_mem(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize, mem);
448
449       /* Compute weighted signal */
450       for (i=0;i<st->lpcSize;i++)
451          mem[i]=st->mem_sp[i];
452       residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, mem);
453       for (i=0;i<st->lpcSize;i++)
454          mem[i]=st->mem_sw[i];
455       syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, mem);
456
457       /* Compute target signal */
458       for (i=0;i<st->subframeSize;i++)
459          target[i]=sw[i]-res[i];
460       {
461          int ind;
462          float gain;
463 #if 0
464          
465          float el=0,eh=0,g;
466          residue_mem(sp, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, st->mem_sp2);
467          
468          for (i=0;i<st->subframeSize;i++)
469             eh+=sqr(exc[i]);
470          overlap_cb_search(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
471                            &stoc[0], 512, &gain, &ind, st->lpcSize,
472                            st->subframeSize);
473          for (i=0;i<st->subframeSize;i++)
474             exc[i]=gain*stoc[ind+i];
475          for (i=0;i<st->subframeSize;i++)
476             el+=sqr(exc[i]);
477          g=sqrt(eh/(el+.001));
478          for (i=0;i<st->subframeSize;i++)
479             exc[i]*=g;
480          
481 #else
482          int k,N=2;
483          float el=0,eh=0,g;
484          int *index;
485          float *gains;
486          gains = PUSH(st->stack, N);
487          index = (int*) PUSH(st->stack, N);
488          residue_mem(sp, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, st->mem_sp2);
489          
490          for (i=0;i<st->subframeSize;i++)
491             eh+=sqr(exc[i]);
492
493          for (i=0;i<st->subframeSize;i++)
494             exc[i]=0;
495          for (k=0;k<N;k++)
496          {
497             int of=k*st->subframeSize/N;
498          overlap_cb_search(target+of, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
499                            &stoc[0], 128, &gains[k], &index[k], st->lpcSize,
500                            st->subframeSize/N);
501          for (i=0;i<st->subframeSize;i++)
502             res[i]=0;
503          for (i=0;i<st->subframeSize/N;i++)
504             res[of+i]=gains[k]*stoc[index[k]+i];
505          residue_zero(res, st->bw_lpc1, res, st->subframeSize, st->lpcSize);
506          syn_filt_zero(res, st->interp_qlpc, res, st->subframeSize, st->lpcSize);
507          syn_filt_zero(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize);
508          for (i=0;i<st->subframeSize;i++)
509             target[i]-=res[i];
510          for (i=0;i<st->subframeSize/N;i++)
511             exc[of+i]+=gains[k]*stoc[index[k]+i];
512          }
513          for (i=0;i<st->subframeSize;i++)
514             el+=sqr(exc[i]);
515          g=sqrt(eh/(el+.001));
516
517          for (i=0;i<st->subframeSize;i++)
518             el+=sqr(st->st_low.exc[offset+i]);
519
520          for (k=0;k<N;k++)
521          {
522             int sign=0;
523             float quant;
524             int of=k*st->subframeSize/N;
525             gains[k]*=g;
526 #if 1
527             if (gains[k]<0)
528             {
529                sign=1;
530                gains[k] = -gains[k];
531             }
532             quant = (1+gains[k])*filter_ratio/(1+sqrt(el/st->subframeSize));
533             if (quant>10)
534               quant=10;
535             if (quant<.0001)
536               quant=.1;
537             gains[k]=quant*(1+sqrt(el/st->subframeSize))/filter_ratio;
538             if (sign)
539                gains[k] = -gains[k];
540             printf ("quant_high_gain: %f\n", quant);
541 #endif
542             for (i=0;i<st->subframeSize/N;i++)
543                exc[of+i]=gains[k]*stoc[index[k]+i];
544          }
545          /*for (i=0;i<st->subframeSize;i++)
546            exc[i]*=g;*/
547          POP(st->stack);
548          POP(st->stack);
549 #endif
550       }
551
552       /*Keep the previous memory*/
553       for (i=0;i<st->lpcSize;i++)
554          mem[i]=st->mem_sp[i];
555       /* Final signal synthesis from excitation */
556       syn_filt_mem(exc, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, st->mem_sp);
557        
558       /* Compute weighted signal again, from synthesized speech (not sure it's the right thing) */
559       residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, mem);
560       syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, st->mem_sw);
561 #endif
562
563       POP(st->stack);
564    }
565
566    /* Up-sample coded low-band and high-band*/
567    for (i=0;i<st->frame_size;i++)
568    {
569       st->x0[(i<<1)]=st->x0d[i];
570       st->x1[(i<<1)]=st->high[i];
571       st->x0[(i<<1)+1]=0;
572       st->x1[(i<<1)+1]=0;
573    }
574    /* Reconstruct the original */
575    fir_mem(st->x0, h0, st->y0, st->full_frame_size, QMF_ORDER, st->g0_mem);
576    fir_mem(st->x1, h1, st->y1, st->full_frame_size, QMF_ORDER, st->g1_mem);
577    for (i=0;i<st->full_frame_size;i++)
578       in[i]=2*(st->y0[i]-st->y1[i]);
579
580    for (i=0;i<st->lpcSize;i++)
581       st->old_lsp[i] = st->lsp[i];
582    for (i=0;i<st->lpcSize;i++)
583       st->old_qlsp[i] = st->qlsp[i];
584
585    st->first=0;
586 }