Multi-pulse seems to work not too bad (but disabled by default)
[speexdsp.git] / libspeex / cb_search.c
1 /*-----------------------------------------------------------------------*\
2
3     FILE........: GAINSHAPE.C
4     TYPE........: C Module
5     AUTHOR......: David Rowe
6     COMPANY.....: Voicetronix
7     DATE CREATED: 19/2/02
8
9     General gain-shape codebook search.
10
11 \*-----------------------------------------------------------------------*/
12
13 /* Modified by Jean-Marc Valin 2002
14
15    This library is free software; you can redistribute it and/or
16    modify it under the terms of the GNU Lesser General Public
17    License as published by the Free Software Foundation; either
18    version 2.1 of the License, or (at your option) any later version.
19    
20    This library is distributed in the hope that it will be useful,
21    but WITHOUT ANY WARRANTY; without even the implied warranty of
22    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
23    Lesser General Public License for more details.
24    
25    You should have received a copy of the GNU Lesser General Public
26    License along with this library; if not, write to the Free Software
27    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
28 */
29
30
31
32 #include <stdlib.h>
33 #include <cb_search.h>
34 #include "filters.h"
35 #include <math.h>
36 #include <stdio.h>
37 #include "stack_alloc.h"
38 #include "vq.h"
39 #include "matrix.h"
40
41 extern float exc_gains_wb2_table[];
42 /*---------------------------------------------------------------------------*\
43                                                                              
44  void overlap_cb_search()                                                             
45                                                                               
46  Searches a gain/shape codebook consisting of overlapping entries for the    
47  closest vector to the target.  Gives identical results to search() above   
48  buts uses fast end correction algorithm for the synthesis of response       
49  vectors.                                                                     
50                                                                              
51 \*---------------------------------------------------------------------------*/
52
53 float overlap_cb_search(
54 float target[],                 /* target vector */
55 float ak[],                     /* LPCs for this subframe */
56 float awk1[],                   /* Weighted LPCs for this subframe */
57 float awk2[],                   /* Weighted LPCs for this subframe */
58 float codebook[],               /* overlapping codebook */
59 int   entries,                  /* number of overlapping entries to search */
60 float *gain,                    /* gain of optimum entry */
61 int   *index,                   /* index of optimum entry */
62 int   p,                        /* number of LPC coeffs */
63 int   nsf                       /* number of samples in subframe */
64 )
65 {
66   float *resp;                  /* zero state response to current entry */
67   float *h;                     /* impulse response of synthesis filter */
68   float *impulse;               /* excitation vector containing one impulse */
69   float d,e,g,score;            /* codebook searching variables */
70   float bscore;                 /* score of "best" vector so far */
71   int i,k;                      /* loop variables */
72
73   /* Initialise */
74   
75   resp = (float*)malloc(sizeof(float)*nsf);
76   h = (float*)malloc(sizeof(float)*nsf);
77   impulse = (float*)malloc(sizeof(float)*nsf);
78
79   for(i=0; i<nsf; i++)
80     impulse[i] = 0.0;
81    
82   *gain = 0.0;
83   *index = 0;
84   bscore = 0.0;
85   impulse[0] = 1.0;
86
87   /* Calculate impulse response of  A(z/g2) / ( A(z)*(z/g1) ) */
88   residue_zero(impulse, awk1, h, nsf, p);
89   syn_filt_zero(h, ak, h, nsf, p);
90   syn_filt_zero(h, awk2, h, nsf,p);
91   
92   /* Calculate codebook zero-response */
93   residue_zero(&codebook[entries-1],awk1,resp,nsf,p);
94   syn_filt_zero(resp,ak,resp,nsf,p);
95   syn_filt_zero(resp,awk2,resp,nsf,p);
96     
97   /* Search codebook backwards using end correction for synthesis */
98   
99   for(k=entries-1; k>=0; k--) {
100
101     d = 0.0; e = 0.0;
102     for(i=0; i<nsf; i++) {
103       d += target[i]*resp[i];
104       e += resp[i]*resp[i];
105     }
106     g = d/(e+.1);
107     score = g*d;
108     /*printf ("score: %f %f %f %f\n", target[0],d,e,score);*/
109     if (score >= bscore) {
110       bscore = score;
111       *gain = g;
112       *index = k;
113     }
114     
115     /* Synthesise next entry */
116     
117     if (k) {
118       for(i=nsf-1; i>=1; i--)
119         resp[i] = resp[i-1] + codebook[k-1]*h[i];
120       resp[0] = codebook[k-1]*h[0];
121     }
122   }
123
124   free(resp);
125   free(h);
126   free(impulse);
127   return bscore;
128 }
129
130
131 void split_cb_search(
132 float target[],                 /* target vector */
133 float ak[],                     /* LPCs for this subframe */
134 float awk1[],                   /* Weighted LPCs for this subframe */
135 float awk2[],                   /* Weighted LPCs for this subframe */
136 void *par,                      /* Codebook/search parameters*/
137 int   p,                        /* number of LPC coeffs */
138 int   nsf,                      /* number of samples in subframe */
139 float *exc,
140 FrameBits *bits,
141 float *stack
142 )
143 {
144    int i,j;
145    float *resp, *E, *Ee;
146    float *t, *r, *e;
147    float *gains;
148    int *ind;
149    float *shape_cb, *gain_cb;
150    int shape_cb_size, gain_cb_size, subvect_size, nb_subvect;
151    split_cb_params *params;
152
153    params = (split_cb_params *) par;
154    subvect_size = params->subvect_size;
155    nb_subvect = params->nb_subvect;
156    shape_cb_size = 1<<params->shape_bits;
157    shape_cb = params->shape_cb;
158    gain_cb_size = 1<<params->gain_bits;
159    gain_cb = params->gain_cb;
160    resp = PUSH(stack, shape_cb_size*subvect_size);
161    E = PUSH(stack, shape_cb_size);
162    Ee = PUSH(stack, shape_cb_size);
163    t = PUSH(stack, nsf);
164    r = PUSH(stack, nsf);
165    e = PUSH(stack, nsf);
166    gains = PUSH(stack, nb_subvect);
167    ind = (int*)PUSH(stack, nb_subvect);
168    
169
170    for (i=0;i<nsf;i++)
171       t[i]=target[i];
172    for (i=0;i<shape_cb_size;i++)
173    {
174       float *res = resp+i*subvect_size;
175       residue_zero(shape_cb+i*subvect_size, awk1, res, subvect_size, p);
176       syn_filt_zero(res, ak, res, subvect_size, p);
177       syn_filt_zero(res, awk2, res, subvect_size,p);
178       E[i]=0;
179       for(j=0;j<subvect_size;j++)
180          E[i]+=res[j]*res[j];
181       Ee[i]=0;
182       for(j=0;j<subvect_size;j++)
183          Ee[i]+=shape_cb[i*subvect_size+j]*shape_cb[i*subvect_size+j];
184       
185    }
186    for (i=0;i<nb_subvect;i++)
187    {
188       int best_index=0;
189       float g, corr, best_gain=0, score, best_score=-1;
190       for (j=0;j<shape_cb_size;j++)
191       {
192          corr=xcorr(resp+j*subvect_size,t+subvect_size*i,subvect_size);
193          score=corr*corr/(.001+E[j]);
194          g = corr/(.001+E[j]);
195          if (score>best_score)
196          {
197             best_index=j;
198             best_score=score;
199             best_gain=corr/(.001+E[j]);
200          }
201       }
202       frame_bits_pack(bits,best_index,params->shape_bits);
203       if (best_gain>0)
204          frame_bits_pack(bits,0,1);
205       else
206           frame_bits_pack(bits,1,1);        
207       ind[i]=best_index;
208       gains[i]=best_gain*Ee[ind[i]];
209
210       for (j=0;j<nsf;j++)
211          e[j]=0;
212       for (j=0;j<subvect_size;j++)
213          e[subvect_size*i+j]=best_gain*shape_cb[best_index*subvect_size+j];
214       residue_zero(e, awk1, r, nsf, p);
215       syn_filt_zero(r, ak, r, nsf, p);
216       syn_filt_zero(r, awk2, r, nsf,p);
217       for (j=0;j<nsf;j++)
218          t[j]-=r[j];
219    }
220
221    {
222       int best_vq_index=0, max_index;
223       float max_gain=0, log_max, min_dist=0, *sign;
224
225       if (gain_cb) /*If no gain codebook, do not quantize (for testing/debugging) */
226       {
227          sign = PUSH(stack, nb_subvect);
228          for (i=0;i<nb_subvect;i++)
229          {
230             if (gains[i]<0)
231             {
232                gains[i]=-gains[i];
233                sign[i]=-1;
234             } else {
235                sign[i]=1;
236             }
237          }
238          for (i=0;i<nb_subvect;i++)
239             if (gains[i]>max_gain)
240                max_gain=gains[i];
241          log_max=log(max_gain+1);
242          max_index = (int)(floor(.5+log_max-3));
243          if (max_index>7)
244             max_index=7;
245          if (max_index<0)
246             max_index=0;
247          max_gain=1/exp(max_index+3.0);
248          for (i=0;i<nb_subvect;i++)
249             gains[i]*=max_gain;
250          frame_bits_pack(bits,max_index,3);
251
252          /*Vector quantize gains[i]*/
253          if (nb_subvect<=5)
254          {
255          best_vq_index = vq_index(gains, gain_cb, nb_subvect, gain_cb_size);
256          frame_bits_pack(bits,best_vq_index,params->gain_bits);
257          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
258          for (i=0;i<nb_subvect;i++)
259             gains[i]= sign[i]*gain_cb[best_vq_index*nb_subvect+i]/max_gain/(Ee[ind[i]]+.001);
260          } else
261          {
262             float tmp[5];
263             int best_vq_index2;
264          best_vq_index = vq_index(gains, gain_cb, nb_subvect/2, gain_cb_size);
265          for (i=0;i<5;i++)
266             tmp[i]=gains[i]-gain_cb[best_vq_index*nb_subvect/2+i];
267          best_vq_index2 = vq_index(tmp, exc_gains_wb2_table, nb_subvect/2, 256);
268
269          frame_bits_pack(bits,best_vq_index,params->gain_bits);
270          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
271          for (i=0;i<nb_subvect/2;i++)
272             gains[i]= sign[i]*(gain_cb[best_vq_index*nb_subvect/2+i]+exc_gains_wb2_table[best_vq_index2*nb_subvect/2+i])/max_gain/(Ee[ind[i]]+.001);
273
274
275          best_vq_index = vq_index(gains+5, gain_cb, nb_subvect/2, gain_cb_size);
276          frame_bits_pack(bits,best_vq_index,params->gain_bits);
277          for (i=0;i<5;i++)
278             tmp[i]=gains[i+5]-gain_cb[best_vq_index*nb_subvect/2+i];
279          best_vq_index2 = vq_index(tmp, exc_gains_wb2_table, nb_subvect/2, 256);
280
281          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
282          for (i=0;i<nb_subvect/2;i++)
283             gains[i+5]= sign[i+5]*(gain_cb[best_vq_index*nb_subvect/2+i]+exc_gains_wb2_table[best_vq_index2*nb_subvect/2+i])/max_gain/(Ee[ind[i+5]]+.001);
284          }
285
286     
287
288          POP(stack);
289       } else {
290          printf ("exc: ");
291          for (i=0;i<nb_subvect;i++)
292             printf ("%f ", gains[i]);
293          printf ("\n");
294          for (i=0;i<nb_subvect;i++)
295             gains[i]= gains[i]/(Ee[ind[i]]+.001);
296       }
297
298       for (i=0;i<nb_subvect;i++)
299          for (j=0;j<subvect_size;j++)
300             exc[subvect_size*i+j]+=gains[i]*shape_cb[ind[i]*subvect_size+j];
301
302    }
303
304    /*TODO: Perform joint optimization of gains*/
305    
306    for (i=0;i<nsf;i++)
307       target[i]=t[i];
308
309    POP(stack);
310    POP(stack);
311    POP(stack);
312    POP(stack);
313    POP(stack);
314    POP(stack);
315    POP(stack);
316    POP(stack);
317 }
318
319
320
321
322 void split_cb_search_wb(
323 float target[],                 /* target vector */
324 float ak[],                     /* LPCs for this subframe */
325 float awk1[],                   /* Weighted LPCs for this subframe */
326 float awk2[],                   /* Weighted LPCs for this subframe */
327 void *par,                      /* Codebook/search parameters*/
328 int   p,                        /* number of LPC coeffs */
329 int   nsf,                      /* number of samples in subframe */
330 float *exc,
331 FrameBits *bits,
332 float *stack
333 )
334 {
335    int i,j;
336    float *resp, *E, *Ee;
337    float *t, *r, *e, *tresp;
338    float *gains;
339    int *ind;
340    float *shape_cb, *gain_cb;
341    int shape_cb_size, gain_cb_size, subvect_size, nb_subvect;
342    split_cb_params *params;
343
344    params = (split_cb_params *) par;
345    subvect_size = params->subvect_size;
346    nb_subvect = params->nb_subvect;
347    shape_cb_size = 1<<params->shape_bits;
348    shape_cb = params->shape_cb;
349    gain_cb_size = 1<<params->gain_bits;
350    gain_cb = params->gain_cb;
351    resp = PUSH(stack, shape_cb_size*subvect_size);
352    tresp = PUSH(stack, shape_cb_size*nsf);
353    E = PUSH(stack, shape_cb_size);
354    Ee = PUSH(stack, shape_cb_size);
355    t = PUSH(stack, nsf);
356    r = PUSH(stack, nsf);
357    e = PUSH(stack, nsf);
358    gains = PUSH(stack, nb_subvect);
359    ind = (int*)PUSH(stack, nb_subvect);
360    
361
362    for (i=0;i<nsf;i++)
363       t[i]=target[i];
364    for (i=0;i<shape_cb_size;i++)
365    {
366       float *res = resp+i*subvect_size;
367       residue_zero(shape_cb+i*subvect_size, awk1, res, subvect_size, p);
368       syn_filt_zero(res, ak, res, subvect_size, p);
369       syn_filt_zero(res, awk2, res, subvect_size,p);
370       E[i]=0;
371       for(j=0;j<subvect_size;j++)
372          E[i]+=res[j]*res[j];
373       Ee[i]=0;
374       for(j=0;j<subvect_size;j++)
375          Ee[i]+=shape_cb[i*subvect_size+j]*shape_cb[i*subvect_size+j];
376       
377    }
378 #if 0
379    {
380       int half,k;
381    for (half=0;half<2;half++)
382    {
383       int nb_half=nb_subvect/2;
384       int max_subvect=0;
385       float max_energy=0;
386       syn_filt_zero(t+half*nsf/2, awk1, r, nsf/2, p);
387       residue_zero(r, ak, r, nsf/2, p);
388       residue_zero(r, awk2, r, nsf/2,p);
389       for (i=0;i<nb_half;i++)
390       {
391          float energy=0;
392          for (k=0;k<subvect_size;k++)
393             energy+=r[subvect_size*i+k]*r[subvect_size*i+k];
394          if (energy>max_energy)
395          {
396             max_subvect=i;
397             max_energy=energy;
398          }
399       }
400       printf ("max_energy: %d %f\n", max_subvect, max_energy);
401       
402       for (i=0;i<nb_half;i++)
403       {
404          int nb_times=1;
405          if (i==max_subvect)
406             nb_times++;
407
408          for (k=0;k<nb_times;k++)
409          {
410          int best_index=0;
411          float g, corr, best_gain=0, score, best_score=-1;
412          for (j=0;j<shape_cb_size;j++)
413          {
414             corr=xcorr(resp+j*subvect_size,t+subvect_size*(i+half*nb_half),subvect_size);
415             score=corr*corr/(.001+E[j]);
416             g = corr/(.001+E[j]);
417             if (score>best_score)
418             {
419                best_index=j;
420                best_score=score;
421                best_gain=corr/(.001+E[j]);
422             }
423          }
424          for (j=0;j<nsf;j++)
425             e[j]=0;
426          for (j=0;j<subvect_size;j++)
427             e[subvect_size*(i+half*nb_half)+j]=best_gain*shape_cb[best_index*subvect_size+j];
428          residue_zero(e, awk1, r, nsf, p);
429          syn_filt_zero(r, ak, r, nsf, p);
430          syn_filt_zero(r, awk2, r, nsf,p);
431          for (j=0;j<nsf;j++)
432             t[j]-=r[j];
433          for (j=0;j<nsf;j++)
434             exc[j]+=e[j];
435          }
436       }
437    }
438    }
439 #else
440    for (i=0;i<nb_subvect;i++)
441    {
442       int best_index=0;
443       float g, corr, best_gain=0, score, best_score=-1;
444       for (j=0;j<shape_cb_size;j++)
445       {
446          corr=xcorr(resp+j*subvect_size,t+subvect_size*i,subvect_size);
447          score=corr*corr/(.001+E[j]);
448          g = corr/(.001+E[j]);
449          if (score>best_score)
450          {
451             best_index=j;
452             best_score=score;
453             best_gain=corr/(.001+E[j]);
454          }
455       }
456       frame_bits_pack(bits,best_index,params->shape_bits);
457       if (best_gain>0)
458          frame_bits_pack(bits,0,1);
459       else
460           frame_bits_pack(bits,1,1);        
461       ind[i]=best_index;
462       gains[i]=best_gain;
463
464       for (j=0;j<nsf;j++)
465          e[j]=0;
466       for (j=0;j<subvect_size;j++)
467          e[subvect_size*i+j]=best_gain*shape_cb[best_index*subvect_size+j];
468       residue_zero(e, awk1, r, nsf, p);
469       syn_filt_zero(r, ak, r, nsf, p);
470       syn_filt_zero(r, awk2, r, nsf,p);
471       for (j=0;j<nsf;j++)
472          tresp[i*nsf+j]=r[j];
473       for (j=0;j<nsf;j++)
474          t[j]-=r[j];
475       /*for (j=0;j<nsf;j++)
476         exc[j]+=e[j];*/
477    }
478    {
479       float A[10][10];
480       float b[10];
481       float c[10];
482       for (i=0;i<10;i++)
483          for (j=0;j<10;j++)
484             A[i][j]=xcorr(tresp+i*nsf, tresp+j*nsf, nsf);
485       for (i=0;i<10;i++)
486          b[i]=xcorr(target,tresp+i*nsf,nsf);
487       for (i=0;i<10;i++)
488          A[i][i]+=.01;
489
490
491       solve(&A[0][0],b,c, 10);
492       for (i=0;i<10;i++)
493         gains[i]*=c[i];
494       
495       for (i=0;i<10;i++)
496          gains[i]*=Ee[ind[i]];
497
498
499
500    {
501       int best_vq_index=0, max_index;
502       float max_gain=0, log_max, min_dist=0, *sign;
503
504       if (gain_cb) /*If no gain codebook, do not quantize (for testing/debugging) */
505       {
506          sign = PUSH(stack, nb_subvect);
507          for (i=0;i<nb_subvect;i++)
508          {
509             if (gains[i]<0)
510             {
511                gains[i]=-gains[i];
512                sign[i]=-1;
513             } else {
514                sign[i]=1;
515             }
516          }
517          for (i=0;i<nb_subvect;i++)
518             if (gains[i]>max_gain)
519                max_gain=gains[i];
520          log_max=log(max_gain+1);
521          max_index = (int)(floor(.5+log_max-3));
522          if (max_index>7)
523             max_index=7;
524          if (max_index<0)
525             max_index=0;
526          max_gain=1/exp(max_index+3.0);
527          for (i=0;i<nb_subvect;i++)
528             gains[i]*=max_gain;
529          frame_bits_pack(bits,max_index,3);
530
531          /*Vector quantize gains[i]*/
532          if (nb_subvect<=5)
533          {
534          best_vq_index = vq_index(gains, gain_cb, nb_subvect, gain_cb_size);
535          frame_bits_pack(bits,best_vq_index,params->gain_bits);
536          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
537          for (i=0;i<nb_subvect;i++)
538             gains[i]= sign[i]*gain_cb[best_vq_index*nb_subvect+i]/max_gain/(Ee[ind[i]]+.001);
539          } else
540          {
541             float tmp[5];
542             int best_vq_index2;
543          best_vq_index = vq_index(gains, gain_cb, nb_subvect/2, gain_cb_size);
544          for (i=0;i<5;i++)
545             tmp[i]=gains[i]-gain_cb[best_vq_index*nb_subvect/2+i];
546          best_vq_index2 = vq_index(tmp, exc_gains_wb2_table, nb_subvect/2, 256);
547
548          frame_bits_pack(bits,best_vq_index,params->gain_bits);
549          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
550          for (i=0;i<nb_subvect/2;i++)
551             gains[i]= sign[i]*(gain_cb[best_vq_index*nb_subvect/2+i]+exc_gains_wb2_table[best_vq_index2*nb_subvect/2+i])/max_gain/(Ee[ind[i]]+.001);
552
553
554          best_vq_index = vq_index(gains+5, gain_cb, nb_subvect/2, gain_cb_size);
555          frame_bits_pack(bits,best_vq_index,params->gain_bits);
556          for (i=0;i<5;i++)
557             tmp[i]=gains[i+5]-gain_cb[best_vq_index*nb_subvect/2+i];
558          best_vq_index2 = vq_index(tmp, exc_gains_wb2_table, nb_subvect/2, 256);
559
560          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
561          for (i=0;i<nb_subvect/2;i++)
562             gains[i+5]= sign[i+5]*(gain_cb[best_vq_index*nb_subvect/2+i]+exc_gains_wb2_table[best_vq_index2*nb_subvect/2+i])/max_gain/(Ee[ind[i+5]]+.001);
563          }
564          
565          
566       } else {
567          
568          for (i=0;i<10;i++)
569             gains[i]/=Ee[ind[i]]+.001;
570          
571       }
572    }
573       for (i=0;i<10;i++)
574          for (j=0;j<subvect_size;j++)
575             e[subvect_size*i+j]=gains[i]*shape_cb[ind[i]*subvect_size+j];
576
577       for (j=0;j<nsf;j++)
578          exc[j]+=e[j];
579       residue_zero(e, awk1, r, nsf, p);
580       syn_filt_zero(r, ak, r, nsf, p);
581       syn_filt_zero(r, awk2, r, nsf,p);
582       for (j=0;j<nsf;j++)
583          target[j]-=r[j];
584
585    }
586 #endif
587
588    /*TODO: Perform joint optimization of gains*/
589    
590    /*for (i=0;i<nsf;i++)
591       target[i]=t[i];
592    */
593    POP(stack);
594    POP(stack);
595    POP(stack);
596    POP(stack);
597    POP(stack);
598    POP(stack);
599    POP(stack);
600    POP(stack);
601 }
602
603
604
605
606 void split_cb_unquant(
607 float *exc,
608 void *par,                      /* non-overlapping codebook */
609 int   nsf,                      /* number of samples in subframe */
610 FrameBits *bits,
611 float *stack
612 )
613 {
614    int i,j;
615    int *ind;
616    float *gains;
617    float *sign;
618    int max_gain_ind, vq_gain_ind;
619    float max_gain, *Ee;
620    float *shape_cb, *gain_cb;
621    int shape_cb_size, gain_cb_size, subvect_size, nb_subvect;
622    split_cb_params *params;
623
624    params = (split_cb_params *) par;
625    subvect_size = params->subvect_size;
626    nb_subvect = params->nb_subvect;
627    shape_cb_size = 1<<params->shape_bits;
628    shape_cb = params->shape_cb;
629    gain_cb_size = 1<<params->gain_bits;
630    gain_cb = params->gain_cb;
631    
632    ind = (int*)PUSH(stack, nb_subvect);
633    gains = PUSH(stack, nb_subvect);
634    sign = PUSH(stack, nb_subvect);
635    Ee=PUSH(stack, nb_subvect);
636
637    for (i=0;i<nb_subvect;i++)
638    {
639       ind[i] = frame_bits_unpack_unsigned(bits, params->shape_bits);
640       if (frame_bits_unpack_unsigned(bits, 1))
641          sign[i]=-1;
642       else
643          sign[i]=1;
644       Ee[i]=.001;
645       for (j=0;j<subvect_size;j++)
646          Ee[i]+=shape_cb[ind[i]*subvect_size+j]*shape_cb[ind[i]*subvect_size+j];
647    }
648    max_gain_ind = frame_bits_unpack_unsigned(bits, 3);
649    vq_gain_ind = frame_bits_unpack_unsigned(bits, params->gain_bits);
650    printf ("unquant gains ind: %d %d\n", max_gain_ind, vq_gain_ind);
651
652    max_gain=exp(max_gain_ind+3.0);
653    for (i=0;i<nb_subvect;i++)
654       gains[i] = sign[i]*gain_cb[vq_gain_ind*nb_subvect+i]*max_gain/Ee[i];
655    
656    printf ("unquant gains: ");
657    for (i=0;i<nb_subvect;i++)
658       printf ("%f ", gains[i]);
659    printf ("\n");
660
661    for (i=0;i<nb_subvect;i++)
662       for (j=0;j<subvect_size;j++)
663          exc[subvect_size*i+j]+=gains[i]*shape_cb[ind[i]*subvect_size+j];
664    
665    POP(stack);
666    POP(stack);
667    POP(stack);
668    POP(stack);
669 }