3c83ec077325f35b11501bb3766fcd4582e9f234
[speexdsp.git] / libspeex / speex.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 #include <stdlib.h>
20 #include <stdio.h>
21 #include <string.h>
22 #include <math.h>
23 #include "speex.h"
24 #include "lpc.h"
25 #include "lsp.h"
26 #include "ltp.h"
27 #include "quant_lsp.h"
28 #include "cb_search.h"
29 #include "filters.h"
30 #include "stack_alloc.h"
31
32 #ifndef M_PI
33 #define M_PI           3.14159265358979323846  /* pi */
34 #endif
35
36 #define sqr(x) ((x)*(x))
37 #define min(a,b) ((a) < (b) ? (a) : (b))
38
39 void *nb_encoder_init(SpeexMode *m)
40 {
41    EncState *st;
42    SpeexNBMode *mode;
43    int i;
44    float tmp;
45
46    mode=m->mode;
47    st = malloc(sizeof(EncState));
48    st->mode=m;
49    /* Codec parameters, should eventually have several "modes"*/
50    st->frameSize = mode->frameSize;
51    st->windowSize = mode->windowSize;
52    st->nbSubframes=mode->frameSize/mode->subframeSize;
53    st->subframeSize=mode->subframeSize;
54    st->lpcSize = mode->lpcSize;
55    st->bufSize = mode->bufSize;
56    st->gamma1=mode->gamma1;
57    st->gamma2=mode->gamma2;
58    st->min_pitch=mode->pitchStart;
59    st->max_pitch=mode->pitchEnd;
60    st->lag_factor=mode->lag_factor;
61    st->lpc_floor = mode->lpc_floor;
62    st->preemph = mode->preemph;
63   
64
65    st->lsp_quant = mode->lsp_quant;
66    st->ltp_quant = mode->ltp_quant;
67    st->ltp_params = mode->ltp_params;
68    st->innovation_quant = mode->innovation_quant;
69    st->innovation_params = mode->innovation_params;
70
71    st->pre_mem=0;
72    st->pre_mem2=0;
73
74    /* Over-sampling filter (fractional pitch)*/
75    st->os_fact=4;
76    st->os_filt_ord2=4*st->os_fact;
77    st->os_filt = malloc((1+2*st->os_filt_ord2)*sizeof(float));
78    st->os_filt[st->os_filt_ord2] = 1;
79    for (i=1;i<=st->os_filt_ord2;i++)
80    {
81       float x=M_PI*i/st->os_fact;
82       st->os_filt[st->os_filt_ord2-i] = st->os_filt[st->os_filt_ord2+i]=sin(x)/x*(.5+.5*cos(M_PI*i/st->os_filt_ord2));
83    }
84    /* Normalizing the over-sampling filter */
85    tmp=0;
86    for (i=0;i<2*st->os_filt_ord2+1;i++)
87       tmp += st->os_filt[i];
88    tmp=1/tmp;
89    for (i=0;i<2*st->os_filt_ord2+1;i++)
90       st->os_filt[i] *= tmp;
91
92    /*for (i=0;i<2*st->os_filt_ord2+1;i++)
93       printf ("%f ", st->os_filt[i]);
94       printf ("\n");*/
95
96    /* Allocating input buffer */
97    st->inBuf = calloc(st->bufSize,sizeof(float));
98    st->frame = st->inBuf + st->bufSize - st->windowSize;
99    /* Allocating excitation buffer */
100    st->excBuf = calloc(st->bufSize,sizeof(float));
101    st->exc = st->excBuf + st->bufSize - st->windowSize;
102    st->swBuf = calloc(st->bufSize,sizeof(float));
103    st->sw = st->swBuf + st->bufSize - st->windowSize;
104
105    st->exc2Buf = calloc(st->bufSize,sizeof(float));
106    st->exc2 = st->exc2Buf + st->bufSize - st->windowSize;
107
108    /* Hanning window */
109    st->window = malloc(st->windowSize*sizeof(float));
110    for (i=0;i<st->windowSize;i++)
111       st->window[i]=.5*(1-cos(2*M_PI*i/st->windowSize));
112
113    /* Create the window for autocorrelation (lag-windowing) */
114    st->lagWindow = malloc((st->lpcSize+1)*sizeof(float));
115    for (i=0;i<st->lpcSize+1;i++)
116       st->lagWindow[i]=exp(-.5*sqr(2*M_PI*st->lag_factor*i));
117
118    st->autocorr = malloc((st->lpcSize+1)*sizeof(float));
119
120    st->stack = calloc(20000, sizeof(float));
121
122    st->buf2 = malloc(st->windowSize*sizeof(float));
123
124    st->lpc = malloc((st->lpcSize+1)*sizeof(float));
125    st->interp_lpc = malloc((st->lpcSize+1)*sizeof(float));
126    st->interp_qlpc = malloc((st->lpcSize+1)*sizeof(float));
127    st->bw_lpc1 = malloc((st->lpcSize+1)*sizeof(float));
128    st->bw_lpc2 = malloc((st->lpcSize+1)*sizeof(float));
129
130    st->lsp = malloc(st->lpcSize*sizeof(float));
131    st->qlsp = malloc(st->lpcSize*sizeof(float));
132    st->old_lsp = malloc(st->lpcSize*sizeof(float));
133    st->old_qlsp = malloc(st->lpcSize*sizeof(float));
134    st->interp_lsp = malloc(st->lpcSize*sizeof(float));
135    st->interp_qlsp = malloc(st->lpcSize*sizeof(float));
136    st->rc = malloc(st->lpcSize*sizeof(float));
137    st->first = 1;
138    
139    st->mem_sp = calloc(st->lpcSize, sizeof(float));
140    st->mem_sw = calloc(st->lpcSize, sizeof(float));
141
142    st->pi_gain = calloc(st->nbSubframes, sizeof(float));
143
144    st->pitch = calloc(st->nbSubframes, sizeof(int));
145    return st;
146 }
147
148 void nb_encoder_destroy(void *state)
149 {
150    EncState *st=state;
151    /* Free all allocated memory */
152    free(st->inBuf);
153    free(st->excBuf);
154    free(st->swBuf);
155    free(st->os_filt);
156    free(st->exc2Buf);
157    free(st->stack);
158
159    free(st->window);
160    free(st->buf2);
161    free(st->lpc);
162    free(st->interp_lpc);
163    free(st->interp_qlpc);
164    
165    free(st->bw_lpc1);
166    free(st->bw_lpc2);
167    free(st->autocorr);
168    free(st->lagWindow);
169    free(st->lsp);
170    free(st->qlsp);
171    free(st->old_lsp);
172    free(st->interp_lsp);
173    free(st->old_qlsp);
174    free(st->interp_qlsp);
175    free(st->rc);
176
177    free(st->mem_sp);
178    free(st->mem_sw);
179    free(st->pi_gain);
180    free(st->pitch);
181    
182    free(st);
183 }
184
185
186 void nb_encode(void *state, float *in, FrameBits *bits)
187 {
188    EncState *st;
189    int i, sub, roots;
190    float error;
191    
192    st=state;
193    /* Copy new data in input buffer */
194    memmove(st->inBuf, st->inBuf+st->frameSize, (st->bufSize-st->frameSize)*sizeof(float));
195    st->inBuf[st->bufSize-st->frameSize] = in[0] - st->preemph*st->pre_mem;
196    for (i=1;i<st->frameSize;i++)
197       st->inBuf[st->bufSize-st->frameSize+i] = in[i] - st->preemph*in[i-1];
198    st->pre_mem = in[st->frameSize-1];
199
200    memmove(st->exc2Buf, st->exc2Buf+st->frameSize, (st->bufSize-st->frameSize)*sizeof(float));
201    memmove(st->excBuf, st->excBuf+st->frameSize, (st->bufSize-st->frameSize)*sizeof(float));
202    memmove(st->swBuf, st->swBuf+st->frameSize, (st->bufSize-st->frameSize)*sizeof(float));
203
204    /* Window for analysis */
205    for (i=0;i<st->windowSize;i++)
206       st->buf2[i] = st->frame[i] * st->window[i];
207
208    /* Compute auto-correlation */
209    autocorr(st->buf2, st->autocorr, st->lpcSize+1, st->windowSize);
210
211    st->autocorr[0] += 1;        /* prevents NANs */
212    st->autocorr[0] *= st->lpc_floor; /* Noise floor in auto-correlation domain */
213    /* Lag windowing: equivalent to filtering in the power-spectrum domain */
214    for (i=0;i<st->lpcSize+1;i++)
215       st->autocorr[i] *= st->lagWindow[i];
216
217    /* Levinson-Durbin */
218    error = wld(st->lpc+1, st->autocorr, st->rc, st->lpcSize);
219    st->lpc[0]=1;
220
221    /* LPC to LSPs (x-domain) transform */
222    roots=lpc_to_lsp (st->lpc, st->lpcSize, st->lsp, 6, 0.002, st->stack);
223    if (roots!=st->lpcSize)
224    {
225       fprintf (stderr, "roots!=st->lpcSize (found only %d roots)\n", roots);
226       exit(1);
227    }
228
229    /* x-domain to angle domain*/
230    for (i=0;i<st->lpcSize;i++)
231       st->lsp[i] = acos(st->lsp[i]);
232    
233    /* LSP Quantization */
234    st->lsp_quant(st->lsp, st->qlsp, st->lpcSize, bits);
235    /*for (i=0;i<st->lpcSize;i++)
236      st->qlsp[i]=st->lsp[i];*/
237    /*printf ("LSP ");
238    for (i=0;i<st->lpcSize;i++)
239       printf ("%f ", st->lsp[i]);
240    printf ("\n");
241    printf ("QLSP ");
242    for (i=0;i<st->lpcSize;i++)
243       printf ("%f ", st->qlsp[i]);
244    printf ("\n");*/
245    /* Special case for first frame */
246    if (st->first)
247    {
248       for (i=0;i<st->lpcSize;i++)
249          st->old_lsp[i] = st->lsp[i];
250       for (i=0;i<st->lpcSize;i++)
251          st->old_qlsp[i] = st->qlsp[i];
252    }
253
254    /* Loop on sub-frames */
255    for (sub=0;sub<st->nbSubframes;sub++)
256    {
257       float esig, enoise, snr, tmp;
258       int   offset;
259       float *sp, *sw, *res, *exc, *target, *mem, *exc2;
260       int pitch;
261
262       /* Offset relative to start of frame */
263       offset = st->subframeSize*sub;
264       /* Original signal */
265       sp=st->frame+offset;
266       /* Excitation */
267       exc=st->exc+offset;
268       /* Weighted signal */
269       sw=st->sw+offset;
270
271       exc2=st->exc2+offset;
272
273       /* Filter response */
274       res = PUSH(st->stack, st->subframeSize);
275       /* Target signal */
276       target = PUSH(st->stack, st->subframeSize);
277       mem = PUSH(st->stack, st->lpcSize);
278
279       /* LSP interpolation (quantized and unquantized) */
280       tmp = (.5 + sub)/st->nbSubframes;
281       for (i=0;i<st->lpcSize;i++)
282          st->interp_lsp[i] = (1-tmp)*st->old_lsp[i] + tmp*st->lsp[i];
283       for (i=0;i<st->lpcSize;i++)
284          st->interp_qlsp[i] = (1-tmp)*st->old_qlsp[i] + tmp*st->qlsp[i];
285
286       /* Compute interpolated LPCs (quantized and unquantized) */
287       for (i=0;i<st->lpcSize;i++)
288          st->interp_lsp[i] = cos(st->interp_lsp[i]);
289       lsp_to_lpc(st->interp_lsp, st->interp_lpc, st->lpcSize,st->stack);
290
291       for (i=0;i<st->lpcSize;i++)
292          st->interp_qlsp[i] = cos(st->interp_qlsp[i]);
293       lsp_to_lpc(st->interp_qlsp, st->interp_qlpc, st->lpcSize, st->stack);
294
295       tmp=1;
296       st->pi_gain[sub]=0;
297       for (i=0;i<=st->lpcSize;i++)
298       {
299          st->pi_gain[sub] += tmp*st->interp_qlpc[i];
300          tmp = -tmp;
301       }
302      
303
304       /* Compute bandwidth-expanded (unquantized) LPCs for perceptual weighting */
305       bw_lpc(st->gamma1, st->interp_lpc, st->bw_lpc1, st->lpcSize);
306       if (st->gamma2>=0)
307          bw_lpc(st->gamma2, st->interp_lpc, st->bw_lpc2, st->lpcSize);
308       else
309       {
310          st->bw_lpc2[0]=1;
311          st->bw_lpc2[1]=-st->preemph;
312          for (i=2;i<=st->lpcSize;i++)
313             st->bw_lpc2[i]=0;
314       }
315 #ifdef DEBUG
316       printf ("\nlpc0 ");
317       for (i=0;i<=st->lpcSize;i++)
318          printf ("%f ", st->interp_lpc[i]);
319       printf ("\nlpc1 ");
320       for (i=0;i<=st->lpcSize;i++)
321          printf ("%f ", st->bw_lpc1[i]);
322       printf ("\nlpc2 ");
323       for (i=0;i<=st->lpcSize;i++)
324          printf ("%f ", st->bw_lpc2[i]);
325       printf ("\n\n");
326 #endif
327       /* Reset excitation */
328       for (i=0;i<st->subframeSize;i++)
329          exc[i]=0;
330       for (i=0;i<st->subframeSize;i++)
331          exc2[i]=0;
332
333       /* Compute zero response of A(z/g1) / ( A(z/g2) * Aq(z) ) */
334       for (i=0;i<st->lpcSize;i++)
335          mem[i]=st->mem_sp[i];
336       syn_filt_mem(exc, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, mem);
337       for (i=0;i<st->lpcSize;i++)
338          mem[i]=st->mem_sp[i];
339       residue_mem(exc, st->bw_lpc1, res, st->subframeSize, st->lpcSize, mem);
340       for (i=0;i<st->lpcSize;i++)
341          mem[i]=st->mem_sw[i];
342       syn_filt_mem(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize, mem);
343
344       /* Compute weighted signal */
345       for (i=0;i<st->lpcSize;i++)
346          mem[i]=st->mem_sp[i];
347       residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, mem);
348       for (i=0;i<st->lpcSize;i++)
349          mem[i]=st->mem_sw[i];
350       syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, mem);
351       
352       esig=0;
353       for (i=0;i<st->subframeSize;i++)
354          esig+=sw[i]*sw[i];
355       
356       /* Compute target signal */
357       for (i=0;i<st->subframeSize;i++)
358          target[i]=sw[i]-res[i];
359
360       for (i=0;i<st->subframeSize;i++)
361          exc[i]=0;
362
363       /* Long-term prediction */
364       pitch = st->ltp_quant(target, sw, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
365                     exc, st->ltp_params, st->min_pitch, st->max_pitch, 
366                     st->lpcSize, st->subframeSize, bits, st->stack, exc2);
367       st->pitch[sub]=pitch;
368
369       /* Update target for adaptive codebook contribution */
370       residue_zero(exc, st->bw_lpc1, res, st->subframeSize, st->lpcSize);
371       syn_filt_zero(res, st->interp_qlpc, res, st->subframeSize, st->lpcSize);
372       syn_filt_zero(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize);
373       for (i=0;i<st->subframeSize;i++)
374         target[i]-=res[i];
375
376       /* Compute noise energy and SNR */
377       enoise=0;
378       for (i=0;i<st->subframeSize;i++)
379          enoise += target[i]*target[i];
380       snr = 10*log10((esig+1)/(enoise+1));
381       /*st->pitch[sub]=(int)snr;*/
382 #ifdef DEBUG
383       printf ("pitch SNR = %f\n", snr);
384 #endif
385
386
387 #if 0 /*If set to 1, compute "real innovation" i.e. cheat to get perfect reconstruction*/
388       syn_filt_zero(target, st->bw_lpc1, res, st->subframeSize, st->lpcSize);
389       residue_zero(res, st->interp_qlpc, st->buf2, st->subframeSize, st->lpcSize);
390       residue_zero(st->buf2, st->bw_lpc2, st->buf2, st->subframeSize, st->lpcSize);
391       if (1||(snr>9 && (rand()%6==0)))
392       {
393          float ener=0;
394          printf ("exc ");
395          for (i=0;i<st->subframeSize;i++)
396          {
397             ener+=st->buf2[i]*st->buf2[i];
398             if (i && i%5==0)
399                printf ("\nexc ");
400             printf ("%f ", st->buf2[i]);
401          }
402          printf ("\n");
403       printf ("innovation_energy = %f\n", ener);
404       }
405       for (i=0;i<st->subframeSize;i++)
406          exc[i]+=st->buf2[i];
407 #else
408       /* Perform a split-codebook search */
409       st->innovation_quant(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,
410                            st->innovation_params, st->lpcSize,
411                            st->subframeSize, exc, bits, st->stack);
412
413 #endif
414       /* Compute weighted noise energy and SNR */
415       enoise=0;
416       for (i=0;i<st->subframeSize;i++)
417          enoise += target[i]*target[i];
418       snr = 10*log10((esig+1)/(enoise+1));
419 #ifdef DEBUG
420       printf ("seg SNR = %f\n", snr);
421 #endif
422       /*Keep the previous memory*/
423       for (i=0;i<st->lpcSize;i++)
424          mem[i]=st->mem_sp[i];
425       /* Final signal synthesis from excitation */
426       syn_filt_mem(exc, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, st->mem_sp);
427
428       /* Compute weighted signal again, from synthesized speech (not sure it's the right thing) */
429       residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, mem);
430       syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, st->mem_sw);
431
432 #if 0
433       /*for (i=0;i<st->subframeSize;i++)
434         exc2[i]=.75*exc[i]+.2*exc[i-pitch]+.05*exc[i-2*pitch];*/
435       {
436          float max_exc=0;
437          for (i=0;i<st->subframeSize;i++)
438             if (fabs(exc[i])>max_exc)
439                max_exc=fabs(exc[i]);
440          max_exc=1/(max_exc+.01);
441          for (i=0;i<st->subframeSize;i++)
442          {
443             float xx=max_exc*exc[i];
444             exc2[i]=exc[i]*(1-exp(-100*xx*xx));
445          }
446       }
447 #else
448       for (i=0;i<st->subframeSize;i++)
449          exc2[i]=exc[i];
450 #endif
451       POP(st->stack);
452       POP(st->stack);
453       POP(st->stack);
454    }
455
456    /* Store the LSPs for interpolation in the next frame */
457    for (i=0;i<st->lpcSize;i++)
458       st->old_lsp[i] = st->lsp[i];
459    for (i=0;i<st->lpcSize;i++)
460       st->old_qlsp[i] = st->qlsp[i];
461
462    /* The next frame will not be the first (Duh!) */
463    st->first = 0;
464
465    /* Replace input by synthesized speech */
466    in[0] = st->frame[0] + st->preemph*st->pre_mem2;
467    for (i=1;i<st->frameSize;i++)
468      in[i]=st->frame[i] + st->preemph*in[i-1];
469    st->pre_mem2=in[st->frameSize-1];
470 }
471
472
473 void *nb_decoder_init(SpeexMode *m)
474 {
475    DecState *st;
476    SpeexNBMode *mode;
477    int i;
478
479    mode=m->mode;
480    st = malloc(sizeof(DecState));
481    st->mode=m;
482
483    st->first=1;
484    /* Codec parameters, should eventually have several "modes"*/
485    st->frameSize = mode->frameSize;
486    st->windowSize = mode->windowSize;
487    st->nbSubframes=mode->frameSize/mode->subframeSize;
488    st->subframeSize=mode->subframeSize;
489    st->lpcSize = mode->lpcSize;
490    st->bufSize = mode->bufSize;
491    st->gamma1=mode->gamma1;
492    st->gamma2=mode->gamma2;
493    st->min_pitch=mode->pitchStart;
494    st->max_pitch=mode->pitchEnd;
495    st->preemph = mode->preemph;
496
497    st->pre_mem=0;
498    st->lsp_unquant = mode->lsp_unquant;
499    st->ltp_unquant = mode->ltp_unquant;
500    st->ltp_params = mode->ltp_params;
501
502    st->innovation_unquant = mode->innovation_unquant;
503    st->innovation_params = mode->innovation_params;
504
505    st->stack = calloc(10000, sizeof(float));
506
507    st->inBuf = malloc(st->bufSize*sizeof(float));
508    st->frame = st->inBuf + st->bufSize - st->windowSize;
509    st->excBuf = malloc(st->bufSize*sizeof(float));
510    st->exc = st->excBuf + st->bufSize - st->windowSize;
511    for (i=0;i<st->bufSize;i++)
512       st->inBuf[i]=0;
513    for (i=0;i<st->bufSize;i++)
514       st->excBuf[i]=0;
515
516    st->interp_qlpc = malloc((st->lpcSize+1)*sizeof(float));
517    st->qlsp = malloc(st->lpcSize*sizeof(float));
518    st->old_qlsp = malloc(st->lpcSize*sizeof(float));
519    st->interp_qlsp = malloc(st->lpcSize*sizeof(float));
520    st->mem_sp = calloc(st->lpcSize, sizeof(float));
521
522    st->pi_gain = calloc(st->nbSubframes, sizeof(float));
523    
524    return st;
525 }
526
527 void nb_decoder_destroy(void *state)
528 {
529    DecState *st;
530    st=state;
531    free(st->inBuf);
532    free(st->excBuf);
533    free(st->interp_qlpc);
534    free(st->qlsp);
535    free(st->old_qlsp);
536    free(st->interp_qlsp);
537    free(st->stack);
538    free(st->mem_sp);
539    free(st->pi_gain);
540    
541    free(state);
542 }
543
544 void nb_decode(void *state, FrameBits *bits, float *out)
545 {
546    DecState *st;
547    int i, sub;
548
549    st=state;
550
551    memmove(st->inBuf, st->inBuf+st->frameSize, (st->bufSize-st->frameSize)*sizeof(float));
552    memmove(st->excBuf, st->excBuf+st->frameSize, (st->bufSize-st->frameSize)*sizeof(float));
553
554
555    st->lsp_unquant(st->qlsp, st->lpcSize, bits);
556    if (st->first)
557    {
558       for (i=0;i<st->lpcSize;i++)
559          st->old_qlsp[i] = st->qlsp[i];
560    }
561
562    /*Loop on subframes */
563    for (sub=0;sub<st->nbSubframes;sub++)
564    {
565       int offset;
566       float *sp, *exc, tmp;
567       
568       /* Offset relative to start of frame */
569       offset = st->subframeSize*sub;
570       /* Original signal */
571       sp=st->frame+offset;
572       /* Excitation */
573       exc=st->exc+offset;
574
575       /* LSP interpolation (quantized and unquantized) */
576       tmp = (.5 + sub)/st->nbSubframes;
577       for (i=0;i<st->lpcSize;i++)
578          st->interp_qlsp[i] = (1-tmp)*st->old_qlsp[i] + tmp*st->qlsp[i];
579
580       /* Compute interpolated LPCs (unquantized) */
581       for (i=0;i<st->lpcSize;i++)
582          st->interp_qlsp[i] = cos(st->interp_qlsp[i]);
583       lsp_to_lpc(st->interp_qlsp, st->interp_qlpc, st->lpcSize, st->stack);
584
585       tmp=1;
586       st->pi_gain[sub]=0;
587       for (i=0;i<=st->lpcSize;i++)
588       {
589          st->pi_gain[sub] += tmp*st->interp_qlpc[i];
590          tmp = -tmp;
591       }
592
593       /* Reset excitation */
594       for (i=0;i<st->subframeSize;i++)
595          exc[i]=0;
596
597       /*Adaptive codebook contribution*/
598       st->ltp_unquant(exc, st->min_pitch, st->max_pitch, st->ltp_params, st->subframeSize, bits, st->stack);
599        
600       /*Fixed codebook contribution*/
601       st->innovation_unquant(exc, st->innovation_params, st->subframeSize, bits, st->stack);
602
603       /*Compute decoded signal*/
604       syn_filt_mem(exc, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, st->mem_sp);
605
606    }
607    
608    /*Copy output signal*/
609    for (i=0;i<st->frameSize;i++)
610       out[i]=st->frame[i];
611
612    out[0] = st->frame[0] + st->preemph*st->pre_mem;
613    for (i=1;i<st->frameSize;i++)
614      out[i]=st->frame[i] + st->preemph*out[i-1];
615    st->pre_mem=out[st->frameSize-1];
616
617
618    /* Store the LSPs for interpolation in the next frame */
619    for (i=0;i<st->lpcSize;i++)
620       st->old_qlsp[i] = st->qlsp[i];
621
622    /* The next frame will not be the first (Duh!) */
623    st->first = 0;
624
625 }