Added joint optimization of excitation gains
[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,k,half;
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    for (half=0;half<2;half++)
380    {
381       int nb_half=nb_subvect/2;
382       int max_subvect=0;
383       float max_energy=0;
384       syn_filt_zero(t+half*nsf/2, awk1, r, nsf/2, p);
385       residue_zero(r, ak, r, nsf/2, p);
386       residue_zero(r, awk2, r, nsf/2,p);
387       for (i=0;i<nb_half;i++)
388       {
389          float energy=0;
390          for (k=0;k<subvect_size;k++)
391             energy+=r[subvect_size*i+k]*r[subvect_size*i+k];
392          if (energy>max_energy)
393          {
394             max_subvect=i;
395             max_energy=energy;
396          }
397       }
398       printf ("max_energy: %d %f\n", max_subvect, max_energy);
399       
400       for (i=0;i<nb_half;i++)
401       {
402          int nb_times=1;
403          if (i==max_subvect)
404             nb_times++;
405
406          for (k=0;k<nb_times;k++)
407          {
408          int best_index=0;
409          float g, corr, best_gain=0, score, best_score=-1;
410          for (j=0;j<shape_cb_size;j++)
411          {
412             corr=xcorr(resp+j*subvect_size,t+subvect_size*(i+half*nb_half),subvect_size);
413             score=corr*corr/(.001+E[j]);
414             g = corr/(.001+E[j]);
415             if (score>best_score)
416             {
417                best_index=j;
418                best_score=score;
419                best_gain=corr/(.001+E[j]);
420             }
421          }
422          for (j=0;j<nsf;j++)
423             e[j]=0;
424          for (j=0;j<subvect_size;j++)
425             e[subvect_size*(i+half*nb_half)+j]=best_gain*shape_cb[best_index*subvect_size+j];
426          residue_zero(e, awk1, r, nsf, p);
427          syn_filt_zero(r, ak, r, nsf, p);
428          syn_filt_zero(r, awk2, r, nsf,p);
429          for (j=0;j<nsf;j++)
430             t[j]-=r[j];
431          for (j=0;j<nsf;j++)
432             exc[j]+=e[j];
433          }
434       }
435    }
436 #else
437    for (i=0;i<nb_subvect;i++)
438    {
439       int best_index=0;
440       float g, corr, best_gain=0, score, best_score=-1;
441       for (j=0;j<shape_cb_size;j++)
442       {
443          corr=xcorr(resp+j*subvect_size,t+subvect_size*i,subvect_size);
444          score=corr*corr/(.001+E[j]);
445          g = corr/(.001+E[j]);
446          if (score>best_score)
447          {
448             best_index=j;
449             best_score=score;
450             best_gain=corr/(.001+E[j]);
451          }
452       }
453       frame_bits_pack(bits,best_index,params->shape_bits);
454       if (best_gain>0)
455          frame_bits_pack(bits,0,1);
456       else
457           frame_bits_pack(bits,1,1);        
458       ind[i]=best_index;
459       gains[i]=best_gain;
460
461       for (j=0;j<nsf;j++)
462          e[j]=0;
463       for (j=0;j<subvect_size;j++)
464          e[subvect_size*i+j]=best_gain*shape_cb[best_index*subvect_size+j];
465       residue_zero(e, awk1, r, nsf, p);
466       syn_filt_zero(r, ak, r, nsf, p);
467       syn_filt_zero(r, awk2, r, nsf,p);
468       for (j=0;j<nsf;j++)
469          tresp[i*nsf+j]=r[j];
470       for (j=0;j<nsf;j++)
471          t[j]-=r[j];
472       /*for (j=0;j<nsf;j++)
473         exc[j]+=e[j];*/
474    }
475    {
476       float A[10][10];
477       float b[10];
478       float c[10];
479       for (i=0;i<10;i++)
480          for (j=0;j<10;j++)
481             A[i][j]=xcorr(tresp+i*nsf, tresp+j*nsf, nsf);
482       for (i=0;i<10;i++)
483          b[i]=xcorr(target,tresp+i*nsf,nsf);
484       for (i=0;i<10;i++)
485          A[i][i]+=.01;
486
487
488       solve(&A[0][0],b,c, 10);
489       for (i=0;i<10;i++)
490          gains[i]*=c[i];
491       
492       for (i=0;i<10;i++)
493          gains[i]*=Ee[ind[i]];
494
495
496
497    {
498       int best_vq_index=0, max_index;
499       float max_gain=0, log_max, min_dist=0, *sign;
500
501       if (gain_cb) /*If no gain codebook, do not quantize (for testing/debugging) */
502       {
503          sign = PUSH(stack, nb_subvect);
504          for (i=0;i<nb_subvect;i++)
505          {
506             if (gains[i]<0)
507             {
508                gains[i]=-gains[i];
509                sign[i]=-1;
510             } else {
511                sign[i]=1;
512             }
513          }
514          for (i=0;i<nb_subvect;i++)
515             if (gains[i]>max_gain)
516                max_gain=gains[i];
517          log_max=log(max_gain+1);
518          max_index = (int)(floor(.5+log_max-3));
519          if (max_index>7)
520             max_index=7;
521          if (max_index<0)
522             max_index=0;
523          max_gain=1/exp(max_index+3.0);
524          for (i=0;i<nb_subvect;i++)
525             gains[i]*=max_gain;
526          frame_bits_pack(bits,max_index,3);
527
528          /*Vector quantize gains[i]*/
529          if (nb_subvect<=5)
530          {
531          best_vq_index = vq_index(gains, gain_cb, nb_subvect, gain_cb_size);
532          frame_bits_pack(bits,best_vq_index,params->gain_bits);
533          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
534          for (i=0;i<nb_subvect;i++)
535             gains[i]= sign[i]*gain_cb[best_vq_index*nb_subvect+i]/max_gain/(Ee[ind[i]]+.001);
536          } else
537          {
538             float tmp[5];
539             int best_vq_index2;
540          best_vq_index = vq_index(gains, gain_cb, nb_subvect/2, gain_cb_size);
541          for (i=0;i<5;i++)
542             tmp[i]=gains[i]-gain_cb[best_vq_index*nb_subvect/2+i];
543          best_vq_index2 = vq_index(tmp, exc_gains_wb2_table, nb_subvect/2, 256);
544
545          frame_bits_pack(bits,best_vq_index,params->gain_bits);
546          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
547          for (i=0;i<nb_subvect/2;i++)
548             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);
549
550
551          best_vq_index = vq_index(gains+5, gain_cb, nb_subvect/2, gain_cb_size);
552          frame_bits_pack(bits,best_vq_index,params->gain_bits);
553          for (i=0;i<5;i++)
554             tmp[i]=gains[i+5]-gain_cb[best_vq_index*nb_subvect/2+i];
555          best_vq_index2 = vq_index(tmp, exc_gains_wb2_table, nb_subvect/2, 256);
556
557          printf ("best_gains_vq_index %d %f %d\n", best_vq_index, min_dist, max_index);
558          for (i=0;i<nb_subvect/2;i++)
559             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);
560          }
561          
562          
563       } else {
564          
565          for (i=0;i<10;i++)
566             gains[i]/=Ee[ind[i]]+.001;
567          
568       }
569    }
570       for (i=0;i<10;i++)
571          for (j=0;j<subvect_size;j++)
572             e[subvect_size*i+j]=gains[i]*shape_cb[ind[i]*subvect_size+j];
573
574       for (j=0;j<nsf;j++)
575          exc[j]+=e[j];
576       residue_zero(e, awk1, r, nsf, p);
577       syn_filt_zero(r, ak, r, nsf, p);
578       syn_filt_zero(r, awk2, r, nsf,p);
579       for (j=0;j<nsf;j++)
580          target[j]-=r[j];
581
582    }
583 #endif
584
585    /*TODO: Perform joint optimization of gains*/
586    
587    /*for (i=0;i<nsf;i++)
588       target[i]=t[i];
589    */
590    POP(stack);
591    POP(stack);
592    POP(stack);
593    POP(stack);
594    POP(stack);
595    POP(stack);
596    POP(stack);
597    POP(stack);
598 }
599
600
601
602
603 void split_cb_unquant(
604 float *exc,
605 void *par,                      /* non-overlapping codebook */
606 int   nsf,                      /* number of samples in subframe */
607 FrameBits *bits,
608 float *stack
609 )
610 {
611    int i,j;
612    int *ind;
613    float *gains;
614    float *sign;
615    int max_gain_ind, vq_gain_ind;
616    float max_gain, *Ee;
617    float *shape_cb, *gain_cb;
618    int shape_cb_size, gain_cb_size, subvect_size, nb_subvect;
619    split_cb_params *params;
620
621    params = (split_cb_params *) par;
622    subvect_size = params->subvect_size;
623    nb_subvect = params->nb_subvect;
624    shape_cb_size = 1<<params->shape_bits;
625    shape_cb = params->shape_cb;
626    gain_cb_size = 1<<params->gain_bits;
627    gain_cb = params->gain_cb;
628    
629    ind = (int*)PUSH(stack, nb_subvect);
630    gains = PUSH(stack, nb_subvect);
631    sign = PUSH(stack, nb_subvect);
632    Ee=PUSH(stack, nb_subvect);
633
634    for (i=0;i<nb_subvect;i++)
635    {
636       ind[i] = frame_bits_unpack_unsigned(bits, params->shape_bits);
637       if (frame_bits_unpack_unsigned(bits, 1))
638          sign[i]=-1;
639       else
640          sign[i]=1;
641       Ee[i]=.001;
642       for (j=0;j<subvect_size;j++)
643          Ee[i]+=shape_cb[ind[i]*subvect_size+j]*shape_cb[ind[i]*subvect_size+j];
644    }
645    max_gain_ind = frame_bits_unpack_unsigned(bits, 3);
646    vq_gain_ind = frame_bits_unpack_unsigned(bits, params->gain_bits);
647    printf ("unquant gains ind: %d %d\n", max_gain_ind, vq_gain_ind);
648
649    max_gain=exp(max_gain_ind+3.0);
650    for (i=0;i<nb_subvect;i++)
651       gains[i] = sign[i]*gain_cb[vq_gain_ind*nb_subvect+i]*max_gain/Ee[i];
652    
653    printf ("unquant gains: ");
654    for (i=0;i<nb_subvect;i++)
655       printf ("%f ", gains[i]);
656    printf ("\n");
657
658    for (i=0;i<nb_subvect;i++)
659       for (j=0;j<subvect_size;j++)
660          exc[subvect_size*i+j]+=gains[i]*shape_cb[ind[i]*subvect_size+j];
661    
662    POP(stack);
663    POP(stack);
664    POP(stack);
665    POP(stack);
666 }