filters.c

一个开源SIP协议栈

void fir_mem2(const spx_sig_t *x, const spx_coef_t *num, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
{
   int i,j;
   spx_word32_t xi,yi;

   for (i=0;i<ord;i++)
      mem[i] = SHR32(mem[i],1);   
   for (i=0;i<N;i++)
   {
      xi=SATURATE(x[i],805306368);
      yi = xi + SHL32(mem[0],2);
      for (j=0;j<ord-1;j++)
      {
         mem[j] = MAC16_32_Q15(mem[j+1], num[j],xi);
      }
      mem[ord-1] = MULT16_32_Q15(num[ord-1],xi);
      y[i] = SATURATE(yi,805306368);
   }
   for (i=0;i<ord;i++)
      mem[i] = SHL32(mem[i],1);   
}
#endif
#endif

#ifdef FIXED_POINT
#ifndef OVERRIDE_FIR_MEM16
void fir_mem16(const spx_word16_t *x, const spx_coef_t *num, spx_word16_t *y, int N, int ord, spx_mem_t *mem, char *stack)
{
   int i,j;
   spx_word16_t xi,yi;

   for (i=0;i<N;i++)
   {
      xi=x[i];
      yi = EXTRACT16(SATURATE(ADD32(EXTEND32(x[i]),PSHR32(mem[0],LPC_SHIFT)),32767));
      for (j=0;j<ord-1;j++)
      {
         mem[j] = MAC16_16(mem[j+1], num[j],xi);
      }
      mem[ord-1] = MULT16_16(num[ord-1],xi);
      y[i] = yi;
   }
}
#endif
#else
void fir_mem16(const spx_word16_t *x, const spx_coef_t *num, spx_word16_t *y, int N, int ord, spx_mem_t *mem, char *stack)
{
   fir_mem2(x, num, y, N, ord, mem);
}
#endif






void syn_percep_zero(const spx_sig_t *xx, const spx_coef_t *ak, const spx_coef_t *awk1, const spx_coef_t *awk2, spx_sig_t *y, int N, int ord, char *stack)
{
   int i;
   VARDECL(spx_mem_t *mem);
   ALLOC(mem, ord, spx_mem_t);
   for (i=0;i<ord;i++)
     mem[i]=0;
   iir_mem2(xx, ak, y, N, ord, mem);
   for (i=0;i<ord;i++)
      mem[i]=0;
   filter_mem2(y, awk1, awk2, y, N, ord, mem);
}

void residue_percep_zero(const spx_sig_t *xx, const spx_coef_t *ak, const spx_coef_t *awk1, const spx_coef_t *awk2, spx_sig_t *y, int N, int ord, char *stack)
{
   int i;
   VARDECL(spx_mem_t *mem);
   ALLOC(mem, ord, spx_mem_t);
   for (i=0;i<ord;i++)
      mem[i]=0;
   filter_mem2(xx, ak, awk1, y, N, ord, mem);
   for (i=0;i<ord;i++)
     mem[i]=0;
   fir_mem2(y, awk2, y, N, ord, mem);
}

#ifndef OVERRIDE_COMPUTE_IMPULSE_RESPONSE
void compute_impulse_response(const spx_coef_t *ak, const spx_coef_t *awk1, const spx_coef_t *awk2, spx_word16_t *y, int N, int ord, char *stack)
{
   int i,j;
   spx_word16_t y1, ny1i, ny2i;
   VARDECL(spx_mem_t *mem1);
   VARDECL(spx_mem_t *mem2);
   ALLOC(mem1, ord, spx_mem_t);
   ALLOC(mem2, ord, spx_mem_t);
   
   y[0] = LPC_SCALING;
   for (i=0;i<ord;i++)
      y[i+1] = awk1[i];
   i++;
   for (;i<N;i++)
      y[i] = VERY_SMALL;
   for (i=0;i<ord;i++)
      mem1[i] = mem2[i] = 0;
   for (i=0;i<N;i++)
   {
      y1 = ADD16(y[i], EXTRACT16(PSHR32(mem1[0],LPC_SHIFT)));
      ny1i = NEG16(y1);
      y[i] = PSHR32(ADD32(SHL32(EXTEND32(y1),LPC_SHIFT+1),mem2[0]),LPC_SHIFT);
      ny2i = NEG16(y[i]);
      for (j=0;j<ord-1;j++)
      {
         mem1[j] = MAC16_16(mem1[j+1], awk2[j],ny1i);
         mem2[j] = MAC16_16(mem2[j+1], ak[j],ny2i);
      }
      mem1[ord-1] = MULT16_16(awk2[ord-1],ny1i);
      mem2[ord-1] = MULT16_16(ak[ord-1],ny2i);
   }
}
#endif

void qmf_decomp(const spx_word16_t *xx, const spx_word16_t *aa, spx_sig_t *y1, spx_sig_t *y2, int N, int M, spx_word16_t *mem, char *stack)
{
   int i,j,k,M2;
   VARDECL(spx_word16_t *a);
   VARDECL(spx_word16_t *x);
   spx_word16_t *x2;
   
   ALLOC(a, M, spx_word16_t);
   ALLOC(x, N+M-1, spx_word16_t);
   x2=x+M-1;
   M2=M>>1;
   for (i=0;i<M;i++)
      a[M-i-1]= aa[i];

   for (i=0;i<M-1;i++)
      x[i]=mem[M-i-2];
   for (i=0;i<N;i++)
      x[i+M-1]=SATURATE(PSHR(xx[i],1),16383);
   for (i=0,k=0;i<N;i+=2,k++)
   {
      y1[k]=0;
      y2[k]=0;
      for (j=0;j<M2;j++)
      {
         y1[k]=ADD32(y1[k],MULT16_16(a[j],ADD16(x[i+j],x2[i-j])));
         y2[k]=SUB32(y2[k],MULT16_16(a[j],SUB16(x[i+j],x2[i-j])));
         j++;
         y1[k]=ADD32(y1[k],MULT16_16(a[j],ADD16(x[i+j],x2[i-j])));
         y2[k]=ADD32(y2[k],MULT16_16(a[j],SUB16(x[i+j],x2[i-j])));
      }
      y1[k] = SHR32(y1[k],1);
      y2[k] = SHR32(y2[k],1);
   }
   for (i=0;i<M-1;i++)
     mem[i]=SATURATE(PSHR(xx[N-i-1],1),16383);
}


/* By segher */
void fir_mem_up(const spx_sig_t *x, const spx_word16_t *a, spx_sig_t *y, int N, int M, spx_word32_t *mem, char *stack)
   /* assumptions:
      all odd x[i] are zero -- well, actually they are left out of the array now
      N and M are multiples of 4 */
{
   int i, j;
   VARDECL(spx_word16_t *xx);
   
   ALLOC(xx, M+N-1, spx_word16_t);

   for (i = 0; i < N/2; i++)
      xx[2*i] = PSHR32(x[N/2-1-i],SIG_SHIFT);
   for (i = 0; i < M - 1; i += 2)
      xx[N+i] = mem[i+1];

   for (i = 0; i < N; i += 4) {
      spx_sig_t y0, y1, y2, y3;
      spx_word16_t x0;

      y0 = y1 = y2 = y3 = 0;
      x0 = xx[N-4-i];

      for (j = 0; j < M; j += 4) {
         spx_word16_t x1;
         spx_word16_t a0, a1;

         a0 = a[j];
         a1 = a[j+1];
         x1 = xx[N-2+j-i];

         y0 = ADD32(y0,SHR32(MULT16_16(a0, x1),2));
         y1 = ADD32(y1,SHR32(MULT16_16(a1, x1),2));
         y2 = ADD32(y2,SHR32(MULT16_16(a0, x0),2));
         y3 = ADD32(y3,SHR32(MULT16_16(a1, x0),2));

         a0 = a[j+2];
         a1 = a[j+3];
         x0 = xx[N+j-i];

         y0 = ADD32(y0,SHR32(MULT16_16(a0, x0),2));
         y1 = ADD32(y1,SHR32(MULT16_16(a1, x0),2));
         y2 = ADD32(y2,SHR32(MULT16_16(a0, x1),2));
         y3 = ADD32(y3,SHR32(MULT16_16(a1, x1),2));
      }
      y[i] = y0;
      y[i+1] = y1;
      y[i+2] = y2;
      y[i+3] = y3;
   }

   for (i = 0; i < M - 1; i += 2)
      mem[i+1] = xx[i];
}

#ifdef FIXED_POINT
#if 0
const spx_word16_t shift_filt[3][7] = {{-33,    1043,   -4551,   19959,   19959,   -4551,    1043},
                                 {-98,    1133,   -4425,   29179,    8895,   -2328,     444},
                                 {444,   -2328,    8895,   29179,   -4425,    1133,     -98}};
#else
const spx_word16_t shift_filt[3][7] = {{-390,    1540,   -4993,   20123,   20123,   -4993,    1540},
                                {-1064,    2817,   -6694,   31589,    6837,    -990,    -209},
                                 {-209,    -990,    6837,   31589,   -6694,    2817,   -1064}};
#endif
#else
#if 0
const float shift_filt[3][7] = {{-9.9369e-04, 3.1831e-02, -1.3889e-01, 6.0910e-01, 6.0910e-01, -1.3889e-01, 3.1831e-02},
                          {-0.0029937, 0.0345613, -0.1350474, 0.8904793, 0.2714479, -0.0710304, 0.0135403},
                          {0.0135403, -0.0710304, 0.2714479, 0.8904793, -0.1350474, 0.0345613,  -0.0029937}};
#else
const float shift_filt[3][7] = {{-0.011915f, 0.046995f, -0.152373f, 0.614108f, 0.614108f, -0.152373f, 0.046995f},
                          {-0.0324855f, 0.0859768f, -0.2042986f, 0.9640297f, 0.2086420f, -0.0302054f, -0.0063646f},
                          {-0.0063646f, -0.0302054f, 0.2086420f, 0.9640297f, -0.2042986f, 0.0859768f, -0.0324855f}};
#endif
#endif

int interp_pitch(
spx_word16_t *exc,          /*decoded excitation*/
spx_word16_t *interp,          /*decoded excitation*/
int pitch,               /*pitch period*/
int len
)
{
   int i,j,k;
   spx_word32_t corr[4][7];
   spx_word32_t maxcorr;
   int maxi, maxj;
   for (i=0;i<7;i++)
   {
      corr[0][i] = inner_prod(exc, exc-pitch-3+i, len);
   }
   for (i=0;i<3;i++)
   {
      for (j=0;j<7;j++)
      {
         int i1, i2;
         spx_word32_t tmp=0;
         i1 = 3-j;
         if (i1<0)
            i1 = 0;
         i2 = 10-j;
         if (i2>7)
            i2 = 7;
         for (k=i1;k<i2;k++)
            tmp += MULT16_32_Q15(shift_filt[i][k],corr[0][j+k-3]);
         corr[i+1][j] = tmp;
      }
   }
   maxi=maxj=0;
   maxcorr = corr[0][0];
   for (i=0;i<4;i++)
   {
      for (j=0;j<7;j++)
      {
         if (corr[i][j] > maxcorr)
         {
            maxcorr = corr[i][j];
            maxi=i;
            maxj=j;
         }
      }
   }
   for (i=0;i<len;i++)
   {
      spx_word32_t tmp = 0;
      if (maxi>0)
      {
         for (k=0;k<7;k++)
         {
            tmp += MULT16_16(exc[i-(pitch-maxj+3)+k-3],shift_filt[maxi-1][k]);
         }
      } else {
         tmp = SHL32(exc[i-(pitch-maxj+3)],15);
      }
      interp[i] = PSHR32(tmp,15);
   }
   return pitch-maxj+3;
}

void multicomb(
spx_word16_t *exc,          /*decoded excitation*/
spx_word16_t *new_exc,      /*enhanced excitation*/
spx_coef_t *ak,           /*LPC filter coefs*/
int p,               /*LPC order*/
int nsf,             /*sub-frame size*/
int pitch,           /*pitch period*/
int max_pitch,
spx_word16_t  comb_gain,    /*gain of comb filter*/
char *stack
)
{
   int i; 
   VARDECL(spx_word16_t *iexc);
   spx_word16_t old_ener, new_ener;
   int corr_pitch;
   
   spx_word16_t iexc0_mag, iexc1_mag, exc_mag;
   spx_word32_t corr0, corr1;
   spx_word16_t gain0, gain1;
   spx_word16_t pgain1, pgain2;
   spx_word16_t c1, c2;
   spx_word16_t g1, g2;
   spx_word16_t ngain;
   spx_word16_t gg1, gg2;
#ifdef FIXED_POINT
   int scaledown=0;
#endif
#if 0 /* Set to 1 to enable full pitch search */
   int nol_pitch[6];
   spx_word16_t nol_pitch_coef[6];
   spx_word16_t ol_pitch_coef;
   open_loop_nbest_pitch(exc, 20, 120, nsf, 
                         nol_pitch, nol_pitch_coef, 6, stack);
   corr_pitch=nol_pitch[0];
   ol_pitch_coef = nol_pitch_coef[0];
   /*Try to remove pitch multiples*/
   for (i=1;i<6;i++)
   {
#ifdef FIXED_POINT
      if ((nol_pitch_coef[i]>MULT16_16_Q15(nol_pitch_coef[0],19661)) && 
#else
      if ((nol_pitch_coef[i]>.6*nol_pitch_coef[0]) && 
#endif
         (ABS(2*nol_pitch[i]-corr_pitch)<=2 || ABS(3*nol_pitch[i]-corr_pitch)<=3 || 
         ABS(4*nol_pitch[i]-corr_pitch)<=4 || ABS(5*nol_pitch[i]-corr_pitch)<=5))
      {
         corr_pitch = nol_pitch[i];
      }
   }
#else
   corr_pitch = pitch;
#endif
   
   ALLOC(iexc, 2*nsf, spx_word16_t);
   
   interp_pitch(exc, iexc, corr_pitch, 80);
   if (corr_pitch>max_pitch)
      interp_pitch(exc, iexc+nsf, 2*corr_pitch, 80);
   else
      interp_pitch(exc, iexc+nsf, -corr_pitch, 80);

#ifdef FIXED_POINT
   for (i=0;i<nsf;i++)
   {
      if (ABS16(exc[i])>16383)
      {
         scaledown = 1;
         break;
      }
   }
   if (scaledown)
   {
      for (i=0;i<nsf;i++)
         exc[i] = SHR16(exc[i],1);
      for (i=0;i<2*nsf;i++)
         iexc[i] = SHR16(iexc[i],1);
   }
#endif
   /*interp_pitch(exc, iexc+2*nsf, 2*corr_pitch, 80);*/
   
   /*printf ("%d %d %f\n", pitch, corr_pitch, max_corr*ener_1);*/
   iexc0_mag = spx_sqrt(1000+inner_prod(iexc,iexc,nsf));
   iexc1_mag = spx_sqrt(1000+inner_prod(iexc+nsf,iexc+nsf,nsf));
   exc_mag = spx_sqrt(1+inner_prod(exc,exc,nsf));
   corr0  = inner_prod(iexc,exc,nsf);
   if (corr0<0)
      corr0=0;
   corr1 = inner_prod(iexc+nsf,exc,nsf);
   if (corr1<0)
      corr1=0;
#ifdef FIXED_POINT
   /* Doesn't cost much to limit the ratio and it makes the rest easier */
   if (SHL32(EXTEND32(iexc0_mag),6) < EXTEND32(exc_mag))
      iexc0_mag = ADD16(1,PSHR16(exc_mag,6));
   if (SHL32(EXTEND32(iexc1_mag),6) < EXTEND32(exc_mag))
      iexc1_mag = ADD16(1,PSHR16(exc_mag,6));
#endif
   if (corr0 > MULT16_16(iexc0_mag,exc_mag))
      pgain1 = QCONST16(1., 14);
   else
      pgain1 = PDIV32_16(SHL32(PDIV32(corr0, exc_mag),14),iexc0_mag);
   if (corr1 > MULT16_16(iexc1_mag,exc_mag))
      pgain2 = QCONST16(1., 14);
   else
      pgain2 = PDIV32_16(SHL32(PDIV32(corr1, exc_mag),14),iexc1_mag);
   gg1 = PDIV32_16(SHL32(EXTEND32(exc_mag),8), iexc0_mag);
   gg2 = PDIV32_16(SHL32(EXTEND32(exc_mag),8), iexc1_mag);
   if (comb_gain>0)
   {
#ifdef FIXED_POINT
      c1 = (MULT16_16_Q15(QCONST16(.4,15),comb_gain)+QCONST16(.07,15));
      c2 = QCONST16(.5,15)+MULT16_16_Q14(QCONST16(1.72,14),(c1-QCONST16(.07,15)));
#else
      c1 = .4*comb_gain+.07;
      c2 = .5+1.72*(c1-.07);
#endif
   } else 
   {
      c1=c2=0;
   }
#ifdef FIXED_POINT
   g1 = 32767 - MULT16_16_Q13(MULT16_16_Q15(c2, pgain1),pgain1);
   g2 = 32767 - MULT16_16_Q13(MULT16_16_Q15(c2, pgain2),pgain2);
#else
   g1 = 1-c2*pgain1*pgain1;
   g2 = 1-c2*pgain2*pgain2;
#endif
   if (g1<c1)
      g1 = c1;
   if (g2<c1)
      g2 = c1;
   g1 = (spx_word16_t)PDIV32_16(SHL32(EXTEND32(c1),14),(spx_word16_t)g1);
   g2 = (spx_word16_t)PDIV32_16(SHL32(EXTEND32(c1),14),(spx_word16_t)g2);
   if (corr_pitch>max_pitch)
   {
      gain0 = MULT16_16_Q15(QCONST16(.7,15),MULT16_16_Q14(g1,gg1));
      gain1 = MULT16_16_Q15(QCONST16(.3,15),MULT16_16_Q14(g2,gg2));
   } else {
      gain0 = MULT16_16_Q15(QCONST16(.6,15),MULT16_16_Q14(g1,gg1));
      gain1 = MULT16_16_Q15(QCONST16(.6,15),MULT16_16_Q14(g2,gg2));
   }
   for (i=0;i<nsf;i++)
      new_exc[i] = ADD16(exc[i], EXTRACT16(PSHR32(ADD32(MULT16_16(gain0,iexc[i]), MULT16_16(gain1,iexc[i+nsf])),8)));
   /* FIXME: compute_rms16 is currently not quite accurate enough (but close) */
   new_ener = compute_rms16(new_exc, nsf);
   old_ener = compute_rms16(exc, nsf);
   
   if (old_ener < 1)
      old_ener = 1;
   if (new_ener < 1)
      new_ener = 1;
   if (old_ener > new_ener)
      old_ener = new_ener;
   ngain = PDIV32_16(SHL32(EXTEND32(old_ener),14),new_ener);
   
   for (i=0;i<nsf;i++)
      new_exc[i] = MULT16_16_Q14(ngain, new_exc[i]);
#ifdef FIXED_POINT
   if (scaledown)
   {
      for (i=0;i<nsf;i++)
         exc[i] = SHL16(exc[i],1);
      for (i=0;i<nsf;i++)
         new_exc[i] = SHL16(SATURATE16(new_exc[i],16383),1);
   }
#endif
}