[xiph-commits] r10567 - trunk/speex/libspeex

[jm at svn.xiph.org]

Author: jm
Date: 2005-12-08 04:16:12 -0800 (Thu, 08 Dec 2005)
New Revision: 10567

Modified:
   trunk/speex/libspeex/mdf.c
   trunk/speex/libspeex/pseudofloat.h
Log:
debugged pseudo-float type, using it for correlations


Modified: trunk/speex/libspeex/mdf.c
===================================================================
--- trunk/speex/libspeex/mdf.c	2005-12-08 09:28:10 UTC (rev 10566)
+++ trunk/speex/libspeex/mdf.c	2005-12-08 12:16:12 UTC (rev 10567)
@@ -44,9 +44,8 @@
 //#include "speex/speex_echo.h"
 #include "smallft.h"
 #include "fftwrap.h"
+#include "pseudofloat.h"
 
-#include <math.h>
-
 #ifndef M_PI
 #define M_PI 3.14159265358979323846
 #endif
@@ -292,7 +291,7 @@
    float Syy=0,See=0;
    float leak_estimate;
    float ss;
-   float adapt_rate;
+   float adapt_rate=0;
    
    N = st->window_size;
    M = st->M;
@@ -361,25 +360,23 @@
       st->power[j] = (1-ss)*st->power[j] + 1 + ss*st->Xf[j];
 
    {
-      float Pey = 0, Pyy=0;
+      spx_float_t Pey = FLOAT_ZERO, Pyy=FLOAT_ZERO;
       float alpha;
-      for (j=0;j<=st->frame_size;j++)
+      for (j=st->frame_size;j>=0;j--)
       {
-         spx_word32_t E, Y, Eh, Yh;
-         E = (st->Rf[j]);
-         Y = (st->Yf[j]);
-         Eh = st->Eh[j] + E;
-         Yh = st->Yh[j] + Y;
-         Pey += Eh*1.f*Yh;
-         Pyy += Yh*1.f*Yh;
-         st->Eh[j] = .95*Eh - E;
-         st->Yh[j] = .95*Yh - Y;
+         spx_float_t Eh, Yh;
+         Eh = PSEUDOFLOAT(st->Rf[j] - st->Eh[j]);
+         Yh = PSEUDOFLOAT(st->Yf[j] - st->Yh[j]);
+         Pey = FLOAT_ADD(Pey,FLOAT_MULT(Eh,Yh));
+         Pyy = FLOAT_ADD(Pyy,FLOAT_MULT(Yh,Yh));
+         st->Eh[j] = .95*st->Eh[j] + .05*st->Rf[j];
+         st->Yh[j] = .95*st->Yh[j] + .05*st->Yf[j];
       }
       alpha = .02*Syy / (1e4+See);
       if (alpha > .02)
          alpha = .02;
-      st->Pey = (1-alpha)*st->Pey + alpha*Pey;
-      st->Pyy = (1-alpha)*st->Pyy + alpha*Pyy;
+      st->Pey = (1-alpha)*st->Pey + alpha*REALFLOAT(Pey);
+      st->Pyy = (1-alpha)*st->Pyy + alpha*REALFLOAT(Pyy);
       if (st->Pey< .001*st->Pyy)
          st->Pey = .001*st->Pyy;
       leak_estimate = st->Pey / (1+st->Pyy);

Modified: trunk/speex/libspeex/pseudofloat.h
===================================================================
--- trunk/speex/libspeex/pseudofloat.h	2005-12-08 09:28:10 UTC (rev 10566)
+++ trunk/speex/libspeex/pseudofloat.h	2005-12-08 12:16:12 UTC (rev 10567)
@@ -36,53 +36,95 @@
 #define PSEUDOFLOAT_H
 
 #include "misc.h"
+#include <math.h>
 
+#ifdef FIXED_POINT
+
 typedef struct {
    spx_int16_t m;
    spx_int16_t e;
 } spx_float_t;
 
-/*#define FLOAT_MULT(a,b) ((spx_float_t) {(spx_int16_t)((spx_int32_t)(a).m*(b).m>>15), (a).e+(b).e+15})
+#define FLOAT_ZERO {0,0}
 
-#define FLOAT_ADD(a,b) ( (a).e > (b).e ? (spx_float_t) {((a).m>>1) + ((a).m>>((a).e-(b).e+1)),(a).e+1} : (spx_float_t) {((b).m>>1) + ((a).m>>((b).e-(a).e+1)),(b).e+1})*/
-
+#define MIN(a,b) ((a)<(b)?(a):(b))
 static inline spx_float_t PSEUDOFLOAT(float x)
 {
    int e=0;
+   int sign=0;
+   if (x<0)
+   {
+      sign = 1;
+      x = -x;
+   }
    if (x==0)
       return (spx_float_t) {0,0};
    while (x>32767)
    {
+      //x >>= 1;
       x *= .5;
       e++;
    }
    while (x<16383)
    {
+      //x <<= 1;
       x *= 2;
       e--;
    }
-   return (spx_float_t) {x,e};
+   if (sign)
+      return (spx_float_t) {-x,e};
+   else      
+      return (spx_float_t) {x,e};
 }
 
+static inline float REALFLOAT(spx_float_t a)
+{
+   return a.m * pow(2,a.e);
+}
+
 static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b)
 {
-   spx_float_t r = (a).e > (b).e ? (spx_float_t) {((a).m>>1) + ((a).m>>((a).e-(b).e+1)),(a).e+1} : (spx_float_t) {((b).m>>1) + ((a).m>>((b).e-(a).e+1)),(b).e+1};
-   if (r.m<16384)
+   if (a.m==0)
+      return b;
+   else if (b.m==0)
+      return a;
+   spx_float_t r = (a).e > (b).e ? (spx_float_t) {((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1)),(a).e+1} : (spx_float_t) {((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1)),(b).e+1};
+   if (r.m>0)
    {
-      r.m<<=1;
-      r.e-=1;
+      if (r.m<16384)
+      {
+         r.m<<=1;
+         r.e-=1;
+      }
+   } else {
+      if (r.m>-16384)
+      {
+         r.m<<=1;
+         r.e-=1;
+      }
    }
+   //printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));
    return r;
 }
 
 static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b)
 {
    spx_float_t r = (spx_float_t) {(spx_int16_t)((spx_int32_t)(a).m*(b).m>>15), (a).e+(b).e+15};
-   if (r.m<16384)
+   if (r.m>0)
    {
-      r.m<<=1;
-      r.e-=1;
+      if (r.m<16384)
+      {
+         r.m<<=1;
+         r.e-=1;
+      }
+   } else {
+      if (r.m>-16384)
+      {
+         r.m<<=1;
+         r.e-=1;
+      }
    }
+   //printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));
    return r;   
 }
 
@@ -144,4 +186,15 @@
    return (spx_float_t) {DIV32_16(a,b),e};
 }
 
+#else
+
+#define spx_float_t float
+#define FLOAT_ZERO 0.f
+#define PSEUDOFLOAT(x) (x)
+#define FLOAT_MULT(a,b) ((a)*(b))
+#define FLOAT_ADD(a,b) ((a)+(b))
+#define REALFLOAT(x) (x)
+
 #endif
+
+#endif