fixed-point: computation of rms values in fp
authorjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 8 Oct 2003 04:52:27 +0000 (04:52 +0000)
committerjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 8 Oct 2003 04:52:27 +0000 (04:52 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@5429 0101bb08-14d6-0310-b084-bc0e0c8e3800

libspeex/cb_search.c
libspeex/filters.c
libspeex/filters.h
libspeex/nb_celp.c
libspeex/sb_celp.c
libspeex/smallft.c

index b48c5a7..0ee86ba 100644 (file)
@@ -314,7 +314,7 @@ int   complexity
       }
 #else
       for (j=0;j<subvect_size;j++)
-         e[subvect_size*i+j]=sign*0.03125*SIG_SCALING*shape_cb[rind*subvect_size+j];
+         e[subvect_size*i+j]=sign*0.03125*shape_cb[rind*subvect_size+j];
 #endif
    }   
    /* Update excitation */
index 13dc96b..28b9f1c 100644 (file)
@@ -49,6 +49,47 @@ void bw_lpc(float gamma, spx_coef_t *lpc_in, spx_coef_t *lpc_out, int order)
 
 #ifdef FIXED_POINT
 
+spx_word16_t compute_rms(spx_sig_t *x, int len)
+{
+   int i;
+   spx_word32_t sum=0;
+   spx_sig_t max_val=1;
+   int sig_shift;
+
+   for (i=0;i<len;i++)
+   {
+      spx_sig_t tmp = x[i];
+      if (tmp<0)
+         tmp = -1;
+      if (tmp > max_val)
+         max_val = tmp;
+   }
+
+   sig_shift=0;
+   while (max_val>2048)
+   {
+      sig_shift++;
+      max_val >>= 1;
+   }
+
+   for (i=0;i<len;i+=4)
+   {
+      spx_word32_t sum2=0;
+      spx_word16_t tmp;
+      tmp = SHR(x[i],sig_shift);
+      sum2 += MULT16_16(tmp,tmp);
+      tmp = SHR(x[i+1],sig_shift);
+      sum2 += MULT16_16(tmp,tmp);
+      tmp = SHR(x[i+2],sig_shift);
+      sum2 += MULT16_16(tmp,tmp);
+      tmp = SHR(x[i+3],sig_shift);
+      sum2 += MULT16_16(tmp,tmp);
+      sum += sum2;
+   }
+   
+   /*FIXME: remove division*/
+   return (1<<sig_shift)*sqrt(1+sum/len)/SIG_SCALING;
+}
 
 #define MUL_16_32_R15(a,bh,bl) ((a)*(bh) + ((a)*(bl)>>15))
 
@@ -123,6 +164,16 @@ void fir_mem2(spx_sig_t *x, spx_coef_t *num, spx_sig_t *y, int N, int ord, spx_m
 #include "filters_sse.h"
 #else
 
+spx_word16_t compute_rms(spx_sig_t *x, int len)
+{
+   int i;
+   float sum=0;
+   for (i=0;i<len;i++)
+   {
+      sum += x[i]*x[i];
+   }
+   return sqrt(.1+sum/len);
+}
 
 void filter_mem2(spx_sig_t *x, spx_coef_t *num, spx_coef_t *den, spx_sig_t *y, int N, int ord,  spx_mem_t *mem)
 {
index 02dd548..a241aa4 100644 (file)
@@ -35,6 +35,8 @@
 
 #include "misc.h"
 
+spx_word16_t compute_rms(spx_sig_t *x, int len);
+
 typedef struct CombFilterMem {
    int   last_pitch;
    float last_pitch_gain[3];
index 5a13688..923dfbd 100644 (file)
@@ -350,10 +350,11 @@ int nb_encode(void *state, float *in, SpeexBits *bits)
       if (st->lbr_48k)
       {
          float ol1=0,ol2=0;
-         for (i=0;i<st->frameSize>>1;i++)
-            ol1 += st->exc[i]*st->exc[i] / (SIG_SCALING*SIG_SCALING);
-         for (i=st->frameSize>>1;i<st->frameSize;i++)
-            ol2 += st->exc[i]*st->exc[i] / (SIG_SCALING*SIG_SCALING);
+         
+         ol1 = compute_rms(st->exc, st->frameSize>>1);
+         ol2 = compute_rms(st->exc+(st->frameSize>>1), st->frameSize>>1);
+         ol1 *= ol1*(st->frameSize>>1);
+         ol2 *= ol2*(st->frameSize>>1);
 
          ol_gain=ol1;
          if (ol2>ol1)
@@ -364,11 +365,7 @@ int nb_encode(void *state, float *in, SpeexBits *bits)
 
       } else {
 #endif
-      ol_gain=0;
-      for (i=0;i<st->frameSize;i++)
-         ol_gain += st->exc[i]*st->exc[i] / (SIG_SCALING*SIG_SCALING);
-      
-      ol_gain=sqrt(1+ol_gain/st->frameSize);
+         ol_gain = compute_rms(st->exc, st->frameSize);
 #ifdef EPIC_48K
       }
 #endif
@@ -783,9 +780,9 @@ int nb_encode(void *state, float *in, SpeexBits *bits)
             innov[i]=0;
          
          residue_percep_zero(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2, st->buf2, st->subframeSize, st->lpcSize, stack);
-         for (i=0;i<st->subframeSize;i++)
-            ener+=st->buf2[i]*st->buf2[i] / (SIG_SCALING*SIG_SCALING);
-         ener=sqrt(.1+ener/st->subframeSize);
+
+         ener = compute_rms(st->buf2, st->subframeSize);
+
          /*for (i=0;i<st->subframeSize;i++)
             printf ("%f\n", st->buf2[i]/ener);
          */
index d92188d..27884d3 100644 (file)
@@ -356,13 +356,11 @@ int sb_encode(void *state, float *in, SpeexBits *bits)
       }
 
 
-      for (i=0;i<st->frame_size;i++)
-      {
-         /*FIXME: Are the two signals (low, high) in sync? */
-         e_low  += st->x0d[i]* st->x0d[i] / (SIG_SCALING*SIG_SCALING);
-         e_high += st->high[i]* st->high[i] / (SIG_SCALING*SIG_SCALING);
-      }
-      ratio = log((1+e_high)/(1+e_low));
+      /*FIXME: Are the two signals (low, high) in sync? */
+      e_low = compute_rms(st->x0d, st->frame_size);
+      e_high = compute_rms(st->high, st->frame_size);
+      ratio = 2*log((1+e_high)/(1+e_low));
+      
       speex_encoder_ctl(st->st_low, SPEEX_GET_RELATIVE_QUALITY, &st->relative_quality);
       if (ratio<-4)
          ratio=-4;
@@ -533,18 +531,16 @@ int sb_encode(void *state, float *in, SpeexBits *bits)
       /* 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]/SIG_SCALING);
+
+      eh = compute_rms(exc, st->subframeSize);
 
       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(low_innov[offset+i]/SIG_SCALING);
+
+         el = compute_rms(low_innov+offset, st->subframeSize);
 
          /* 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;
          /*print_vec(&g, 1, "gain factor");*/
@@ -565,11 +561,10 @@ int sb_encode(void *state, float *in, SpeexBits *bits)
       } else {
          float gc, scale, scale_1;
 
-         for (i=0;i<st->subframeSize;i++)
-            el+=sqr(low_exc[offset+i]/SIG_SCALING);
-         /*speex_bits_pack(bits, 0, 1);*/
+         el = compute_rms(low_exc+offset, st->subframeSize);
+         /*FIXME: cleanup the "historical" mess with sqrt(st->subframeSize) */
+         gc = (1+eh)*filter_ratio/(1+el)/sqrt(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)
@@ -580,7 +575,7 @@ int sb_encode(void *state, float *in, SpeexBits *bits)
             gc = exp((1/3.7)*qgc-2);
          }
 
-         scale = gc*sqrt(1+el)/filter_ratio;
+         scale = gc*(1+el*sqrt(st->subframeSize))/filter_ratio;
          scale_1 = 1/scale;
 
          for (i=0;i<st->subframeSize;i++)
@@ -1024,8 +1019,6 @@ int sb_decode(void *state, SpeexBits *bits, float *out)
          float g;
          int quant;
 
-         for (i=0;i<st->subframeSize;i++)
-            el+=sqr(low_innov[offset+i]/SIG_SCALING);
          quant = speex_bits_unpack_unsigned(bits, 5);
          g= exp(((float)quant-10)/8.0);
          
@@ -1039,13 +1032,12 @@ int sb_decode(void *state, SpeexBits *bits, float *out)
       } else {
          float gc, scale;
          int qgc = speex_bits_unpack_unsigned(bits, 4);
-         for (i=0;i<st->subframeSize;i++)
-            el+=sqr(low_exc[offset+i]/SIG_SCALING);
 
+         el = compute_rms(low_exc+offset, st->subframeSize);
 
          gc = exp((1/3.7)*qgc-2);
 
-         scale = gc*sqrt(1+el)/filter_ratio;
+         scale = gc*(1+el*sqrt(st->subframeSize))/filter_ratio;
 
          SUBMODE(innovation_unquant)(exc, SUBMODE(innovation_params), st->subframeSize, 
                                 bits, stack);
index 75b5276..de37eaf 100644 (file)
@@ -11,7 +11,7 @@
  ********************************************************************
 
  function: *unnormalized* fft transform
- last mod: $Id: smallft.c,v 1.9 2003/10/08 04:50:44 jm Exp $
+ last mod: $Id: smallft.c,v 1.10 2003/10/08 04:52:27 jm Exp $
 
  ********************************************************************/