converted theta and prior_ratio to fixed-point.
authorjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 1 Nov 2006 14:12:57 +0000 (14:12 +0000)
committerjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Wed, 1 Nov 2006 14:12:57 +0000 (14:12 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@11984 0101bb08-14d6-0310-b084-bc0e0c8e3800

libspeex/arch.h
libspeex/fixed_generic.h
libspeex/math_approx.c
libspeex/math_approx.h
libspeex/preprocess.c

index 8d40fac..501cd2a 100644 (file)
@@ -151,6 +151,7 @@ typedef float spx_word32_t;
 #define MULT16_32_Q13(a,b)     ((a)*(b))
 #define MULT16_32_Q14(a,b)     ((a)*(b))
 #define MULT16_32_Q15(a,b)     ((a)*(b))
+#define MULT16_32_P15(a,b)     ((a)*(b))
 
 #define MAC16_32_Q11(c,a,b)     ((c)+(a)*(b))
 #define MAC16_32_Q15(c,a,b)     ((c)+(a)*(b))
index 375050c..2345e36 100644 (file)
@@ -77,6 +77,7 @@
 #define MULT16_32_Q11(a,b) ADD32(MULT16_16((a),SHR((b),11)), SHR(MULT16_16((a),((b)&0x000007ff)),11))
 #define MAC16_32_Q11(c,a,b) ADD32(c,ADD32(MULT16_16((a),SHR((b),11)), SHR(MULT16_16((a),((b)&0x000007ff)),11)))
 
+#define MULT16_32_P15(a,b) ADD32(MULT16_16((a),SHR((b),15)), PSHR(MULT16_16((a),((b)&0x00007fff)),15))
 #define MULT16_32_Q15(a,b) ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15))
 #define MAC16_32_Q15(c,a,b) ADD32(c,ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15)))
 
index c008ec2..9bd5f88 100644 (file)
@@ -235,6 +235,7 @@ static spx_word32_t spx_exp2(spx_word16_t x)
       return SHR32(EXTEND32(frac), -integer-2);
 }
 
+/* Input in Q11 format, output in Q16 */
 spx_word32_t spx_exp(spx_word16_t x)
 {
    if (x>21290)
index d8cb6d1..a3c3d6c 100644 (file)
@@ -44,9 +44,11 @@ spx_int16_t spx_ilog4(spx_uint32_t x);
 #ifdef FIXED_POINT
 spx_word16_t spx_sqrt(spx_word32_t x);
 spx_word16_t spx_acos(spx_word16_t x);
+spx_word32_t spx_exp(spx_word16_t x);
 #else
 #define spx_sqrt sqrt
 #define spx_acos acos
+#define spx_exp exp
 #endif
 
 #endif
index c33315f..b072f3b 100644 (file)
@@ -151,6 +151,11 @@ static inline spx_word16_t DIV32_16_Q15(spx_word32_t a, spx_word32_t b)
 #define FRAC_SCALING_1 3.0518e-05
 #define FRAC_SHIFT 1
 
+#define EXPIN_SCALING 2048.f
+#define EXPIN_SCALING_1 0.00048828f
+#define EXPIN_SHIFT 11
+#define EXPOUT_SCALING_1 1.5259e-05
+
 #define NOISE_SHIFT 7
 
 #else
@@ -165,6 +170,10 @@ static inline spx_word16_t DIV32_16_Q15(spx_word32_t a, spx_word32_t b)
 #define FRAC_SHIFT 0
 #define NOISE_SHIFT 0
 
+#define EXPIN_SCALING 1.f
+#define EXPIN_SCALING_1 1.f
+#define EXPOUT_SCALING_1 1.f
+
 #endif
 
 /** Speex pre-processor state. */
@@ -717,33 +726,36 @@ int speex_preprocess_run(SpeexPreprocessState *st, spx_int16_t *x)
    
    noise_floor = exp(.2302585f*st->noise_suppress);
    echo_floor = exp(.2302585f* (st->echo_suppress*(1-Pframe) + st->echo_suppress_active*Pframe));
-   /*print_vec(&Pframe, 1, "");*/
-   /*for (i=N;i<N+M;i++)
-      st->gain2[i] = qcurve (st->zeta[i]);
-   filterbank_compute_psd(st->bank,st->gain2+N, st->gain2);*/
    
+   /* Compute Ephraim & Malah gain speech probability of presence for each critical band (Bark scale) 
+      Technically this is actually wrong because the EM gaim assumes a slightly different probability 
+      distribution */
    for (i=N;i<N+M;i++)
    {
-      float theta, MM;
-      float prior_ratio;
+      spx_word32_t theta;
+      float MM;
+      spx_word16_t prior_ratio;
       float q;
       float P1;
 
       /* Compute the gain floor based on different floors for the background noise and residual echo */
       st->gain_floor[i] = FRAC_SCALING*sqrt((noise_floor*PSHR32(st->noise[i],NOISE_SHIFT) + echo_floor*st->echo_noise[i])/(1+PSHR32(st->noise[i],NOISE_SHIFT) + st->echo_noise[i]));
-      prior_ratio = st->prior[i]/(SNR_SCALING+st->prior[i]);
-      theta = (1.f+SNR_SCALING_1*st->post[i])*prior_ratio;
+      prior_ratio = FRAC_SCALING*st->prior[i]/(SNR_SCALING+st->prior[i]);
+      theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
 
-      MM = hypergeom_gain(theta);
+      MM = hypergeom_gain(EXPIN_SCALING_1*theta);
       /* Gain with bound */
-      st->gain[i] = MIN16(FRAC_SCALING, FRAC_SCALING*prior_ratio * MM);
+      st->gain[i] = MIN16(FRAC_SCALING, prior_ratio * MM);
       
       /* Save old Bark power spectrum */
       st->old_ps[i] = .2*st->old_ps[i] + .8*FRAC_SCALING_1*FRAC_SCALING_1*st->gain[i]*st->gain[i]*ps[i];
 
       P1 = .2+.8*qcurve (st->zeta[i]);
       q = 1-Pframe*P1;
-      st->gain2[i]=FRAC_SCALING/(1.f + (q/(1.f-q))*(1.f+SNR_SCALING_1*st->prior[i])*exp(-theta));
+#ifdef FIXED_POINT
+      theta = MIN32(theta, 32767);
+#endif
+      st->gain2[i]=FRAC_SCALING/(1.f + (q/(1.f-q))*(1.f+SNR_SCALING_1*st->prior[i])*EXPOUT_SCALING_1*spx_exp(-EXTRACT16(theta)));
    }
    filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
    filterbank_compute_psd16(st->bank,st->gain+N, st->gain);