printf ("gain correction %f\n", f);
void mpulse_search(float target[],			/* target vector */float ak[],			/* LPCs for this subframe */float awk1[],			/* Weighted LPCs for this subframe */float awk2[],			/* Weighted LPCs for this subframe */void *par,                      /* Codebook/search parameters*/int   p,                        /* number of LPC coeffs */int   nsf,                      /* number of samples in subframe */float *exc,FrameBits *bits,float *stack){   int i,j, nb_pulse;   float *resp, *resp2, *energy, *t, *e, *pulses;   float te=0,ee=0;   float g, gain_coef;   int nb_tracks, track_ind_bits;   int *tracks, *signs, *tr, *nb;   mpulse_params *params;   int pulses_per_track;   params = (mpulse_params *) par;   nb_pulse=params->nb_pulse;   nb_tracks=params->nb_tracks;   pulses_per_track=nb_pulse/nb_tracks;   track_ind_bits=params->track_ind_bits;   gain_coef=params->gain_coef;   tracks = (int*)PUSH(stack,nb_pulse);   signs = (int*)PUSH(stack,nb_pulse);   tr = (int*)PUSH(stack,pulses_per_track);   nb = (int*)PUSH(stack,nb_tracks);   resp=PUSH(stack, nsf);   resp2=PUSH(stack, nsf);   energy=PUSH(stack, nsf);   t=PUSH(stack, nsf);   e=PUSH(stack, nsf);   pulses=PUSH(stack, nsf);      /*Compute optimal (real) excitation from target*/   syn_filt_zero(target, awk1, e, nsf, p);   residue_zero(e, ak, e, nsf, p);   residue_zero(e, awk2, e, nsf, p);   for (i=0;i<nsf;i++)   {      pulses[i]=0;      te+=target[i]*target[i];      ee+=e[i]*e[i];   }   /*Compute global gain (coef found from linear regression and tweaking)*/   g=gain_coef/sqrt(nb_pulse)*exp(0.18163*log(te+1)+0.17293*log(ee+1));      e[0]=1;   for (i=1;i<nsf;i++)      e[i]=0;   /*Impulse response of W(z)/A(z)*/   residue_zero(e, awk1, resp, nsf, p);   syn_filt_zero(resp, ak, resp, nsf, p);   syn_filt_zero(resp, awk2, resp, nsf, p);      /*Impulse response * gain*/   for (i=0;i<nsf;i++)      resp2[i]=g*resp[i];   for (i=0;i<nsf;i++)      e[i]=0;   for (i=0;i<nsf;i++)      t[i]=target[i];   for (i=0;i<nb_tracks;i++)      nb[i]=0;   /*For all pulses*/   for (i=0;i<nb_pulse;i++)   {      float best_score=1e30, best_gain=0;      int best_ind=0;      /*For all positions*/      energy[0]=0;      for (j=1;j<nsf;j++)         energy[j]=energy[j-1]+t[j-1]*t[j-1];      /*For each position*/      for (j=0;j<nsf;j++)      {         int k;         float dist;         float *base=t+j;         /*Fill any track until it's full*/         if (nb[j%nb_tracks]==pulses_per_track)              continue;         /*Constrain search in alternating tracks*/         /*FIXME: Should get rid of this, it's *really* slow*/         /*if ((i%nb_tracks) != (j%nb_tracks))           continue;*/         /*Try for positive sign*/         if (pulses[j]>=0)         {            dist=energy[j];            for (k=0;k<nsf-j;k++)            {               float tmp=(base[k]-resp2[k]);               dist+=tmp*tmp;            }            if (dist<best_score || j==0)            {               best_score=dist;               best_gain=g;               best_ind=j;            }         }         /*Try again for negative sign*/         if (pulses[j]<=0)         {            dist=energy[j];            for (k=0;k<nsf-j;k++)            {               float tmp=(base[k]+resp2[k]);               dist+=tmp*tmp;            }            if (dist<best_score)            {               best_score=dist;               best_gain=-g;               best_ind=j;            }         }      }#ifdef DEBUG      printf ("best pulse: %d %d %f %f %f %f\n", i, best_ind, best_gain, te, ee, g);#endif      /*Remove pulse contribution from target*/      for (j=best_ind;j<nsf;j++)         t[j] -= best_gain * resp[j-best_ind];      e[best_ind]+=best_gain;      if (best_gain>0)         pulses[best_ind]+=1;      else         pulses[best_ind]-=1;      {         int t=best_ind%nb_tracks;         tracks[t*pulses_per_track+nb[t]] = best_ind/nb_tracks;         signs[t*pulses_per_track+nb[t]]  = best_gain >= 0 ? 0 : 1;         nb[t]++;      }   }   /*Global gain re-estimation*/   if (1) {      float f;      int quant_gain;      residue_zero(e, awk1, resp, nsf, p);      syn_filt_zero(resp, ak, resp, nsf, p);      syn_filt_zero(resp, awk2, resp, nsf, p);      f=((.1+(xcorr(resp,target,nsf)))/(.1+xcorr(resp,resp,nsf)));      printf ("gain correction %f\n", f);      /*for (i=0;i<nsf;i++)        e[i]*=f;*/      g *= f;      if (g<0)         g=0;            quant_gain=(int)floor(.5+8*(log(1+fabs(g))-1));      if (quant_gain<0)         quant_gain=0;      if (quant_gain>127)         quant_gain=127;      frame_bits_pack(bits,quant_gain,7);      g=exp((quant_gain/8.0)+1);            for (i=0;i<nsf;i++)         e[i]=g*pulses[i];#ifdef DEBUG      printf ("global gain = %f\n", g);#endif      for (i=0;i<nsf;i++)         t[i]=target[i]-f*resp[i];   }#ifdef DEBUG   for (i=0;i<nsf;i++)      printf ("%f ", e[i]);   printf ("\n");#endif   for (i=0;i<nsf;i++)      exc[i]+=e[i];   for (i=0;i<nsf;i++)      target[i]=t[i];      for (i=0;i<nb_tracks;i++)   {      int bit1, ind=0;      bit1=porder(tracks+i*pulses_per_track, signs+i*pulses_per_track,tr,pulses_per_track);      frame_bits_pack(bits,bit1,1);      for (j=0;j<pulses_per_track;j++)      {         ind*=nsf/nb_tracks;         ind+=tr[j];         /*printf ("%d ", ind);*/      }            frame_bits_pack(bits,ind,track_ind_bits);      /*printf ("track %d %d:", i, ind);      for (j=0;j<pulses_per_track;j++)        printf ("%d ", tr[j]);        printf ("\n");*/   }   POP(stack);   POP(stack);   POP(stack);   POP(stack);   POP(stack);   POP(stack);   POP(stack);   POP(stack);   POP(stack);   POP(stack);}