[xiph-commits] r10559 - in trunk/speex: include/speex libspeex
jm at svn.xiph.org
jm at svn.xiph.org
Wed Dec 7 06:16:29 PST 2005
Author: jm
Date: 2005-12-07 06:16:24 -0800 (Wed, 07 Dec 2005)
New Revision: 10559
Modified:
trunk/speex/include/speex/speex_echo.h
trunk/speex/libspeex/mdf.c
Log:
Most exchanged variables are now integers. Still need to convert the functions
and test it more.
Modified: trunk/speex/include/speex/speex_echo.h
===================================================================
--- trunk/speex/include/speex/speex_echo.h 2005-12-07 04:31:05 UTC (rev 10558)
+++ trunk/speex/include/speex/speex_echo.h 2005-12-07 14:16:24 UTC (rev 10559)
@@ -39,42 +39,9 @@
#endif
/*struct drft_lookup;*/
+typedef SpeexEchoState;
-/** Speex echo cancellation state. */
-typedef struct SpeexEchoState {
- int frame_size; /**< Number of samples processed each time */
- int window_size;
- int M;
- int cancel_count;
- int adapted;
- float sum_adapt;
- float *e;
- float *x;
- float *X;
- float *d;
- float *y;
- float *last_y;
- float *Yps;
- float *Y;
- float *E;
- float *PHI;
- float *W;
- float *power;
- float *power_1;
- float *Rf;
- float *Yf;
- float *Xf;
- float *Eh;
- float *Yh;
- float Pey;
- float Pyy;
- /*struct drft_lookup *fft_lookup;*/
- void *fft_table;
-
-} SpeexEchoState;
-
-
/** Creates a new echo canceller state */
SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length);
Modified: trunk/speex/libspeex/mdf.c
===================================================================
--- trunk/speex/libspeex/mdf.c 2005-12-07 04:31:05 UTC (rev 10558)
+++ trunk/speex/libspeex/mdf.c 2005-12-07 14:16:24 UTC (rev 10559)
@@ -38,8 +38,10 @@
#include "config.h"
#endif
+#define MDF_C
+
#include "misc.h"
-#include "speex/speex_echo.h"
+//#include "speex/speex_echo.h"
#include "smallft.h"
#include "fftwrap.h"
@@ -56,59 +58,132 @@
#define max(a,b) ((a)>(b) ? (a) : (b))
#ifdef FIXED_POINT
-#define WEIGHT_SCALING 128.f
-#define WEIGHT_SCALING_1 0.0078125f
+#define WEIGHT_SHIFT 11
+#define WEIGHT_SCALING 2048
+#define WEIGHT_SCALING_1 0.00048828f
+//#define WEIGHT_SCALING (100*16*128.f)
+//#define WEIGHT_SCALING_1 (.0625f*0.0078125f)
#else
#define WEIGHT_SCALING 1.f
#define WEIGHT_SCALING_1 1.f
#endif
+/** Speex echo cancellation state. */
+typedef struct {
+ int frame_size; /**< Number of samples processed each time */
+ int window_size;
+ int M;
+ int cancel_count;
+ int adapted;
+ float sum_adapt;
+ spx_word16_t *e;
+ spx_word16_t *x;
+ spx_word16_t *X;
+ spx_word16_t *d;
+ spx_word16_t *y;
+ spx_word16_t *last_y;
+ spx_word32_t *Yps;
+ spx_word16_t *Y;
+ spx_word16_t *E;
+ spx_word16_t *PHI;
+ spx_word16_t *W;
+ spx_word32_t *power;
+ float *power_1;
+ spx_word32_t *Rf;
+ spx_word32_t *Yf;
+ spx_word32_t *Xf;
+ spx_word32_t *Eh;
+ spx_word32_t *Yh;
+ float Pey;
+ float Pyy;
+ /*struct drft_lookup *fft_lookup;*/
+ void *fft_table;
+
+
+} SpeexEchoState;
+
+
/** Compute inner product of two real vectors */
-static inline float inner_prod(float *x, float *y, int N)
+static inline float inner_prod(spx_word16_t *x, spx_word16_t *y, int N)
{
int i;
float ret=0;
for (i=0;i<N;i++)
- ret += x[i]*y[i];
+ ret += (1.f*x[i])*y[i];
return ret;
}
/** Compute power spectrum of a half-complex (packed) vector */
-static inline void power_spectrum(float *X, float *ps, int N)
+static inline void power_spectrum(spx_word16_t *X, spx_word32_t *ps, int N)
{
int i, j;
- ps[0]=X[0]*X[0];
+ ps[0]=MULT16_16(X[0],X[0]);
for (i=1,j=1;i<N-1;i+=2,j++)
{
- ps[j] = X[i]*X[i] + X[i+1]*X[i+1];
+ ps[j] = MULT16_16(X[i],X[i]) + MULT16_16(X[i+1],X[i+1]);
}
- ps[j]=X[i]*X[i];
+ ps[j]=MULT16_16(X[i],X[i]);
}
/** Compute cross-power spectrum of a half-complex (packed) vectors and add to acc */
-static inline void spectral_mul_accum(float *X, float *Y, float *acc, int N)
+#ifdef FIXED_POINT
+static inline void spectral_mul_accum(spx_word16_t *X, spx_word16_t *Y, spx_word16_t *acc, int N, int M)
{
- int i;
- acc[0] += X[0]*Y[0];
+ int i,j;
+ spx_word32_t tmp1=0,tmp2=0;
+ for (j=0;j<M;j++)
+ {
+ tmp1 = MAC16_16(tmp1, X[j*N],Y[j*N]);
+ }
+ acc[0] = SHR32(tmp1,WEIGHT_SHIFT);
for (i=1;i<N-1;i+=2)
{
- acc[i] += (X[i]*Y[i] - X[i+1]*Y[i+1]);
- acc[i+1] += (X[i+1]*Y[i] + X[i]*Y[i+1]);
+ tmp1 = tmp2 = 0;
+ for (j=0;j<M;j++)
+ {
+ tmp1 = SUB32(MAC16_16(tmp1, X[j*N+i],Y[j*N+i]), MULT16_16(X[j*N+i+1],Y[j*N+i+1]));
+ tmp2 = MAC16_16(MAC16_16(tmp2, X[j*N+i+1],Y[j*N+i]), X[j*N+i], Y[j*N+i+1]);
+ }
+ acc[i] = SHR32(tmp1,WEIGHT_SHIFT);
+ acc[i+1] = SHR32(tmp2,WEIGHT_SHIFT);
}
- acc[i] += X[i]*Y[i];
+ tmp1 = tmp2 = 0;
+ for (j=0;j<M;j++)
+ {
+ tmp1 = MAC16_16(tmp1, X[(j+1)*N-1],Y[(j+1)*N-1]);
+ }
+ acc[N-1] = SHR32(tmp1,WEIGHT_SHIFT);
}
+#else
+static inline void spectral_mul_accum(spx_word16_t *X, spx_word16_t *Y, spx_word16_t *acc, int N, int M)
+{
+ int i,j;
+ for (i=0;i<N;i++)
+ acc[i] = 0;
+ for (j=0;j<M;j++)
+ {
+ acc[0] += X[0]*Y[0];
+ for (i=1;i<N-1;i+=2)
+ {
+ acc[i] += (X[i]*Y[i] - X[i+1]*Y[i+1]);
+ acc[i+1] += (X[i+1]*Y[i] + X[i]*Y[i+1]);
+ }
+ acc[i] += X[i]*Y[i];
+ }
+}
+#endif
/** Compute weighted cross-power spectrum of a half-complex (packed) vector with conjugate */
-static inline void weighted_spectral_mul_conj(float *w, float *X, float *Y, float *prod, int N)
+static inline void weighted_spectral_mul_conj(float *w, spx_word16_t *X, spx_word16_t *Y, spx_word16_t *prod, int N)
{
int i, j;
- prod[0] = w[0]*X[0]*Y[0];
+ prod[0] = w[0]*MULT16_16(X[0],Y[0]);
for (i=1,j=1;i<N-1;i+=2,j++)
{
- prod[i] = w[j]*(X[i]*Y[i] + X[i+1]*Y[i+1]);
- prod[i+1] = w[j]*(-X[i+1]*Y[i] + X[i]*Y[i+1]);
+ prod[i] = w[j]*MAC16_16(MULT16_16(X[i],Y[i]), X[i+1],Y[i+1]);
+ prod[i+1] = w[j]*MAC16_16(MULT16_16(-X[i+1],Y[i]), X[i],Y[i+1]);
}
- prod[i] = w[j]*X[i]*Y[i];
+ prod[i] = w[j]*MULT16_16(X[i],Y[i]);
}
@@ -127,24 +202,24 @@
st->fft_table = spx_fft_init(N);
- st->e = (float*)speex_alloc(N*sizeof(float));
- st->x = (float*)speex_alloc(N*sizeof(float));
- st->d = (float*)speex_alloc(N*sizeof(float));
- st->y = (float*)speex_alloc(N*sizeof(float));
- st->Yps = (float*)speex_alloc(N*sizeof(float));
- st->last_y = (float*)speex_alloc(N*sizeof(float));
- st->Yf = (float*)speex_alloc((st->frame_size+1)*sizeof(float));
- st->Rf = (float*)speex_alloc((st->frame_size+1)*sizeof(float));
- st->Xf = (float*)speex_alloc((st->frame_size+1)*sizeof(float));
- st->Yh = (float*)speex_alloc((st->frame_size+1)*sizeof(float));
- st->Eh = (float*)speex_alloc((st->frame_size+1)*sizeof(float));
+ st->e = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->x = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->d = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->y = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->Yps = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
+ st->last_y = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->Yf = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t));
+ st->Rf = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t));
+ st->Xf = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t));
+ st->Yh = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t));
+ st->Eh = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t));
- st->X = (float*)speex_alloc(M*N*sizeof(float));
- st->Y = (float*)speex_alloc(N*sizeof(float));
- st->E = (float*)speex_alloc(N*sizeof(float));
- st->W = (float*)speex_alloc(M*N*sizeof(float));
- st->PHI = (float*)speex_alloc(M*N*sizeof(float));
- st->power = (float*)speex_alloc((frame_size+1)*sizeof(float));
+ st->X = (spx_word16_t*)speex_alloc(M*N*sizeof(spx_word16_t));
+ st->Y = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->E = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
+ st->W = (spx_word16_t*)speex_alloc(M*N*sizeof(spx_word16_t));
+ st->PHI = (spx_word16_t*)speex_alloc(M*N*sizeof(spx_word16_t));
+ st->power = (spx_word32_t*)speex_alloc((frame_size+1)*sizeof(spx_word32_t));
st->power_1 = (float*)speex_alloc((frame_size+1)*sizeof(float));
for (i=0;i<N*M;i++)
@@ -206,6 +281,7 @@
speex_free(st);
}
+extern int fixed_point;
/** Performs echo cancellation on a frame */
void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out, float *Yout)
{
@@ -238,26 +314,21 @@
st->X[i]=st->X[i+N];
/* Convert x (echo input) to frequency domain */
- spx_fft_float(st->fft_table, st->x, &st->X[(M-1)*N]);
+ spx_fft(st->fft_table, st->x, &st->X[(M-1)*N]);
+
/* Compute filter response Y */
- for (i=0;i<N;i++)
- st->Y[i] = 0;
- for (j=0;j<M;j++)
- spectral_mul_accum(&st->X[j*N], &st->W[j*N], st->Y, N);
+ spectral_mul_accum(st->X, st->W, st->Y, N, M);
- /* Convert Y (filter response) to time domain */
- for (i=0;i<N;i++)
- st->Y[i] *= WEIGHT_SCALING_1;
-
- spx_ifft_float(st->fft_table, st->Y, st->y);
+ spx_ifft(st->fft_table, st->Y, st->y);
/* Compute error signal (signal with echo removed) */
for (i=0;i<st->frame_size;i++)
{
- float tmp_out;
- tmp_out = (float)ref[i] - st->y[i+st->frame_size];
+ spx_word32_t tmp_out;
+ tmp_out = SUB32(EXTEND32(ref[i]), EXTEND32(st->y[i+st->frame_size]));
st->e[i] = 0;
+ /* Do we need saturation? */
st->e[i+st->frame_size] = tmp_out;
/* Saturation */
@@ -273,10 +344,10 @@
Syy = inner_prod(st->y+st->frame_size, st->y+st->frame_size, st->frame_size);
/* Convert error to frequency domain */
- spx_fft_float(st->fft_table, st->e, st->E);
+ spx_fft(st->fft_table, st->e, st->E);
for (i=0;i<st->frame_size;i++)
st->y[i] = 0;
- spx_fft_float(st->fft_table, st->y, st->Y);
+ spx_fft(st->fft_table, st->y, st->Y);
/* Compute power spectrum of echo (X), error (E) and filter response (Y) */
power_spectrum(st->E, st->Rf, N);
@@ -285,20 +356,20 @@
/* Smooth echo energy estimate over time */
for (j=0;j<=st->frame_size;j++)
- st->power[j] = (1-ss)*st->power[j] + ss*st->Xf[j];
+ st->power[j] = (1-ss)*st->power[j] + 1 + ss*st->Xf[j];
{
float Pey = 0, Pyy=0;
float alpha;
for (j=0;j<=st->frame_size;j++)
{
- float E, Y, Eh, Yh;
+ 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*Yh;
- Pyy += Yh*Yh;
+ Pey += Eh*1.f*Yh;
+ Pyy += Yh*1.f*Yh;
st->Eh[j] = .95*Eh - E;
st->Yh[j] = .95*Yh - Y;
}
@@ -341,16 +412,16 @@
{
float r;
/* Compute frequency-domain adaptation mask */
- r = leak_estimate*st->Yf[i] / (1+st->Rf[i]);
+ r = leak_estimate*st->Yf[i] / (1.f+st->Rf[i]);
if (r>1)
r = 1;
st->power_1[i] = WEIGHT_SCALING*adapt_rate*r/(1.f+st->power[i]);
+ /*printf ("%f ", st->power_1[i]);*/
}
} else {
for (i=0;i<=st->frame_size;i++)
- st->power_1[i] = WEIGHT_SCALING*adapt_rate/(1.f+st->power[i]);
+ st->power_1[i] = WEIGHT_SCALING*adapt_rate/(1.f+st->power[i]);
}
-
/* Compute weight gradient */
for (j=0;j<M;j++)
{
@@ -359,24 +430,28 @@
/* Gradient descent */
for (i=0;i<M*N;i++)
+ {
st->W[i] += st->PHI[i];
+ //printf ("%f ", st->PHI[i]);
+ }
+ /*if (st->cancel_count==1100)
+ for (i=0;i<M*N;i++)
+ printf ("%f ", st->W[i]);*/
+ /*printf ("\n");*/
/* AUMDF weight constraint */
for (j=0;j<M;j++)
{
/* Remove the "if" to make this an MDF filter */
+ //if(1)
if (j==M-1 || st->cancel_count%(M-1) == j)
{
float w[N];
+#ifdef FIXED_POINT
float w2[N];
- spx_ifft_float(st->fft_table, &st->W[j*N], w);
-#if 0
- for (i=st->frame_size;i<N;i++)
- {
- w[i]=0;
- }
- spx_fft_float(st->fft_table, w, &st->W[j*N]);
-#else
+ for (i=0;i<N;i++)
+ w2[i] = .03125*st->W[j*N+i];
+ spx_ifft_float(st->fft_table, w2, w);
for (i=0;i<st->frame_size;i++)
{
w[i]=0;
@@ -391,7 +466,17 @@
w2[i]*=.25;
}
for (i=0;i<N;i++)
- st->W[j*N+i] -= w2[i];
+ st->W[j*N+i] -= 32*w2[i];
+#else
+ float w2[N];
+ fixed_point = 0;
+ spx_ifft_float(st->fft_table, &st->W[j*N], w);
+ for (i=st->frame_size;i<N;i++)
+ {
+ w[i]=0;
+ }
+ spx_fft_float(st->fft_table, w, &st->W[j*N]);
+ fixed_point=1;
#endif
}
}
@@ -417,8 +502,8 @@
st->y[i] = (.5-.5*cos(2*M_PI*i/N))*st->last_y[i];
/* Compute power spectrum of the echo */
- spx_fft_float(st->fft_table, st->y, st->Yps);
- power_spectrum(st->Yps, st->Yps, N);
+ spx_fft(st->fft_table, st->y, st->Y);
+ power_spectrum(st->Y, st->Yps, N);
/* Estimate residual echo */
for (i=0;i<=st->frame_size;i++)
More information about the commits
mailing list