Added a wrapper around the FFT so any FFT can be used
authorjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Fri, 2 Dec 2005 10:08:44 +0000 (10:08 +0000)
committerjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Fri, 2 Dec 2005 10:08:44 +0000 (10:08 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@10523 0101bb08-14d6-0310-b084-bc0e0c8e3800

include/speex/speex_echo.h
libspeex/Makefile.am
libspeex/fftwrap.c [new file with mode: 0644]
libspeex/fftwrap.h [new file with mode: 0644]
libspeex/mdf.c

index 6f1af75..0ade2bd 100644 (file)
@@ -38,7 +38,7 @@
 extern "C" {
 #endif
 
 extern "C" {
 #endif
 
-struct drft_lookup;
+/*struct drft_lookup;*/
 
 /** Speex echo cancellation state. */
 typedef struct SpeexEchoState {
 
 /** Speex echo cancellation state. */
 typedef struct SpeexEchoState {
@@ -48,7 +48,7 @@ typedef struct SpeexEchoState {
    int cancel_count;
    int adapted;
    float sum_adapt;
    int cancel_count;
    int adapted;
    float sum_adapt;
-   
+   float *e;
    float *x;
    float *X;
    float *d;
    float *x;
    float *X;
    float *d;
@@ -68,7 +68,8 @@ typedef struct SpeexEchoState {
    float *Yh;
    float Pey;
    float Pyy;
    float *Yh;
    float Pey;
    float Pyy;
-   struct drft_lookup *fft_lookup;
+   /*struct drft_lookup *fft_lookup;*/
+   void *fft_table;
 
 
 } SpeexEchoState;
 
 
 } SpeexEchoState;
index da29c9e..da3b29d 100644 (file)
@@ -10,19 +10,20 @@ lib_LTLIBRARIES = libspeex.la
 
 # Sources for compilation in the library
 libspeex_la_SOURCES = nb_celp.c        sb_celp.c       lpc.c   ltp.c   lsp.c   quant_lsp.c \
 
 # Sources for compilation in the library
 libspeex_la_SOURCES = nb_celp.c        sb_celp.c       lpc.c   ltp.c   lsp.c   quant_lsp.c \
-               lsp_tables_nb.c         gain_table.c    gain_table_lbr.c        cb_search.c     filters.c       bits.c \
-               modes.c         speex.c         vq.c    high_lsp_tables.c       vbr.c   hexc_table.c \
-               exc_5_256_table.c       exc_5_64_table.c        exc_8_128_table.c       exc_10_32_table.c \
-               exc_10_16_table.c       exc_20_32_table.c       hexc_10_32_table.c      misc.c  speex_header.c \
-               speex_callbacks.c       math_approx.c   stereo.c        preprocess.c    smallft.c       lbr_48k_tables.c \
-               jitter.c        mdf.c vorbis_psy.c
+                       lsp_tables_nb.c         gain_table.c    gain_table_lbr.c        cb_search.c     filters.c       bits.c \
+                       modes.c         speex.c         vq.c    high_lsp_tables.c       vbr.c   hexc_table.c \
+                       exc_5_256_table.c       exc_5_64_table.c        exc_8_128_table.c       exc_10_32_table.c \
+                       exc_10_16_table.c       exc_20_32_table.c       hexc_10_32_table.c      misc.c  speex_header.c \
+                       speex_callbacks.c       math_approx.c   stereo.c        preprocess.c    smallft.c       lbr_48k_tables.c \
+                       jitter.c        mdf.c vorbis_psy.c fftwrap.c
 
 noinst_HEADERS = lsp.h         nb_celp.h       lpc.h   lpc_bfin.h      ltp.h   quant_lsp.h \
 
 noinst_HEADERS = lsp.h         nb_celp.h       lpc.h   lpc_bfin.h      ltp.h   quant_lsp.h \
-               cb_search.h     filters.h       stack_alloc.h   vq.h    vq_sse.h        vq_arm4.h       vq_bfin.h \
-               modes.h         sb_celp.h       vbr.h   misc.h  misc_bfin.h     ltp_sse.h       ltp_arm4.h \
-               ltp_bfin.h      filters_sse.h   filters_arm4.h  filters_bfin.h  math_approx.h \
-               smallft.h       arch.h  fixed_arm4.h    fixed_arm5e.h   fixed_bfin.h    fixed_debug.h \
-               fixed_generic.h         cb_search_sse.h         cb_search_arm4.h        cb_search_bfin.h vorbis_psy.h
+                       cb_search.h     filters.h       stack_alloc.h   vq.h    vq_sse.h        vq_arm4.h       vq_bfin.h \
+                       modes.h         sb_celp.h       vbr.h   misc.h  misc_bfin.h     ltp_sse.h       ltp_arm4.h \
+                       ltp_bfin.h      filters_sse.h   filters_arm4.h  filters_bfin.h  math_approx.h \
+                       smallft.h       arch.h  fixed_arm4.h    fixed_arm5e.h   fixed_bfin.h    fixed_debug.h \
+                       fixed_generic.h         cb_search_sse.h         cb_search_arm4.h        cb_search_bfin.h vorbis_psy.h \
+       fftwrap.h
 
 
 libspeex_la_LDFLAGS = -version-info @SPEEX_LT_CURRENT@:@SPEEX_LT_REVISION@:@SPEEX_LT_AGE@
 
 
 libspeex_la_LDFLAGS = -version-info @SPEEX_LT_CURRENT@:@SPEEX_LT_REVISION@:@SPEEX_LT_AGE@
diff --git a/libspeex/fftwrap.c b/libspeex/fftwrap.c
new file mode 100644 (file)
index 0000000..53becb6
--- /dev/null
@@ -0,0 +1,107 @@
+/* Copyright (C) 2005 Jean-Marc Valin 
+   File: fftwrap.c
+
+   Wrapper for various FFTs 
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+   
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+   
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+   
+   - Neither the name of the Xiph.org Foundation nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+   
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#define USE_SMALLFT
+
+
+#include "misc.h"
+
+
+#ifdef USE_SMALLFT
+
+#include "smallft.h"
+#include <math.h>
+
+void *spx_fft_init(int size)
+{
+   struct drft_lookup *table;
+   table = speex_alloc(sizeof(struct drft_lookup));
+   spx_drft_init((struct drft_lookup *)table, size);
+   return (void*)table;
+}
+
+void spx_fft_destroy(void *table)
+{
+   spx_drft_clear(table);
+   speex_free(table);
+}
+
+void spx_fft(void *table, float *in, float *out)
+{
+   if (in==out)
+   {
+      int i;
+      speex_warning("FFT should not be done in-place");
+      float scale = 1./((struct drft_lookup *)table)->n;
+      for (i=0;i<((struct drft_lookup *)table)->n;i++)
+         out[i] = scale*in[i];
+   } else {
+      int i;
+      float scale = 1./((struct drft_lookup *)table)->n;
+      for (i=0;i<((struct drft_lookup *)table)->n;i++)
+         out[i] = scale*in[i];
+   }
+   spx_drft_forward((struct drft_lookup *)table, out);
+}
+
+void spx_ifft(void *table, float *in, float *out)
+{
+   if (in==out)
+   {
+      int i;
+      speex_warning("FFT should not be done in-place");
+   } else {
+      int i;
+      for (i=0;i<((struct drft_lookup *)table)->n;i++)
+         out[i] = in[i];
+   }
+   spx_drft_backward((struct drft_lookup *)table, out);
+}
+
+#else
+
+#error No other FFT implemented
+
+#endif
+
+
+
+void spx_fft_float(void *table, float *in, float *out)
+{
+   spx_fft(table, in, out);
+}
+void spx_ifft_float(void *table, float *in, float *out)
+{
+   spx_ifft(table, in, out);
+}
diff --git a/libspeex/fftwrap.h b/libspeex/fftwrap.h
new file mode 100644 (file)
index 0000000..826b38e
--- /dev/null
@@ -0,0 +1,58 @@
+/* Copyright (C) 2005 Jean-Marc Valin 
+   File: fftwrap.h
+
+   Wrapper for various FFTs 
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+   
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+   
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+   
+   - Neither the name of the Xiph.org Foundation nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+   
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#ifndef FFTWRAP_H
+#define FFTWRAP_H
+
+#include "misc.h"
+
+/** Compute tables for an FFT */
+void *spx_fft_init(int size);
+
+/** Destroy tables for an FFT */
+void spx_fft_destroy(void *table);
+
+/** Forward (real to half-complex) transform */
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out);
+
+/** Backward (half-complex to real) transform */
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out);
+
+/** Forward (real to half-complex) transform of float data */
+void spx_fft_float(void *table, float *in, float *out);
+
+/** Backward (half-complex to real) transform of float data */
+void spx_ifft_float(void *table, float *in, float *out);
+
+#endif
index 278bef7..d0ac3da 100644 (file)
@@ -41,6 +41,8 @@
 #include "misc.h"
 #include "speex/speex_echo.h"
 #include "smallft.h"
 #include "misc.h"
 #include "speex/speex_echo.h"
 #include "smallft.h"
+#include "fftwrap.h"
+
 #include <math.h>
 
 #ifndef M_PI
 #include <math.h>
 
 #ifndef M_PI
@@ -114,10 +116,10 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length)
    M = st->M = (filter_length+st->frame_size-1)/frame_size;
    st->cancel_count=0;
    st->sum_adapt = 0;
    M = st->M = (filter_length+st->frame_size-1)/frame_size;
    st->cancel_count=0;
    st->sum_adapt = 0;
-         
-   st->fft_lookup = (struct drft_lookup*)speex_alloc(sizeof(struct drft_lookup));
-   spx_drft_init(st->fft_lookup, N);
+
+   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->x = (float*)speex_alloc(N*sizeof(float));
    st->d = (float*)speex_alloc(N*sizeof(float));
    st->y = (float*)speex_alloc(N*sizeof(float));
@@ -171,8 +173,9 @@ void speex_echo_state_reset(SpeexEchoState *st)
 /** Destroys an echo canceller state */
 void speex_echo_state_destroy(SpeexEchoState *st)
 {
 /** Destroys an echo canceller state */
 void speex_echo_state_destroy(SpeexEchoState *st)
 {
-   spx_drft_clear(st->fft_lookup);
-   speex_free(st->fft_lookup);
+   spx_fft_destroy(st->fft_table);
+
+   speex_free(st->e);
    speex_free(st->x);
    speex_free(st->d);
    speex_free(st->y);
    speex_free(st->x);
    speex_free(st->d);
    speex_free(st->y);
@@ -195,21 +198,18 @@ void speex_echo_state_destroy(SpeexEchoState *st)
    speex_free(st);
 }
 
    speex_free(st);
 }
 
-
 /** Performs echo cancellation on a frame */
 void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out, float *Yout)
 {
    int i,j;
    int N,M;
 /** Performs echo cancellation on a frame */
 void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out, float *Yout)
 {
    int i,j;
    int N,M;
-   float scale;
    float Syy=0,See=0;
    float leak_estimate;
    float ss;
    float adapt_rate;
    float Syy=0,See=0;
    float leak_estimate;
    float ss;
    float adapt_rate;
-
+   
    N = st->window_size;
    M = st->M;
    N = st->window_size;
    M = st->M;
-   scale = 1.0f/N;
    st->cancel_count++;
    ss = 1.0f/st->cancel_count;
    if (ss < .4/M)
    st->cancel_count++;
    ss = 1.0f/st->cancel_count;
    if (ss < .4/M)
@@ -229,25 +229,15 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
    for (i=0;i<N*(M-1);i++)
       st->X[i]=st->X[i+N];
 
    for (i=0;i<N*(M-1);i++)
       st->X[i]=st->X[i+N];
 
-   /* Copy new echo frame */
-   for (i=0;i<N;i++)
-      st->X[(M-1)*N+i]=st->x[i];
-
    /* Convert x (echo input) to frequency domain */
    /* Convert x (echo input) to frequency domain */
-   spx_drft_forward(st->fft_lookup, &st->X[(M-1)*N]);
-
+   spx_fft_float(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);
    /* 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);
-   
-   /* Convert Y (filter response) to time domain */
-   for (i=0;i<N;i++)
-      st->y[i] = st->Y[i];
-   spx_drft_backward(st->fft_lookup, st->y);
-   for (i=0;i<N;i++)
-      st->y[i] *= scale;
+      
+   spx_ifft_float(st->fft_table, st->Y, st->y);
 
    /* Compute error signal (signal with echo removed) */ 
    for (i=0;i<st->frame_size;i++)
 
    /* Compute error signal (signal with echo removed) */ 
    for (i=0;i<st->frame_size;i++)
@@ -255,8 +245,8 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
       float tmp_out;
       tmp_out = (float)ref[i] - st->y[i+st->frame_size];
       
       float tmp_out;
       tmp_out = (float)ref[i] - st->y[i+st->frame_size];
       
-      st->E[i] = 0;
-      st->E[i+st->frame_size] = tmp_out;
+      st->e[i] = 0;
+      st->e[i+st->frame_size] = tmp_out;
       
       /* Saturation */
       if (tmp_out>32767)
       
       /* Saturation */
       if (tmp_out>32767)
@@ -267,17 +257,15 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
    }
    
    /* Compute a bunch of correlations */
    }
    
    /* Compute a bunch of correlations */
-   See = inner_prod(st->E+st->frame_size, st->E+st->frame_size, st->frame_size);
+   See = inner_prod(st->e+st->frame_size, st->e+st->frame_size, st->frame_size);
    Syy = inner_prod(st->y+st->frame_size, st->y+st->frame_size, st->frame_size);
    
    /* Convert error to frequency domain */
    Syy = inner_prod(st->y+st->frame_size, st->y+st->frame_size, st->frame_size);
    
    /* Convert error to frequency domain */
-   spx_drft_forward(st->fft_lookup, st->E);
+   spx_fft_float(st->fft_table, st->e, st->E);
    for (i=0;i<st->frame_size;i++)
       st->y[i] = 0;
    for (i=0;i<st->frame_size;i++)
       st->y[i] = 0;
-   for (i=0;i<N;i++)
-      st->Y[i] = st->y[i];
-   spx_drft_forward(st->fft_lookup, st->Y);
-   
+   spx_fft_float(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);
    power_spectrum(st->Y, st->Yf, N);
    /* 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);
@@ -302,7 +290,7 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
          st->Eh[j] = .95*Eh - E;
          st->Yh[j] = .95*Yh - Y;
       }
          st->Eh[j] = .95*Eh - E;
          st->Yh[j] = .95*Yh - Y;
       }
-      alpha = .02*Syy / (1+See);
+      alpha = .02*Syy / (1e4+See);
       if (alpha > .02)
          alpha = .02;
       st->Pey = (1-alpha)*st->Pey + alpha*Pey;
       if (alpha > .02)
          alpha = .02;
       st->Pey = (1-alpha)*st->Pey + alpha*Pey;
@@ -348,7 +336,7 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
       }
    } else {
       for (i=0;i<=st->frame_size;i++)
       }
    } else {
       for (i=0;i<=st->frame_size;i++)
-         st->power_1[i] = adapt_rate/(1.f+st->power[i]);      
+         st->power_1[i] = adapt_rate/(1.f+st->power[i]);
    }
 
    /* Compute weight gradient */
    }
 
    /* Compute weight gradient */
@@ -367,14 +355,13 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
       /* Remove the "if" to make this an MDF filter */
       if (j==M-1 || st->cancel_count%(M-1) == j)
       {
       /* Remove the "if" to make this an MDF filter */
       if (j==M-1 || st->cancel_count%(M-1) == j)
       {
-         spx_drft_backward(st->fft_lookup, &st->W[j*N]);
-         for (i=0;i<N;i++)
-            st->W[j*N+i]*=scale;
+         float w[N];
+         spx_ifft_float(st->fft_table, &st->W[j*N], w);
          for (i=st->frame_size;i<N;i++)
          {
          for (i=st->frame_size;i<N;i++)
          {
-            st->W[j*N+i]=0;
+            w[i]=0;
          }
          }
-         spx_drft_forward(st->fft_lookup, &st->W[j*N]);
+         spx_fft_float(st->fft_table, w, &st->W[j*N]);
       }
    }
 
       }
    }
 
@@ -396,10 +383,10 @@ void speex_echo_cancel(SpeexEchoState *st, short *ref, short *echo, short *out,
       
       /* Apply hanning window (should pre-compute it)*/
       for (i=0;i<N;i++)
       
       /* Apply hanning window (should pre-compute it)*/
       for (i=0;i<N;i++)
-         st->Yps[i] = (.5-.5*cos(2*M_PI*i/N))*st->last_y[i];
+         st->y[i] = (.5-.5*cos(2*M_PI*i/N))*st->last_y[i];
       
       /* Compute power spectrum of the echo */
       
       /* Compute power spectrum of the echo */
-      spx_drft_forward(st->fft_lookup, st->Yps);
+      spx_fft_float(st->fft_table, st->y, st->Yps);
       power_spectrum(st->Yps, st->Yps, N);
       
       /* Estimate residual echo */
       power_spectrum(st->Yps, st->Yps, N);
       
       /* Estimate residual echo */