From dd007a35405e4908d456782260a3c5583fd5d8b9 Mon Sep 17 00:00:00 2001 From: Jean-Marc Valin Date: Thu, 17 Aug 2006 16:27:15 +1000 Subject: some multi-channel conversion (lots more to do) --- libspeex/mdf.c | 216 +++++++++++++++++++++++++++++++++------------------------ 1 file changed, 125 insertions(+), 91 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index afadb7b..cf9316f 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -111,6 +111,8 @@ struct SpeexEchoState_ { int cancel_count; int adapted; int saturated; + int C; /** Number of input channels (microphones) */ + int K; /** Number of output channels (loudspeakers) */ spx_int32_t sampling_rate; spx_word16_t spec_average; spx_word16_t beta0; @@ -143,7 +145,7 @@ struct SpeexEchoState_ { spx_word16_t *window; spx_word16_t *prop; void *fft_table; - spx_word16_t memX, memD, memE; + spx_word16_t *memX, *memD, *memE; spx_word16_t preemph; spx_word16_t notch_radius; spx_mem_t notch_mem[2]; @@ -153,7 +155,7 @@ struct SpeexEchoState_ { int play_buf_pos; }; -static inline void filter_dc_notch16(const spx_int16_t *in, spx_word16_t radius, spx_word16_t *out, int len, spx_mem_t *mem) +static inline void filter_dc_notch16(const spx_int16_t *in, spx_word16_t radius, spx_word16_t *out, int len, spx_mem_t *mem, int stride) { int i; spx_word16_t den2; @@ -165,7 +167,7 @@ static inline void filter_dc_notch16(const spx_int16_t *in, spx_word16_t radius, /*printf ("%d %d %d %d %d %d\n", num[0], num[1], num[2], den[0], den[1], den[2]);*/ for (i=0;iK = nb_speakers; + st->C = nb_mic; + C=st->C; + K=st->K; st->frame_size = frame_size; st->window_size = 2*frame_size; N = st->window_size; @@ -294,25 +300,25 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) st->fft_table = spx_fft_init(N); - 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->e = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); + st->x = (spx_word16_t*)speex_alloc(K*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 = (spx_word16_t*)speex_alloc((M+1)*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_word32_t*)speex_alloc(M*N*sizeof(spx_word32_t)); + st->y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); + st->Yps = (spx_word32_t*)speex_alloc(C*N*sizeof(spx_word32_t)); + st->last_y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); + st->Yf = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); + st->Rf = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); + st->Xf = (spx_word32_t*)speex_alloc(K*(st->frame_size+1)*sizeof(spx_word32_t)); + st->Yh = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); + st->Eh = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); + + st->X = (spx_word16_t*)speex_alloc(K*(M+1)*N*sizeof(spx_word16_t)); + st->Y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); + st->E = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); + st->W = (spx_word32_t*)speex_alloc(C*K*M*N*sizeof(spx_word32_t)); st->PHI = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); - st->power = (spx_word32_t*)speex_alloc((frame_size+1)*sizeof(spx_word32_t)); - st->power_1 = (spx_float_t*)speex_alloc((frame_size+1)*sizeof(spx_float_t)); + st->power = (spx_word32_t*)speex_alloc(K*(frame_size+1)*sizeof(spx_word32_t)); + st->power_1 = (spx_float_t*)speex_alloc(K*(frame_size+1)*sizeof(spx_float_t)); st->window = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); st->prop = (spx_word16_t*)speex_alloc(M*sizeof(spx_word16_t)); st->wtmp = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); @@ -350,7 +356,9 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) } } - st->memX=st->memD=st->memE=0; + st->memX = (spx_word16_t*)speex_alloc(K*sizeof(spx_word16_t)); + st->memD = (spx_word16_t*)speex_alloc(C*sizeof(spx_word16_t)); + st->memE = (spx_word16_t*)speex_alloc(C*sizeof(spx_word16_t)); st->preemph = QCONST16(.9,15); if (st->sampling_rate<12000) st->notch_radius = QCONST16(.9, 15); @@ -363,7 +371,7 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) st->adapted = 0; st->Pey = st->Pyy = FLOAT_ONE; - st->play_buf = (spx_int16_t*)speex_alloc(2*st->frame_size*sizeof(spx_int16_t)); + st->play_buf = (spx_int16_t*)speex_alloc(K*2*st->frame_size*sizeof(spx_int16_t)); st->play_buf_pos = 0; return st; @@ -395,7 +403,7 @@ void speex_echo_state_reset(SpeexEchoState *st) } /** Destroys an echo canceller state */ -void speex_echo_state_destroy(SpeexEchoState *st) +void mc_echo_state_destroy(SpeexEchoState *st) { spx_fft_destroy(st->fft_table); @@ -428,12 +436,12 @@ void speex_echo_state_destroy(SpeexEchoState *st) speex_free(st); } -void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out, spx_int32_t *Yout) +void mc_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out, spx_int32_t *Yout) { int i; if (st->play_buf_pos>=st->frame_size) { - speex_echo_cancel(st, rec, st->play_buf, out, Yout); + mc_echo_cancel(st, rec, st->play_buf, out, Yout); st->play_buf_pos -= st->frame_size; for (i=0;iframe_size;i++) st->play_buf[i] = st->play_buf[i+st->frame_size]; @@ -449,7 +457,7 @@ void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t } } -void speex_echo_playback(SpeexEchoState *st, const spx_int16_t *play) +void mc_echo_playback(SpeexEchoState *st, const spx_int16_t *play) { if (st->play_buf_pos<=st->frame_size) { @@ -463,10 +471,10 @@ void speex_echo_playback(SpeexEchoState *st, const spx_int16_t *play) } /** Performs echo cancellation on a frame */ -void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_t *echo, spx_int16_t *out, spx_int32_t *Yout) +void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_t *echo, spx_int16_t *out, spx_int32_t *Yout) { - int i,j; - int N,M; + int i,j, chan, speak; + int N,M, C, K; spx_word32_t Syy,See,Sxx; spx_word16_t leak_estimate; spx_word16_t ss, ss_1; @@ -477,6 +485,8 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int N = st->window_size; M = st->M; + C = st->C; + K = st->K; st->cancel_count++; #ifdef FIXED_POINT ss=DIV32_16(11469,M); @@ -486,73 +496,94 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int ss_1 = 1-ss; #endif - filter_dc_notch16(ref, st->notch_radius, st->d, st->frame_size, st->notch_mem); - /* Copy input data to buffer */ - for (i=0;iframe_size;i++) + for (chan = 0; chan < C; chan++) { - spx_word16_t tmp; - spx_word32_t tmp32; - st->x[i] = st->x[i+st->frame_size]; - tmp32 = SUB32(EXTEND32(echo[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX))); -#ifdef FIXED_POINT - /*FIXME: If saturation occurs here, we need to freeze adaptation for M frames (not just one) */ - if (tmp32 > 32767) - { - tmp32 = 32767; - st->saturated = 1; - } - if (tmp32 < -32767) + filter_dc_notch16(ref+chan, st->notch_radius, st->d, st->frame_size, st->notch_mem, C); + /* Copy input data to buffer */ + for (i=0;iframe_size;i++) { - tmp32 = -32767; - st->saturated = 1; - } -#endif - st->x[i+st->frame_size] = EXTRACT16(tmp32); - st->memX = echo[i]; - - tmp = st->d[i]; - st->d[i] = st->d[i+st->frame_size]; - tmp32 = SUB32(EXTEND32(tmp), EXTEND32(MULT16_16_P15(st->preemph, st->memD))); + spx_word16_t tmp; + spx_word32_t tmp32; + tmp = st->d[chan*N+i]; + st->d[chan*N+i] = st->d[chan*N+i+st->frame_size]; + tmp32 = SUB32(EXTEND32(tmp), EXTEND32(MULT16_16_P15(st->preemph, st->memD[chan]))); #ifdef FIXED_POINT - if (tmp32 > 32767) - { - tmp32 = 32767; - st->saturated = 1; - } - if (tmp32 < -32767) - { - tmp32 = -32767; - st->saturated = 1; - } + if (tmp32 > 32767) + { + tmp32 = 32767; + st->saturated = 1; + } + if (tmp32 < -32767) + { + tmp32 = -32767; + st->saturated = 1; + } #endif - st->d[i+st->frame_size] = tmp32; - st->memD = tmp; + st->d[chan*N+i+st->frame_size] = tmp32; + st->memD[chan] = tmp; + } } - /* Shift memory: this could be optimized eventually*/ - for (j=M-1;j>=0;j--) + for (speak = 0; speak < K; speak++) { - for (i=0;iX[(j+1)*N+i] = st->X[j*N+i]; + for (i=0;iframe_size;i++) + { + spx_word16_t tmp; + spx_word32_t tmp32; + st->x[speak+N+i] = st->x[speak+N+i+st->frame_size]; + tmp32 = SUB32(EXTEND32(echo[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX[speak]))); +#ifdef FIXED_POINT + /*FIXME: If saturation occurs here, we need to freeze adaptation for M frames (not just one) */ + if (tmp32 > 32767) + { + tmp32 = 32767; + st->saturated = 1; + } + if (tmp32 < -32767) + { + tmp32 = -32767; + st->saturated = 1; + } +#endif + st->x[speak+N+i+st->frame_size] = EXTRACT16(tmp32); + st->memX[speak] = echo[i]; + } + } + + for (speak = 0; speak < K; speak++) + { + /* Shift memory: this could be optimized eventually*/ + for (j=M-1;j>=0;j--) + { + for (i=0;iX[(j+1)*N*K+speak*N+i] = st->X[j*N*K+speak*N+i]; + } + /* Convert x (echo input) to frequency domain */ + spx_fft(st->fft_table, st->x+speak*N, &st->X[speak*N]); } - - /* Convert x (echo input) to frequency domain */ - spx_fft(st->fft_table, st->x, &st->X[0]); + for (chan = 0; chan < C; chan++) + { #ifdef SMOOTH_BLOCKS - spectral_mul_accum(st->X, st->W, st->Y, N, M); - spx_ifft(st->fft_table, st->Y, st->e); + spectral_mul_accum(st->X, st->W+chan*N*K*M, st->Y+chan*N, N, M*K); + spx_ifft(st->fft_table, st->Y+chan*N, st->e+chan*N); #endif - + } + /* Compute weight gradient */ if (!st->saturated) { - for (j=M-1;j>=0;j--) + for (chan = 0; chan < C; chan++) { - weighted_spectral_mul_conj(st->power_1, &st->X[(j+1)*N], st->E, st->PHI, N); - for (i=0;iW[j*N+i] += MULT16_32_Q15(st->prop[j], st->PHI[i]); - + for (speak = 0; speak < K; speak++) + { + for (j=M-1;j>=0;j--) + { + weighted_spectral_mul_conj(st->power_1+K*N, &st->X[(j+1)*N*K+speak*N], st->E+chan*N, st->PHI, N); + for (i=0;iW[chan*N*K*M + j*N*K + speak*N +i] += MULT16_32_Q15(st->prop[j], st->PHI[i]); + } + } } } @@ -592,9 +623,12 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int } } - /* Compute filter response Y */ - spectral_mul_accum(st->X, st->W, st->Y, N, M); - spx_ifft(st->fft_table, st->Y, st->y); + for (chan = 0; chan < C; chan++) + { + /* Compute filter response Y */ + spectral_mul_accum(st->X, st->W+chan*N*K*M, st->Y+chan*N, N, M*K); + spx_ifft(st->fft_table, st->Y+chan*N, st->y+chan*N); + } /* Compute error signal (for the output with de-emphasis) */ @@ -776,17 +810,17 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int if (Yout) { spx_word16_t leak2; + for (i=0;iframe_size;i++) + st->last_y[i] = st->last_y[st->frame_size+i]; if (st->adapted) { /* If the filter is adapted, take the filtered echo */ - for (i=0;iframe_size;i++) - st->last_y[i] = st->last_y[st->frame_size+i]; for (i=0;iframe_size;i++) st->last_y[st->frame_size+i] = ref[i]-out[i]; } else { /* If filter isn't adapted yet, all we can do is take the echo signal directly */ - for (i=0;ilast_y[i] = st->x[i]; + for (i=0;iframe_size;i++) + st->last_y[st->frame_size+i] = echo[i]; } /* Apply hanning window (should pre-compute it)*/ @@ -815,7 +849,7 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int } -int speex_echo_ctl(SpeexEchoState *st, int request, void *ptr) +int mc_echo_ctl(SpeexEchoState *st, int request, void *ptr) { switch(request) { -- cgit v1.2.3 From f865d74a61e4c3ffcd985c3927054a092013feb7 Mon Sep 17 00:00:00 2001 From: val058 Date: Thu, 17 Aug 2006 07:27:03 +0000 Subject: more MC work git-svn-id: svn://lasagne.centie.net.au/trunk/audio/aec@179 630fbc32-c412-0410-8d9e-966872dfccd7 --- libspeex/mdf.c | 86 +++++++++++++++++++++++++--------------------------------- 1 file changed, 37 insertions(+), 49 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index cf9316f..8f2fcff 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -475,7 +475,6 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ { int i,j, chan, speak; int N,M, C, K; - spx_word32_t Syy,See,Sxx; spx_word16_t leak_estimate; spx_word16_t ss, ss_1; spx_float_t Pey = FLOAT_ONE, Pyy=FLOAT_ONE; @@ -487,6 +486,8 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ M = st->M; C = st->C; K = st->K; + spx_word32_t Syy[C],See[C],Sxx[K]; + st->cancel_count++; #ifdef FIXED_POINT ss=DIV32_16(11469,M); @@ -628,48 +629,50 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ /* Compute filter response Y */ spectral_mul_accum(st->X, st->W+chan*N*K*M, st->Y+chan*N, N, M*K); spx_ifft(st->fft_table, st->Y+chan*N, st->y+chan*N); - } - /* Compute error signal (for the output with de-emphasis) */ - for (i=0;iframe_size;i++) - { - spx_word32_t tmp_out; + /* Compute error signal (for the output with de-emphasis) */ + for (i=0;iframe_size;i++) + { + spx_word32_t tmp_out; #ifdef SMOOTH_BLOCKS - spx_word16_t y = MULT16_16_Q15(st->window[i+st->frame_size],st->e[i+st->frame_size]) + MULT16_16_Q15(st->window[i],st->y[i+st->frame_size]); - tmp_out = SUB32(EXTEND32(st->d[i+st->frame_size]), EXTEND32(y)); + spx_word16_t y = MULT16_16_Q15(st->window[i+st->frame_size],st->e[chan*N+i+st->frame_size]) + MULT16_16_Q15(st->window[i],st->y[chan*N+i+st->frame_size]); + tmp_out = SUB32(EXTEND32(st->d[chan*N+i+st->frame_size]), EXTEND32(y)); #else - tmp_out = SUB32(EXTEND32(st->d[i+st->frame_size]), EXTEND32(st->y[i+st->frame_size])); + tmp_out = SUB32(EXTEND32(st->d[chan*N+i+st->frame_size]), EXTEND32(st->y[chan*N+i+st->frame_size])); #endif - /* Saturation */ - if (tmp_out>32767) - tmp_out = 32767; - else if (tmp_out<-32768) - tmp_out = -32768; - tmp_out = ADD32(tmp_out, EXTEND32(MULT16_16_P15(st->preemph, st->memE))); - /* This is an arbitrary test for saturation */ - if (ref[i] <= -32000 || ref[i] >= 32000) + /* Saturation */ + if (tmp_out>32767) + tmp_out = 32767; + else if (tmp_out<-32768) + tmp_out = -32768; + tmp_out = ADD32(tmp_out, EXTEND32(MULT16_16_P15(st->preemph, st->memE[chan]))); + /* This is an arbitrary test for saturation */ + if (ref[i*C+chan] <= -32000 || ref[i*C+chan] >= 32000) + { + tmp_out = 0; + st->saturated = 1; + } + out[i*C+chan] = (spx_int16_t)tmp_out; + st->memE[chan] = tmp_out; + } + + /* Compute error signal (filter update version) */ + for (i=0;iframe_size;i++) { - tmp_out = 0; - st->saturated = 1; + st->e[chan*N+i] = 0; + st->e[chan*N+i+st->frame_size] = st->d[chan*N+i+st->frame_size] - st->y[chan*N+i+st->frame_size]; } - out[i] = (spx_int16_t)tmp_out; - st->memE = tmp_out; - } - - /* Compute error signal (filter update version) */ - for (i=0;iframe_size;i++) - { - st->e[i] = 0; - st->e[i+st->frame_size] = st->d[i+st->frame_size] - st->y[i+st->frame_size]; + + /* Compute a bunch of correlations */ + See[chan] = mdf_inner_prod(st->e+chan*N+st->frame_size, st->e+chan*N+st->frame_size, st->frame_size); + See[chan] = ADD32(See[chan], SHR32(EXTEND32(10000),6)); + Syy[chan] = mdf_inner_prod(st->y+chan*N+st->frame_size, st->y+chan*N+st->frame_size, st->frame_size); } - - /* Compute a bunch of correlations */ - See = mdf_inner_prod(st->e+st->frame_size, st->e+st->frame_size, st->frame_size); - See = ADD32(See, SHR32(EXTEND32(10000),6)); - Syy = mdf_inner_prod(st->y+st->frame_size, st->y+st->frame_size, st->frame_size); - Sxx = mdf_inner_prod(st->x+st->frame_size, st->x+st->frame_size, st->frame_size); + + for (speak = 0; speak < K; speak++) + Sxx[speak] = mdf_inner_prod(st->x+speak*N+st->frame_size, st->x+speak*N+st->frame_size, st->frame_size); /* Convert error to frequency domain */ spx_fft(st->fft_table, st->e, st->E); @@ -685,21 +688,6 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ /* Smooth echo energy estimate over time */ for (j=0;j<=st->frame_size;j++) st->power[j] = MULT16_32_Q15(ss_1,st->power[j]) + 1 + MULT16_32_Q15(ss,st->Xf[j]); - - /* Enable this to compute the power based only on the tail (would need to compute more - efficiently to make this really useful */ - if (0) - { - float scale2 = .5f/M; - for (j=0;j<=st->frame_size;j++) - st->power[j] = 0; - for (i=0;iX[i*N], st->Xf, N); - for (j=0;j<=st->frame_size;j++) - st->power[j] += scale2*st->Xf[j]; - } - } /* Compute filtered spectra and (cross-)correlations */ for (j=st->frame_size;j>=0;j--) -- cgit v1.2.3 From d8944185c1f8d479417e3b2e9a2689b76729b764 Mon Sep 17 00:00:00 2001 From: val058 Date: Fri, 18 Aug 2006 01:08:28 +0000 Subject: MC work -- continued git-svn-id: svn://lasagne.centie.net.au/trunk/audio/aec@182 630fbc32-c412-0410-8d9e-966872dfccd7 --- libspeex/mdf.c | 48 +++++++++++++++++++++++++++--------------------- libspeex/testecho.c | 8 ++++---- 2 files changed, 31 insertions(+), 25 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index 8f2fcff..9d5b413 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -306,19 +306,19 @@ SpeexEchoState *mc_echo_state_init(int frame_size, int filter_length, int nb_mic st->y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); st->Yps = (spx_word32_t*)speex_alloc(C*N*sizeof(spx_word32_t)); st->last_y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); - st->Yf = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); - st->Rf = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); - st->Xf = (spx_word32_t*)speex_alloc(K*(st->frame_size+1)*sizeof(spx_word32_t)); - st->Yh = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_t)); - st->Eh = (spx_word32_t*)speex_alloc(C*(st->frame_size+1)*sizeof(spx_word32_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 = (spx_word16_t*)speex_alloc(K*(M+1)*N*sizeof(spx_word16_t)); st->Y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); st->E = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); st->W = (spx_word32_t*)speex_alloc(C*K*M*N*sizeof(spx_word32_t)); st->PHI = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); - st->power = (spx_word32_t*)speex_alloc(K*(frame_size+1)*sizeof(spx_word32_t)); - st->power_1 = (spx_float_t*)speex_alloc(K*(frame_size+1)*sizeof(spx_float_t)); + st->power = (spx_word32_t*)speex_alloc((frame_size+1)*sizeof(spx_word32_t)); + st->power_1 = (spx_float_t*)speex_alloc((frame_size+1)*sizeof(spx_float_t)); st->window = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); st->prop = (spx_word16_t*)speex_alloc(M*sizeof(spx_word16_t)); st->wtmp = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); @@ -486,7 +486,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ M = st->M; C = st->C; K = st->K; - spx_word32_t Syy[C],See[C],Sxx[K]; + spx_word32_t Syy,See,Sxx; st->cancel_count++; #ifdef FIXED_POINT @@ -580,9 +580,9 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ { for (j=M-1;j>=0;j--) { - weighted_spectral_mul_conj(st->power_1+K*N, &st->X[(j+1)*N*K+speak*N], st->E+chan*N, st->PHI, N); + weighted_spectral_mul_conj(st->power_1, &st->X[(j+1)*N*K+speak*N], st->E+chan*N, st->PHI, N); for (i=0;iW[chan*N*K*M + j*N*K + speak*N +i] += MULT16_32_Q15(st->prop[j], st->PHI[i]); + st->W[chan*N*K*M + j*N*K + speak*N + i] += MULT16_32_Q15(st->prop[j], st->PHI[i]); } } } @@ -590,6 +590,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ st->saturated = 0; + /* FIXME: MC conversion required */ /* Update weight to prevent circular convolution (MDF / AUMDF) */ for (j=0;jx+speak*N+st->frame_size, st->x+speak*N+st->frame_size, st->frame_size); + + for (chan = 0; chan < C; chan++) { /* Compute filter response Y */ @@ -666,25 +671,25 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ } /* Compute a bunch of correlations */ - See[chan] = mdf_inner_prod(st->e+chan*N+st->frame_size, st->e+chan*N+st->frame_size, st->frame_size); - See[chan] = ADD32(See[chan], SHR32(EXTEND32(10000),6)); - Syy[chan] = mdf_inner_prod(st->y+chan*N+st->frame_size, st->y+chan*N+st->frame_size, st->frame_size); + See += mdf_inner_prod(st->e+chan*N+st->frame_size, st->e+chan*N+st->frame_size, st->frame_size); + Syy += mdf_inner_prod(st->y+chan*N+st->frame_size, st->y+chan*N+st->frame_size, st->frame_size); + + /* Convert error to frequency domain */ + spx_fft(st->fft_table, st->e+chan*N, st->E+chan*N); + for (i=0;iframe_size;i++) + st->y[i+chan*N] = 0; + spx_fft(st->fft_table, st->y+chan*N, st->Y+chan*N); } + See = ADD32(See, SHR32(EXTEND32(10000),6)); - for (speak = 0; speak < K; speak++) - Sxx[speak] = mdf_inner_prod(st->x+speak*N+st->frame_size, st->x+speak*N+st->frame_size, st->frame_size); - - /* Convert error to frequency domain */ - spx_fft(st->fft_table, st->e, st->E); - for (i=0;iframe_size;i++) - st->y[i] = 0; - spx_fft(st->fft_table, st->y, st->Y); + /* FIXME: MC conversion required */ /* Compute power spectrum of echo (X), error (E) and filter response (Y) */ power_spectrum(st->E, st->Rf, N); power_spectrum(st->Y, st->Yf, N); power_spectrum(st->X, st->Xf, N); + /* FIXME: MC conversion required */ /* Smooth echo energy estimate over time */ for (j=0;j<=st->frame_size;j++) st->power[j] = MULT16_32_Q15(ss_1,st->power[j]) + 1 + MULT16_32_Q15(ss,st->Xf[j]); @@ -794,6 +799,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ st->sum_adapt = ADD32(st->sum_adapt,adapt_rate); } + /* FIXME: MC conversion required */ /* Compute spectrum of estimated echo for use in an echo post-filter (if necessary)*/ if (Yout) { diff --git a/libspeex/testecho.c b/libspeex/testecho.c index f78499e..fc5bf34 100644 --- a/libspeex/testecho.c +++ b/libspeex/testecho.c @@ -32,19 +32,19 @@ int main(int argc, char **argv) ref_fd = open (argv[1], O_RDONLY); e_fd = open (argv[3], O_WRONLY | O_CREAT | O_TRUNC, 0644); - st = speex_echo_state_init(NN, TAIL); + st = mc_echo_state_init(NN, TAIL, 1, 1); den = speex_preprocess_state_init(NN, 8000); int tmp = 8000; - speex_echo_ctl(st, SPEEX_ECHO_SET_SAMPLING_RATE, &tmp); + mc_echo_ctl(st, SPEEX_ECHO_SET_SAMPLING_RATE, &tmp); while (read(ref_fd, ref_buf, NN*2)) { read(echo_fd, echo_buf, NN*2); - speex_echo_cancel(st, ref_buf, echo_buf, e_buf, noise); + mc_echo_cancel(st, ref_buf, echo_buf, e_buf, noise); /*speex_preprocess(den, e_buf, noise);*/ write(e_fd, e_buf, NN*2); } - speex_echo_state_destroy(st); + mc_echo_state_destroy(st); speex_preprocess_state_destroy(den); close(e_fd); close(echo_fd); -- cgit v1.2.3 From 8bc9da44c54efe8570662a40668965b83b423b3d Mon Sep 17 00:00:00 2001 From: val058 Date: Fri, 18 Aug 2006 01:22:15 +0000 Subject: seems to at least be working for single channel git-svn-id: svn://lasagne.centie.net.au/trunk/audio/aec@183 630fbc32-c412-0410-8d9e-966872dfccd7 --- libspeex/mdf.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index 9d5b413..cf9354a 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -486,7 +486,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ M = st->M; C = st->C; K = st->K; - spx_word32_t Syy,See,Sxx; + spx_word32_t Syy=0,See=0,Sxx=0; st->cancel_count++; #ifdef FIXED_POINT @@ -531,7 +531,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ { spx_word16_t tmp; spx_word32_t tmp32; - st->x[speak+N+i] = st->x[speak+N+i+st->frame_size]; + st->x[speak*N+i] = st->x[speak*N+i+st->frame_size]; tmp32 = SUB32(EXTEND32(echo[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX[speak]))); #ifdef FIXED_POINT /*FIXME: If saturation occurs here, we need to freeze adaptation for M frames (not just one) */ @@ -546,7 +546,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ st->saturated = 1; } #endif - st->x[speak+N+i+st->frame_size] = EXTRACT16(tmp32); + st->x[speak*N+i+st->frame_size] = EXTRACT16(tmp32); st->memX[speak] = echo[i]; } } -- cgit v1.2.3 From 4f1285784657fd0f045948ad7d873244e2c4f07c Mon Sep 17 00:00:00 2001 From: val058 Date: Fri, 18 Aug 2006 01:45:53 +0000 Subject: most of conversion done (but MC untested so far) git-svn-id: svn://lasagne.centie.net.au/trunk/audio/aec@184 630fbc32-c412-0410-8d9e-966872dfccd7 --- libspeex/mdf.c | 40 +++++++++++++++++++++++++++------------- 1 file changed, 27 insertions(+), 13 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index cf9354a..8b5d90d 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -206,6 +206,18 @@ static inline void power_spectrum(const spx_word16_t *X, spx_word32_t *ps, int N ps[j]=MULT16_16(X[i],X[i]); } +/** Compute power spectrum of a half-complex (packed) vector and accumulate */ +static inline void power_spectrum_accum(const spx_word16_t *X, spx_word32_t *ps, int N) +{ + int i, j; + ps[0]+=MULT16_16(X[0],X[0]); + for (i=1,j=1;ix[speak*N+i] = st->x[speak*N+i+st->frame_size]; - tmp32 = SUB32(EXTEND32(echo[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX[speak]))); + tmp32 = SUB32(EXTEND32(echo[i*K+speak]), EXTEND32(MULT16_16_P15(st->preemph, st->memX[speak]))); #ifdef FIXED_POINT /*FIXME: If saturation occurs here, we need to freeze adaptation for M frames (not just one) */ if (tmp32 > 32767) @@ -547,7 +559,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ } #endif st->x[speak*N+i+st->frame_size] = EXTRACT16(tmp32); - st->memX[speak] = echo[i]; + st->memX[speak] = echo[i*K+speak]; } } @@ -625,10 +637,10 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ } } - for (speak = 0; speak < K; speak++) - Sxx += mdf_inner_prod(st->x+speak*N+st->frame_size, st->x+speak*N+st->frame_size, st->frame_size); - - + /* So we can use power_spectrum_accum */ + for (i=0;i<=st->frame_size;i++) + st->Rf[i] = st->Yf[i] = st->Xf[i] = 0; + for (chan = 0; chan < C; chan++) { /* Compute filter response Y */ @@ -679,17 +691,19 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ for (i=0;iframe_size;i++) st->y[i+chan*N] = 0; spx_fft(st->fft_table, st->y+chan*N, st->Y+chan*N); + + /* Compute power spectrum of echo (X), error (E) and filter response (Y) */ + power_spectrum_accum(st->E+chan*N, st->Rf, N); + power_spectrum_accum(st->Y+chan*N, st->Yf, N); } See = ADD32(See, SHR32(EXTEND32(10000),6)); - - /* FIXME: MC conversion required */ - /* Compute power spectrum of echo (X), error (E) and filter response (Y) */ - power_spectrum(st->E, st->Rf, N); - power_spectrum(st->Y, st->Yf, N); - power_spectrum(st->X, st->Xf, N); + for (speak = 0; speak < K; speak++) + { + Sxx += mdf_inner_prod(st->x+speak*N+st->frame_size, st->x+speak*N+st->frame_size, st->frame_size); + power_spectrum_accum(st->X+speak*N, st->Xf, N); + } - /* FIXME: MC conversion required */ /* Smooth echo energy estimate over time */ for (j=0;j<=st->frame_size;j++) st->power[j] = MULT16_32_Q15(ss_1,st->power[j]) + 1 + MULT16_32_Q15(ss,st->Xf[j]); -- cgit v1.2.3 From 0f8e7a7dc2a54dd005a92c0b8d5ce716379ffcd4 Mon Sep 17 00:00:00 2001 From: val058 Date: Fri, 18 Aug 2006 01:54:30 +0000 Subject: 2-input doesn't crash anymore (1st chan is good, 2nd needs fixing) git-svn-id: svn://lasagne.centie.net.au/trunk/audio/aec@185 630fbc32-c412-0410-8d9e-966872dfccd7 --- libspeex/mdf.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index 8b5d90d..ede56ab 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -314,7 +314,7 @@ SpeexEchoState *mc_echo_state_init(int frame_size, int filter_length, int nb_mic st->e = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); st->x = (spx_word16_t*)speex_alloc(K*N*sizeof(spx_word16_t)); - st->d = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); + st->d = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); st->y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); st->Yps = (spx_word32_t*)speex_alloc(C*N*sizeof(spx_word32_t)); st->last_y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t)); @@ -511,7 +511,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ for (chan = 0; chan < C; chan++) { - filter_dc_notch16(ref+chan, st->notch_radius, st->d, st->frame_size, st->notch_mem, C); + filter_dc_notch16(ref+chan, st->notch_radius, st->d+chan*N, st->frame_size, st->notch_mem, C); /* Copy input data to buffer */ for (i=0;iframe_size;i++) { -- cgit v1.2.3 From 6bd022014a21ecca9c27d6041397009a5933ac39 Mon Sep 17 00:00:00 2001 From: val058 Date: Fri, 18 Aug 2006 02:15:25 +0000 Subject: Fixed weight update and notch filter memory. Now works for multiple mics (multiple speakers untested) git-svn-id: svn://lasagne.centie.net.au/trunk/audio/aec@186 630fbc32-c412-0410-8d9e-966872dfccd7 --- libspeex/mdf.c | 73 +++++++++++++++++++++++++++++++++------------------------- 1 file changed, 41 insertions(+), 32 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index ede56ab..3d79383 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -148,7 +148,7 @@ struct SpeexEchoState_ { spx_word16_t *memX, *memD, *memE; spx_word16_t preemph; spx_word16_t notch_radius; - spx_mem_t notch_mem[2]; + spx_mem_t *notch_mem; /* NOTE: If you only use speex_echo_cancel() and want to save some memory, remove this */ spx_int16_t *play_buf; @@ -347,7 +347,7 @@ SpeexEchoState *mc_echo_state_init(int frame_size, int filter_length, int nb_mic #endif for (i=0;i<=st->frame_size;i++) st->power_1[i] = FLOAT_ONE; - for (i=0;iW[i] = 0; for (i=0;iPHI[i] = 0; @@ -379,7 +379,7 @@ SpeexEchoState *mc_echo_state_init(int frame_size, int filter_length, int nb_mic else st->notch_radius = QCONST16(.992, 15); - st->notch_mem[0] = st->notch_mem[1] = 0; + st->notch_mem = (spx_mem_t*)speex_alloc(2*C*sizeof(spx_mem_t)); st->adapted = 0; st->Pey = st->Pyy = FLOAT_ONE; @@ -392,10 +392,12 @@ SpeexEchoState *mc_echo_state_init(int frame_size, int filter_length, int nb_mic /** Resets echo canceller state */ void speex_echo_state_reset(SpeexEchoState *st) { - int i, M, N; + int i, M, N, C, K; st->cancel_count=0; N = st->window_size; M = st->M; + C=st->C; + K=st->K; for (i=0;iW[i] = 0; for (i=0;ipower[i] = 0; for (i=0;iE[i] = 0; - st->notch_mem[0] = st->notch_mem[1] = 0; + for (i=0;i<2*C;i++) + st->notch_mem[i] = 0; st->saturated = 0; st->adapted = 0; @@ -511,7 +514,7 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ for (chan = 0; chan < C; chan++) { - filter_dc_notch16(ref+chan, st->notch_radius, st->d+chan*N, st->frame_size, st->notch_mem, C); + filter_dc_notch16(ref+chan, st->notch_radius, st->d+chan*N, st->frame_size, st->notch_mem+2*chan, C); /* Copy input data to buffer */ for (i=0;iframe_size;i++) { @@ -604,39 +607,45 @@ void mc_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_ /* FIXME: MC conversion required */ /* Update weight to prevent circular convolution (MDF / AUMDF) */ - for (j=0;jcancel_count%(M-1) == j-1) + for (speak = 0; speak < K; speak++) { -#ifdef FIXED_POINT - for (i=0;iwtmp2[i] = EXTRACT16(PSHR32(st->W[j*N+i],NORMALIZE_SCALEDOWN+16)); - spx_ifft(st->fft_table, st->wtmp2, st->wtmp); - for (i=0;iframe_size;i++) + for (j=0;jwtmp[i]=0; - } - for (i=st->frame_size;iwtmp[i]=SHL16(st->wtmp[i],NORMALIZE_SCALEUP); - } - spx_fft(st->fft_table, st->wtmp, st->wtmp2); - /* The "-1" in the shift is a sort of kludge that trades less efficient update speed for decrease noise */ - for (i=0;iW[j*N+i] -= SHL32(EXTEND32(st->wtmp2[i]),16+NORMALIZE_SCALEDOWN-NORMALIZE_SCALEUP-1); + /* This is a variant of the Alternatively Updated MDF (AUMDF) */ + /* Remove the "if" to make this an MDF filter */ + if (j==0 || st->cancel_count%(M-1) == j-1) + { +#ifdef FIXED_POINT + for (i=0;iwtmp2[i] = EXTRACT16(PSHR32(st->W[chan*N*K*M + j*N*K + speak*N + i],NORMALIZE_SCALEDOWN+16)); + spx_ifft(st->fft_table, st->wtmp2, st->wtmp); + for (i=0;iframe_size;i++) + { + st->wtmp[i]=0; + } + for (i=st->frame_size;iwtmp[i]=SHL16(st->wtmp[i],NORMALIZE_SCALEUP); + } + spx_fft(st->fft_table, st->wtmp, st->wtmp2); + /* The "-1" in the shift is a sort of kludge that trades less efficient update speed for decrease noise */ + for (i=0;iW[chan*N*K*M + j*N*K + speak*N + i] -= SHL32(EXTEND32(st->wtmp2[i]),16+NORMALIZE_SCALEDOWN-NORMALIZE_SCALEUP-1); #else - spx_ifft(st->fft_table, &st->W[j*N], st->wtmp); - for (i=st->frame_size;iwtmp[i]=0; - } - spx_fft(st->fft_table, st->wtmp, &st->W[j*N]); + spx_ifft(st->fft_table, &st->W[chan*N*K*M + j*N*K + speak*N], st->wtmp); + for (i=st->frame_size;iwtmp[i]=0; + } + spx_fft(st->fft_table, st->wtmp, &st->W[chan*N*K*M + j*N*K + speak*N]); #endif + } + } } } - + /* So we can use power_spectrum_accum */ for (i=0;i<=st->frame_size;i++) st->Rf[i] = st->Yf[i] = st->Xf[i] = 0; -- cgit v1.2.3 From d794a6ab626ab6d3bbe6689185b22a1df993b8e0 Mon Sep 17 00:00:00 2001 From: Jean-Marc Valin Date: Fri, 30 May 2008 21:22:28 +1000 Subject: Made multi-channel AEC API compatible with the previous one. --- libspeex/mdf.c | 7 ++++++- libspeex/testecho.c | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) (limited to 'libspeex') diff --git a/libspeex/mdf.c b/libspeex/mdf.c index b9945db..d8bacaf 100644 --- a/libspeex/mdf.c +++ b/libspeex/mdf.c @@ -399,7 +399,12 @@ static void dump_audio(const spx_int16_t *rec, const spx_int16_t *play, const sp #endif /** Creates a new echo canceller state */ -EXPORT SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length, int nb_mic, int nb_speakers) +EXPORT SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) +{ + return speex_echo_state_init_mc(frame_size, filter_length, 1, 1); +} + +EXPORT SpeexEchoState *speex_echo_state_init_mc(int frame_size, int filter_length, int nb_mic, int nb_speakers) { int i,N,M, C, K; SpeexEchoState *st = (SpeexEchoState *)speex_alloc(sizeof(SpeexEchoState)); diff --git a/libspeex/testecho.c b/libspeex/testecho.c index 01b4539..5ae855f 100644 --- a/libspeex/testecho.c +++ b/libspeex/testecho.c @@ -31,7 +31,7 @@ int main(int argc, char **argv) ref_fd = fopen(argv[1], "rb"); e_fd = fopen(argv[3], "wb"); - st = speex_echo_state_init(NN, TAIL, 1, 1); + st = speex_echo_state_init(NN, TAIL); den = speex_preprocess_state_init(NN, sampleRate); speex_echo_ctl(st, SPEEX_ECHO_SET_SAMPLING_RATE, &sampleRate); speex_preprocess_ctl(den, SPEEX_PREPROCESS_SET_ECHO_STATE, st); -- cgit v1.2.3