fprintf (stderr, "roots!=st->lpcSize (found only %d roots)\n", roots); exit(1);
roots = lpc_to_lsp (st->lpc, st->lpcSize, st->lsp, 11, 0.02, st->stack); if (roots!=st->lpcSize) {    for (i=0;i<st->lpcSize;i++) { st->lsp[i]=cos(M_PI*((float)(i+1))/(st->lpcSize+1)); } }
void sb_encode(void *state, float *in, SpeexBits *bits){   SBEncState *st;   int i, roots, sub;   st = state;   /* Compute the two sub-bands by filtering with h0 and h1*/   fir_mem(in, h0, st->x0, st->full_frame_size, QMF_ORDER, st->h0_mem);   fir_mem(in, h1, st->x1, st->full_frame_size, QMF_ORDER, st->h1_mem);   /* Down-sample x0 and x1 */   for (i=0;i<st->frame_size;i++)      st->x1d[i]=st->x1[i<<1];   for (i=0;i<st->frame_size;i++)      st->x0d[i]=st->x0[i<<1];   /* Encode the narrowband part*/   nb_encode(st->st_low, st->x0d, bits);   speex_bits_pack(bits, 1, 1);   speex_bits_pack(bits, st->submodeID, SB_SUBMODE_BITS);   /* High-band buffering / sync with low band */   for (i=0;i<st->windowSize-st->frame_size;i++)      st->high[i] = st->high[st->frame_size+i];   for (i=0;i<st->frame_size;i++)      st->high[st->windowSize-st->frame_size+i]=st->x1d[i];   speex_move(st->excBuf, st->excBuf+st->frame_size, (st->bufSize-st->frame_size)*sizeof(float));   /* Start encoding the high-band */   for (i=0;i<st->windowSize;i++)      st->buf[i] = st->high[i] * st->window[i];   /* Compute auto-correlation */   autocorr(st->buf, st->autocorr, st->lpcSize+1, st->windowSize);   st->autocorr[0] += 1;        /* prevents NANs */   st->autocorr[0] *= st->lpc_floor; /* Noise floor in auto-correlation domain */   /* Lag windowing: equivalent to filtering in the power-spectrum domain */   for (i=0;i<st->lpcSize+1;i++)      st->autocorr[i] *= st->lagWindow[i];   /* Levinson-Durbin */   wld(st->lpc+1, st->autocorr, st->rc, st->lpcSize);   st->lpc[0]=1;   /* LPC to LSPs (x-domain) transform */   roots=lpc_to_lsp (st->lpc, st->lpcSize, st->lsp, 6, 0.01, st->stack);   if (roots!=st->lpcSize)   {      fprintf (stderr, "roots!=st->lpcSize (found only %d roots)\n", roots);      exit(1);   }   /* x-domain to angle domain*/   for (i=0;i<st->lpcSize;i++)      st->lsp[i] = acos(st->lsp[i]);   /* If null mode (no transmission), just set a couple things to zero*/   if (st->submodes[st->submodeID] == NULL)   {      for (i=0;i<st->frame_size;i++)         st->exc[i]=st->sw[i]=0;      for (i=0;i<st->lpcSize;i++)         st->mem_sw[i]=0;      st->first=1;      /* Final signal synthesis from excitation */      iir_mem2(st->exc, st->interp_qlpc, st->high, st->subframeSize, st->lpcSize, st->mem_sp);#ifndef RELEASE      /* Up-sample coded low-band and high-band*/      for (i=0;i<st->frame_size;i++)      {         st->x0[(i<<1)]=st->x0d[i];         st->x1[(i<<1)]=st->high[i];         st->x0[(i<<1)+1]=0;         st->x1[(i<<1)+1]=0;      }      /* Reconstruct the original */      fir_mem(st->x0, h0, st->y0, st->full_frame_size, QMF_ORDER, st->g0_mem);      fir_mem(st->x1, h1, st->y1, st->full_frame_size, QMF_ORDER, st->g1_mem);      for (i=0;i<st->full_frame_size;i++)         in[i]=2*(st->y0[i]-st->y1[i]);#endif      return;   }   /* LSP quantization */   SUBMODE(lsp_quant)(st->lsp, st->qlsp, st->lpcSize, bits);      /*printf ("high_lsp:");   for (i=0;i<st->lpcSize;i++)      printf (" %f", st->lsp[i]);      printf ("\n");*/   /*for (i=0;i<st->lpcSize;i++)     st->qlsp[i]=st->lsp[i];*/      if (st->first)   {      for (i=0;i<st->lpcSize;i++)         st->old_lsp[i] = st->lsp[i];      for (i=0;i<st->lpcSize;i++)         st->old_qlsp[i] = st->qlsp[i];   }      for (sub=0;sub<st->nbSubframes;sub++)   {      float *exc, *sp, *mem, *res, *target, *sw, tmp, filter_ratio;      int offset;      float rl, rh, eh=0, el=0;      int fold;      offset = st->subframeSize*sub;      sp=st->high+offset;      exc=st->excBuf+offset;      res=st->res+offset;      target=st->target+offset;      sw=st->sw+offset;      mem=PUSH(st->stack, st->lpcSize);            /* LSP interpolation (quantized and unquantized) */      tmp = (1.0 + sub)/st->nbSubframes;      for (i=0;i<st->lpcSize;i++)         st->interp_lsp[i] = (1-tmp)*st->old_lsp[i] + tmp*st->lsp[i];      for (i=0;i<st->lpcSize;i++)         st->interp_qlsp[i] = (1-tmp)*st->old_qlsp[i] + tmp*st->qlsp[i];            /* Compute interpolated LPCs (quantized and unquantized) */      for (i=0;i<st->lpcSize;i++)         st->interp_lsp[i] = cos(st->interp_lsp[i]);      for (i=0;i<st->lpcSize;i++)         st->interp_qlsp[i] = cos(st->interp_qlsp[i]);      lsp_enforce_margin(st->interp_lsp, st->lpcSize, .002);      lsp_enforce_margin(st->interp_qlsp, st->lpcSize, .002);      lsp_to_lpc(st->interp_lsp, st->interp_lpc, st->lpcSize,st->stack);      lsp_to_lpc(st->interp_qlsp, st->interp_qlpc, st->lpcSize, st->stack);      bw_lpc(st->gamma1, st->interp_lpc, st->bw_lpc1, st->lpcSize);      bw_lpc(st->gamma2, st->interp_lpc, st->bw_lpc2, st->lpcSize);      /* Compute mid-band (4000 Hz for wideband) response of low-band and high-band         filters */      rl=rh=0;      tmp=1;      for (i=0;i<=st->lpcSize;i++)      {         rh += tmp*st->interp_qlpc[i];         tmp = -tmp;      }      rl = ((EncState*)st->st_low)->pi_gain[sub];      rl=1/(fabs(rl)+.01);      rh=1/(fabs(rh)+.01);      /* Compute ratio, will help predict the gain */      filter_ratio=fabs(.01+rh)/(.01+fabs(rl));      fold = filter_ratio<5;      /*printf ("filter_ratio %f\n", filter_ratio);*/      fold=0;      /* Compute "real excitation" */      fir_mem2(sp, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, st->mem_sp2);      /* Compute energy of low-band and high-band excitation */      for (i=0;i<st->subframeSize;i++)         eh+=sqr(exc[i]);      if (!SUBMODE(innovation_quant)) {/* 1 for spectral folding excitation, 0 for stochastic */         float g;         /*speex_bits_pack(bits, 1, 1);*/         for (i=0;i<st->subframeSize;i++)            el+=sqr(((EncState*)st->st_low)->innov[offset+i]);         /* Gain to use if we want to use the low-band excitation for high-band */         g=eh/(.01+el);         g=sqrt(g);         g *= filter_ratio;         /* Gain quantization */         {            int quant = (int) floor(.5 + 27 + 8.0 * log((g+.0001)));            if (quant<0)               quant=0;            if (quant>31)               quant=31;            speex_bits_pack(bits, quant, 5);            g= .1*exp(quant/9.4);         }         /*printf ("folding gain: %f\n", g);*/         g /= filter_ratio;      } else {         float gc, scale, scale_1;         float *innov;         for (i=0;i<st->subframeSize;i++)            el+=sqr(((EncState*)st->st_low)->exc[offset+i]);         /*speex_bits_pack(bits, 0, 1);*/         innov = PUSH(st->stack, st->subframeSize);         gc = sqrt(1+eh)*filter_ratio/sqrt((1+el)*st->subframeSize);         {            int qgc = (int)floor(.5+3.7*(log(gc)+2));            if (qgc<0)               qgc=0;            if (qgc>15)               qgc=15;            speex_bits_pack(bits, qgc, 4);            gc = exp((1/3.7)*qgc-2);         }         scale = gc*sqrt(1+el)/filter_ratio;         scale_1 = 1/scale;         if (0 && rand()%5==0)         {            float sc = 1/sqrt(.1+eh/st->subframeSize);            if (rand()&1)               sc=-sc;            for (i=0;i<st->subframeSize;i++)            {               float tmp=exc[i]*sc;               if (i%8==0)                  printf ("\nhexc");               printf (" %f", tmp);            }         }         /* Reset excitation */         for (i=0;i<st->subframeSize;i++)            exc[i]=0;                  /* Compute zero response (ringing) of A(z/g1) / ( A(z/g2) * Aq(z) ) */         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sp[i];         iir_mem2(exc, st->interp_qlpc, exc, st->subframeSize, st->lpcSize, mem);         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sw[i];         filter_mem2(exc, st->bw_lpc1, st->bw_lpc2, res, st->subframeSize, st->lpcSize, mem);         /* Compute weighted signal */         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sw[i];         filter_mem2(sp, st->bw_lpc1, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, mem);#if 0         /*for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sp[i];         residue_mem(exc, st->bw_lpc1, res, st->subframeSize, st->lpcSize, mem);         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sw[i];            syn_filt_mem(res, st->bw_lpc2, res, st->subframeSize, st->lpcSize, mem);*/                  /* Compute weighted signal */         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sp[i];         residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, mem);         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sw[i];         syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, mem);#endif         /* Compute target signal */         for (i=0;i<st->subframeSize;i++)            target[i]=sw[i]-res[i];         for (i=0;i<st->subframeSize;i++)           exc[i]=0;         for (i=0;i<st->subframeSize;i++)            target[i]*=scale_1;                  /* Reset excitation */         for (i=0;i<st->subframeSize;i++)            innov[i]=0;         /*print_vec(target, st->subframeSize, "\ntarget");*/         SUBMODE(innovation_quant)(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,                                 SUBMODE(innovation_params), st->lpcSize, st->subframeSize,                                 innov, bits, st->stack, st->complexity);         /*print_vec(target, st->subframeSize, "after");*/         for (i=0;i<st->subframeSize;i++)            exc[i] += innov[i]*scale;         if (SUBMODE(double_codebook)) {            float *innov2 = PUSH(st->stack, st->subframeSize);            for (i=0;i<st->subframeSize;i++)               innov2[i]=0;            for (i=0;i<st->subframeSize;i++)               target[i]*=2.5;            SUBMODE(innovation_quant)(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2,                                       SUBMODE(innovation_params), st->lpcSize, st->subframeSize,                                       innov2, bits, st->stack, st->complexity);            for (i=0;i<st->subframeSize;i++)               innov2[i]*=scale*(1/2.5);            for (i=0;i<st->subframeSize;i++)               exc[i] += innov2[i];            POP(st->stack);         }         if (0) {            float en=0;            for (i=0;i<st->subframeSize;i++)               en+=exc[i]*exc[i];            en=sqrt(eh/(1+en));            for (i=0;i<st->subframeSize;i++)               exc[i]*=en;            printf ("correction high: %f\n", en);         }         POP(st->stack);      }#if 1         /*Keep the previous memory*/         for (i=0;i<st->lpcSize;i++)            mem[i]=st->mem_sp[i];         /* Final signal synthesis from excitation */         iir_mem2(exc, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, st->mem_sp);                  /* Compute weighted signal again, from synthesized speech (not sure it's the right thing) */         /*residue_mem(sp, st->bw_lpc1, sw, st->subframeSize, st->lpcSize, mem);         syn_filt_mem(sw, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, st->mem_sw);         */         filter_mem2(sp, st->bw_lpc1, st->bw_lpc2, sw, st->subframeSize, st->lpcSize, st->mem_sw);#endif                  POP(st->stack);   }#ifndef RELEASE   /* Up-sample coded low-band and high-band*/   for (i=0;i<st->frame_size;i++)   {      st->x0[(i<<1)]=st->x0d[i];      st->x1[(i<<1)]=st->high[i];      st->x0[(i<<1)+1]=0;      st->x1[(i<<1)+1]=0;   }   /* Reconstruct the original */   fir_mem(st->x0, h0, st->y0, st->full_frame_size, QMF_ORDER, st->g0_mem);   fir_mem(st->x1, h1, st->y1, st->full_frame_size, QMF_ORDER, st->g1_mem);   for (i=0;i<st->full_frame_size;i++)      in[i]=2*(st->y0[i]-st->y1[i]);#endif   for (i=0;i<st->lpcSize;i++)      st->old_lsp[i] = st->lsp[i];   for (i=0;i<st->lpcSize;i++)      st->old_qlsp[i] = st->qlsp[i];   st->first=0;}