Changeset 6129

Timestamp:

Jan 9, 2020 9:05:50 AM (5 years ago)

Author:

ming

Message:

Closed #589: Update Speex AEC to the latest version to get multichannel EC

Location:

pjproject/trunk

Files:

• pjmedia/src/pjmedia/echo_common.c (modified) (1 diff)
• pjmedia/src/pjmedia/echo_speex.c (modified) (12 diffs)
• third_party/build/speex/Makefile (modified) (1 diff)
• third_party/speex/include/speex/speex_buffer.h (modified) (1 diff)
• third_party/speex/include/speex/speex_echo.h (modified) (1 diff)
• third_party/speex/include/speex/speex_jitter.h (modified) (1 diff)
• third_party/speex/include/speex/speex_preprocess.h (modified) (1 diff)
• third_party/speex/include/speex/speex_resampler.h (modified) (3 diffs)
• third_party/speex/include/speex/speexdsp_types.h (added)
• third_party/speex/libspeex/arch.h (modified) (8 diffs)
• third_party/speex/libspeex/bfin.h (added)
• third_party/speex/libspeex/buffer.c (modified) (4 diffs)
• third_party/speex/libspeex/fftwrap.c (modified) (5 diffs)
• third_party/speex/libspeex/fixed_bfin.h (modified) (10 diffs)
• third_party/speex/libspeex/fixed_generic.h (modified) (2 diffs)
• third_party/speex/libspeex/jitter.c (modified) (3 diffs)
• third_party/speex/libspeex/math_approx.h (modified) (3 diffs)
• third_party/speex/libspeex/mdf.c (modified) (45 diffs)
• third_party/speex/libspeex/misc_bfin.h (modified) (4 diffs)
• third_party/speex/libspeex/os_support.h (modified) (4 diffs)
• third_party/speex/libspeex/resample.c (modified) (55 diffs)
• third_party/speex/libspeex/resample_neon.h (added)
• third_party/speex/libspeex/resample_sse.h (modified) (4 diffs)
• third_party/speex/libspeex/scal.c (modified) (2 diffs)
• third_party/speex/libspeex/testresample2.c (added)

pjproject/trunk/pjmedia/src/pjmedia/echo_common.c

@@ -100 +100 @@
 static struct ec_operations speex_aec_op = 
 {
-    "AEC",
+    "Speex AEC",
     &speex_aec_create,
     &speex_aec_destroy,

pjproject/trunk/pjmedia/src/pjmedia/echo_speex.c

@@ -36 +36 @@
 {
     SpeexEchoState       *state;
-    SpeexPreprocessState *preprocess;
+    SpeexDecorrState     *decorr;
+    SpeexPreprocessState **preprocess;
 
     unsigned              samples_per_frame;
-    unsigned              prefetch;
+    unsigned              channel_count;
+    unsigned              spf_per_channel;
     unsigned              options;
     pj_int16_t           *tmp_frame;
@@ -58 +60 @@
 {
     speex_ec *echo;
-    int sampling_rate;
+    int i, sampling_rate;
 
     *p_echo = NULL;
@@ -65 +67 @@
     PJ_ASSERT_RETURN(echo != NULL, PJ_ENOMEM);
 
+    echo->channel_count = channel_count;
     echo->samples_per_frame = samples_per_frame;
+    echo->spf_per_channel = samples_per_frame / channel_count;
     echo->options = options;
 
-#if 0
-    echo->state = speex_echo_state_init_mc(echo->samples_per_frame,
+#if 1
+    echo->state = speex_echo_state_init_mc(echo->spf_per_channel,
                                            clock_rate * tail_ms / 1000,
                                            channel_count, channel_count);
@@ -80 +84 @@
                                         clock_rate * tail_ms / 1000);
 #endif
-    if (echo->state == NULL) {
+    if (echo->state == NULL)
         return PJ_ENOMEM;
-    }
+
+    echo->decorr = speex_decorrelate_new(clock_rate, channel_count,
+                                         echo->spf_per_channel);
+    if (echo->decorr == NULL)
+        return PJ_ENOMEM;
 
     /* Set sampling rate */
@@ -89 +97 @@
                    &sampling_rate);
 
-    echo->preprocess = speex_preprocess_state_init(echo->samples_per_frame,
-                                                   clock_rate);
-    if (echo->preprocess == NULL) {
-        speex_echo_state_destroy(echo->state);
+    /* We need to create one state per channel processed. */
+    echo->preprocess = PJ_POOL_ZALLOC_T(pool, SpeexPreprocessState *);
+    for (i = 0; i < channel_count; i++) {
+        spx_int32_t enabled;
+
+        echo->preprocess[i] = speex_preprocess_state_init(
+                                  echo->spf_per_channel, clock_rate);
+        if (echo->preprocess[i] == NULL) {
+            speex_aec_destroy(echo);
+            return PJ_ENOMEM;
+        }
+
+        /* Enable/disable AGC & denoise */
+        enabled = PJMEDIA_SPEEX_AEC_USE_AGC;
+        speex_preprocess_ctl(echo->preprocess[i], SPEEX_PREPROCESS_SET_AGC,
+                             &enabled);
+
+        enabled = PJMEDIA_SPEEX_AEC_USE_DENOISE;
+        speex_preprocess_ctl(echo->preprocess[i],
+                             SPEEX_PREPROCESS_SET_DENOISE, &enabled);
+
+        /* Currently, VAD and dereverb are set at default setting. */
+        /*
+        enabled = 1;
+        speex_preprocess_ctl(echo->preprocess[i], SPEEX_PREPROCESS_SET_VAD,
+                             &enabled);
+        speex_preprocess_ctl(echo->preprocess[i],
+                             SPEEX_PREPROCESS_SET_DEREVERB,
+                             &enabled);
+        */
+
+        /* Control echo cancellation in the preprocessor */
+        speex_preprocess_ctl(echo->preprocess[i],
+                             SPEEX_PREPROCESS_SET_ECHO_STATE, echo->state);
+    }
+
+    /* Create temporary frame for echo cancellation */
+    echo->tmp_frame = (pj_int16_t*) pj_pool_zalloc(pool, sizeof(pj_int16_t) *
+                                                   channel_count *
+                                                   samples_per_frame);
+    if (!echo->tmp_frame) {
+        speex_aec_destroy(echo);
         return PJ_ENOMEM;
     }
-
-    /* Disable all preprocessing, we only want echo cancellation */
-#if 0
-    disabled = 0;
-    enabled = 1;
-    speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_DENOISE, 
-                         &enabled);
-    speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_AGC, 
-                         &disabled);
-    speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_VAD, 
-                         &disabled);
-    speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_DEREVERB, 
-                         &enabled);
-#endif
-
-    /* Enable/disable AGC & denoise */
-    {
-        spx_int32_t enabled;
-
-        enabled = PJMEDIA_SPEEX_AEC_USE_AGC;
-        speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_AGC, 
-                             &enabled);
-
-        enabled = PJMEDIA_SPEEX_AEC_USE_DENOISE;
-        speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_DENOISE, 
-                             &enabled);
-    }
-
-    /* Control echo cancellation in the preprocessor */
-   speex_preprocess_ctl(echo->preprocess, SPEEX_PREPROCESS_SET_ECHO_STATE, 
-                        echo->state);
-
-
-    /* Create temporary frame for echo cancellation */
-    echo->tmp_frame = (pj_int16_t*) pj_pool_zalloc(pool, 2*samples_per_frame);
-    PJ_ASSERT_RETURN(echo->tmp_frame != NULL, PJ_ENOMEM);
 
     /* Done */
@@ -145 +155 @@
 {
     speex_ec *echo = (speex_ec*) state;
+    unsigned i;
 
     PJ_ASSERT_RETURN(echo && echo->state, PJ_EINVAL);
@@ -153 +164 @@
     }
 
+    if (echo->decorr) {
+        speex_decorrelate_destroy(echo->decorr);
+        echo->decorr = NULL;
+    }
+
     if (echo->preprocess) {
-        speex_preprocess_state_destroy(echo->preprocess);
+        for (i = 0; i < echo->channel_count; i++) {
+            if (echo->preprocess[i]) {
+                speex_preprocess_state_destroy(echo->preprocess[i]);
+                echo->preprocess[i] = NULL;
+            }
+        }
         echo->preprocess = NULL;
     }
@@ -182 +203 @@
 {
     speex_ec *echo = (speex_ec*) state;
+    unsigned i;
 
     /* Sanity checks */
@@ -193 +215 @@
 
 
-    /* Preprocess output */
-    speex_preprocess_run(echo->preprocess, (spx_int16_t*)echo->tmp_frame);
+    /* Preprocess output per channel */
+    for (i = 0; i < echo->channel_count; i++) {
+        spx_int16_t *buf = (spx_int16_t*)echo->tmp_frame;
+        unsigned j;
+
+        /* De-interleave each channel. */
+        if (echo->channel_count > 1) {
+            for (j = 0; j < echo->spf_per_channel; j++) {
+                rec_frm[j] = echo->tmp_frame[j * echo->channel_count + i];
+            }
+            buf = (spx_int16_t*)rec_frm;
+        }
+
+        speex_preprocess_run(echo->preprocess[i], buf);
+
+        if (echo->channel_count > 1) {
+            for (j = 0; j < echo->spf_per_channel; j++) {
+                echo->tmp_frame[j * echo->channel_count + i] = rec_frm[j];
+            }
+        }
+    }
 
     /* Copy temporary buffer back to original rec_frm */
@@ -213 +254 @@
     /* Sanity checks */
     PJ_ASSERT_RETURN(echo && play_frm, PJ_EINVAL);
+
+    /* Channel decorrelation algorithm is useful for multi-channel echo
+     * cancellation only .
+     */
+    if (echo->channel_count > 1) {
+        pjmedia_copy_samples(echo->tmp_frame, play_frm, echo->samples_per_frame);
+        speex_decorrelate(echo->decorr, (spx_int16_t*)echo->tmp_frame,
+                          (spx_int16_t*)play_frm, 100);
+    }
 
     speex_echo_playback(echo->state, (spx_int16_t*)play_frm);
@@ -228 +278 @@
 {
     speex_ec *echo = (speex_ec*) state;
+    unsigned i;
 
     /* Sanity checks */
@@ -235 +286 @@
 
     /* Cancel echo */
-    pjmedia_copy_samples(echo->tmp_frame, rec_frm, echo->samples_per_frame);
     speex_echo_capture(echo->state,
-                       (spx_int16_t*)echo->tmp_frame,
-                       (spx_int16_t*)rec_frm);
-
-    /* Apply preprocessing */
-    speex_preprocess_run(echo->preprocess, (spx_int16_t*)rec_frm);
+                       (spx_int16_t*)rec_frm,
+                       (spx_int16_t*)echo->tmp_frame);
+
+    /* Apply preprocessing per channel. */
+    for (i = 0; i < echo->channel_count; i++) {
+        spx_int16_t *buf = (spx_int16_t*)echo->tmp_frame;
+        unsigned j;
+
+        /* De-interleave each channel. */
+        if (echo->channel_count > 1) {
+            for (j = 0; j < echo->spf_per_channel; j++) {
+                rec_frm[j] = echo->tmp_frame[j * echo->channel_count + i];
+            }
+            buf = (spx_int16_t*)rec_frm;
+        }
+
+        speex_preprocess_run(echo->preprocess[i], buf);
+
+        if (echo->channel_count > 1) {
+            for (j = 0; j < echo->spf_per_channel; j++) {
+                echo->tmp_frame[j * echo->channel_count + i] = rec_frm[j];
+            }
+        }
+    }
+
+    pjmedia_copy_samples(rec_frm, echo->tmp_frame, echo->samples_per_frame);
 
     return PJ_SUCCESS;

pjproject/trunk/third_party/build/speex/Makefile

@@ -40 +40 @@
                         quant_lsp.o resample.o sb_celp.o smallft.o \
                         speex.o speex_callbacks.o speex_header.o \
-                        stereo.o vbr.o vq.o window.o
+                        scal.o stereo.o vbr.o vq.o window.o
 
 export SPEEX_CFLAGS = -DHAVE_CONFIG_H $(_CFLAGS)

pjproject/trunk/third_party/speex/include/speex/speex_buffer.h

@@ -35 +35 @@
 #define SPEEX_BUFFER_H
 
-#include "speex/speex_types.h"
+#include "speexdsp_types.h"
 
 #ifdef __cplusplus

pjproject/trunk/third_party/speex/include/speex/speex_echo.h

@@ -38 +38 @@
  *  @{
  */
-#include "speex/speex_types.h"
+#include "speexdsp_types.h"
 
 #ifdef __cplusplus

pjproject/trunk/third_party/speex/include/speex/speex_jitter.h

@@ -42 +42 @@
  */
 
-#include "speex/speex_types.h"
+#include "speexdsp_types.h"
 
 #ifdef __cplusplus

pjproject/trunk/third_party/speex/include/speex/speex_preprocess.h

@@ -44 +44 @@
  */
 
-#include "speex/speex_types.h"
+#include "speexdsp_types.h"
 
 #ifdef __cplusplus

pjproject/trunk/third_party/speex/include/speex/speex_resampler.h

@@ -83 +83 @@
 #define spx_uint32_t unsigned int
       
+#define speex_assert(cond)
+
 #else /* OUTSIDE_SPEEX */
 
-#include "speex/speex_types.h"
+#include "speexdsp_types.h"
 
 #endif /* OUTSIDE_SPEEX */
@@ -105 +107 @@
    RESAMPLER_ERR_INVALID_ARG     = 3,
    RESAMPLER_ERR_PTR_OVERLAP     = 4,
+   RESAMPLER_ERR_OVERFLOW        = 5,
    
    RESAMPLER_ERR_MAX_ERROR
@@ -303 +306 @@
                                       spx_uint32_t *stride);
 
-/** Get the latency in input samples introduced by the resampler.
+/** Get the latency introduced by the resampler measured in input samples.
  * @param st Resampler state
  */
 int speex_resampler_get_input_latency(SpeexResamplerState *st);
 
-/** Get the latency in output samples introduced by the resampler.
+/** Get the latency introduced by the resampler measured in output samples.
  * @param st Resampler state
  */

pjproject/trunk/third_party/speex/libspeex/arch.h

@@ -8 +8 @@
    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
@@ -50 +50 @@
 #error You cannot compile as floating point and fixed point at the same time
 #endif
-#ifdef _USE_SSE
+#ifdef USE_SSE
 #error SSE is only for floating-point
 #endif
@@ -76 +76 @@
 
 #ifndef OUTSIDE_SPEEX
-#include "speex/speex_types.h"
+#include "speex/speexdsp_types.h"
 #endif
 
@@ -90 +90 @@
 
 typedef spx_int16_t spx_word16_t;
-typedef spx_int32_t   spx_word32_t;
+typedef spx_int32_t spx_word32_t;
 typedef spx_word32_t spx_mem_t;
 typedef spx_word16_t spx_coef_t;
@@ -109 +109 @@
 #define SIG_SHIFT    14
 #define GAIN_SHIFT   6
+
+#define WORD2INT(x) ((x) < -32767 ? -32768 : ((x) > 32766 ? 32767 : (x)))
 
 #define VERY_SMALL 0
@@ -172 +174 @@
 #define SATURATE16(x,a) (x)
 #define SATURATE32(x,a) (x)
+#define SATURATE32PSHR(x,shift,a) (x)
 
 #define PSHR(a,shift)       (a)
@@ -211 +214 @@
 #define PDIV32(a,b)     (((spx_word32_t)(a))/(spx_word32_t)(b))
 
-
+#define WORD2INT(x) ((x) < -32767.5f ? -32768 : \
+                    ((x) > 32766.5f ? 32767 : (spx_int16_t)floor(.5 + (x))))
 #endif
 
@@ -218 +222 @@
 
 /* 2 on TI C5x DSP */
-#define BYTES_PER_CHAR 2 
+#define BYTES_PER_CHAR 2
 #define BITS_PER_CHAR 16
 #define LOG2_BITS_PER_CHAR 4
 
-#else 
+#else
 
 #define BYTES_PER_CHAR 1

pjproject/trunk/third_party/speex/libspeex/buffer.c

@@ -1 +1 @@
 /* Copyright (C) 2007 Jean-Marc Valin
-      
+
    File: buffer.c
    This is a very simple ring buffer implementation. It is not thread-safe
@@ -39 +39 @@
 #include "os_support.h"
 #include "arch.h"
-#include <speex/speex_buffer.h>
+#include "speex/speex_buffer.h"
 
 struct SpeexBuffer_ {
@@ -100 +100 @@
 EXPORT int speex_buffer_writezeros(SpeexBuffer *st, int len)
 {
-   /* This is almost the same as for speex_buffer_write() but using 
+   /* This is almost the same as for speex_buffer_write() but using
    SPEEX_MEMSET() instead of SPEEX_COPY(). Update accordingly. */
    int end;
@@ -136 +136 @@
    if (len > st->available)
    {
-      SPEEX_MEMSET(data+st->available, 0, st->size-st->available);
+      SPEEX_MEMSET(data+st->available, 0, len - st->available);
       len = st->available;
    }

pjproject/trunk/third_party/speex/libspeex/fftwrap.c

@@ -1 @@
-/* Copyright (C) 2005-2006 Jean-Marc Valin 
+/* Copyright (C) 2005-2006 Jean-Marc Valin
    File: fftwrap.c
 
-   Wrapper for various FFTs 
+   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
@@ -63 +63 @@
    {
       out[i] = SHL16(in[i], shift);
-   }   
+   }
    return shift;
 }
@@ -166 +166 @@
 }
 
+#elif defined(USE_INTEL_IPP)
+
+#include <ipps.h>
+
+struct ipp_fft_config
+{
+  IppsDFTSpec_R_32f *dftSpec;
+  Ipp8u *buffer;
+};
+
+void *spx_fft_init(int size)
+{
+  int bufferSize = 0;
+  int hint;
+  struct ipp_fft_config *table;
+
+  table = (struct ipp_fft_config *)speex_alloc(sizeof(struct ipp_fft_config));
+
+  /* there appears to be no performance difference between ippAlgHintFast and
+     ippAlgHintAccurate when using the with the floating point version
+     of the fft. */
+  hint = ippAlgHintAccurate;
+
+  ippsDFTInitAlloc_R_32f(&table->dftSpec, size, IPP_FFT_DIV_FWD_BY_N, hint);
+
+  ippsDFTGetBufSize_R_32f(table->dftSpec, &bufferSize);
+  table->buffer = ippsMalloc_8u(bufferSize);
+
+  return table;
+}
+
+void spx_fft_destroy(void *table)
+{
+  struct ipp_fft_config *t = (struct ipp_fft_config *)table;
+  ippsFree(t->buffer);
+  ippsDFTFree_R_32f(t->dftSpec);
+  speex_free(t);
+}
+
+void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+  struct ipp_fft_config *t = (struct ipp_fft_config *)table;
+  ippsDFTFwd_RToPack_32f(in, out, t->dftSpec, t->buffer);
+}
+
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
+{
+  struct ipp_fft_config *t = (struct ipp_fft_config *)table;
+  ippsDFTInv_PackToR_32f(in, out, t->dftSpec, t->buffer);
+}
+
 #elif defined(USE_GPL_FFTW3)
 
@@ -220 +271 @@
 }
 
-void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 
+void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
 {
   int i;
@@ -235 +286 @@
 
   fftwf_execute(t->ifft);
-  
+
   for(i=0;i<N;++i)
     out[i] = optr[i];

pjproject/trunk/third_party/speex/libspeex/fixed_bfin.h

@@ -9 +9 @@
    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
@@ -36 +36 @@
 #ifndef FIXED_BFIN_H
 #define FIXED_BFIN_H
+
+#include "bfin.h"
 
 #undef PDIV32_16
@@ -58 +60 @@
    : "=m" (res)
    : "m" (a), "m" (bb)
-   : "P0", "R0", "R1", "cc");
+   : "P0", "R0", "R1", "ASTAT" BFIN_HWLOOP0_REGS);
    return res;
 }
@@ -67 +69 @@
    spx_word32_t res, bb;
    bb = b;
-   /* Make the roundinf consistent with the C version 
+   /* Make the roundinf consistent with the C version
       (do we need to do that?)*/
-   if (a<0) 
+   if (a<0)
       a += (b-1);
    __asm__  (
@@ -85 +87 @@
    : "=m" (res)
    : "m" (a), "m" (bb)
-   : "P0", "R0", "R1", "cc");
+   : "P0", "R0", "R1", "ASTAT" BFIN_HWLOOP0_REGS);
    return res;
 }
@@ -99 +101 @@
    : "=d" (res)
    : "%d" (a), "d" (b)
+   : "ASTAT"
    );
    return res;
@@ -114 +117 @@
    : "=&W" (res), "=&d" (b)
    : "d" (a), "1" (b)
-   : "A1"
+   : "A1", "ASTAT"
    );
    return res;
@@ -131 +134 @@
    : "=&W" (res), "=&d" (b)
    : "d" (a), "1" (b), "d" (c)
-   : "A1"
+   : "A1", "ASTAT"
          );
    return res;
@@ -148 +151 @@
    : "=W" (res), "=d" (a), "=d" (b)
    : "1" (a), "2" (b)
-   : "A1"
+   : "A1", "ASTAT"
          );
    return res;
@@ -166 +169 @@
    : "=&W" (res), "=&d" (b)
    : "d" (a), "1" (b), "d" (c)
-   : "A1"
+   : "A1", "ASTAT"
          );
    return res;

pjproject/trunk/third_party/speex/libspeex/fixed_generic.h

@@ -8 +8 @@
    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
@@ -52 +52 @@
 #define SATURATE16(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x)))
 #define SATURATE32(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x)))
+
+#define SATURATE32PSHR(x,shift,a) (((x)>=(SHL32(a,shift))) ? (a) : \
+                                   (x)<=-(SHL32(a,shift)) ? -(a) : \
+                                   (PSHR32(x, shift)))
 
 #define SHR(a,shift) ((a) >> (shift))

pjproject/trunk/third_party/speex/libspeex/jitter.c

@@ -56 +56 @@
 
 #include "arch.h"
-#include <speex/speex.h>
-#include <speex/speex_bits.h>
 #include <speex/speex_jitter.h>
 #include "os_support.h"
@@ -467 +465 @@
    int i;
    unsigned int j;
-   int incomplete = 0;
    spx_int16_t opt;
    
@@ -572 +569 @@
       {
          i=besti;
-         incomplete = 1;
          /*fprintf (stderr, "incomplete: %d %d %d %d\n", jitter->packets[i].timestamp, jitter->pointer_timestamp, chunk_size, jitter->packets[i].span);*/
       }

pjproject/trunk/third_party/speex/libspeex/math_approx.h

@@ -8 +8 @@
    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
@@ -169 +169 @@
    }
    x = SUB16(16384,x);
-   
+
    x = x >> 1;
    sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3))))));
    ret = spx_sqrt(SHL32(EXTEND32(sq),13));
-   
+
    /*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/
    if (s)
@@ -209 +209 @@
 {
    spx_word16_t x2;
-   
+
    x2 = MULT16_16_P15(x,x);
    return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2))))))));

pjproject/trunk/third_party/speex/libspeex/mdf.c

@@ -1 @@
-/* Copyright (C) 2003-2006 Jean-Marc Valin
+/* Copyright (C) 2003-2008 Jean-Marc Valin
 
    File: mdf.c
@@ -34 +34 @@
    The echo canceller is based on the MDF algorithm described in:
 
-   J. S. Soo, K. K. Pang Multidelay block frequency adaptive filter, 
-   IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, No. 2, 
+   J. S. Soo, K. K. Pang Multidelay block frequency adaptive filter,
+   IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, No. 2,
    February 1990.
-   
-   We use the Alternatively Updated MDF (AUMDF) variant. Robustness to 
+
+   We use the Alternatively Updated MDF (AUMDF) variant. Robustness to
    double-talk is achieved using a variable learning rate as described in:
-   
-   Valin, J.-M., On Adjusting the Learning Rate in Frequency Domain Echo 
+
+   Valin, J.-M., On Adjusting the Learning Rate in Frequency Domain Echo
    Cancellation With Double-Talk. IEEE Transactions on Audio,
    Speech and Language Processing, Vol. 15, No. 3, pp. 1030-1034, 2007.
    http://people.xiph.org/~jm/papers/valin_taslp2006.pdf
-   
+
    There is no explicit double-talk detection, but a continuous variation
    in the learning rate based on residual echo, double-talk and background
    noise.
-   
+
    About the fixed-point version:
-   All the signals are represented with 16-bit words. The filter weights 
+   All the signals are represented with 16-bit words. The filter weights
    are represented with 32-bit words, but only the top 16 bits are used
    in most cases. The lower 16 bits are completely unreliable (due to the
@@ -57 +57 @@
    adaptation -- probably by removing a "threshold effect" due to
    quantization (rounding going to zero) when the gradient is small.
-   
+
    Another kludge that seems to work good: when performing the weight
    update, we only move half the way toward the "goal" this seems to
@@ -63 +63 @@
    can be seen as applying a gradient descent on a "soft constraint"
    instead of having a hard constraint.
-   
+
 */
 
@@ -132 +132 @@
    int saturated;
    int screwed_up;
+   int C;                    /** Number of input channels (microphones) */
+   int K;                    /** Number of output channels (loudspeakers) */
    spx_int32_t sampling_rate;
    spx_word16_t spec_average;
@@ -138 +140 @@
    spx_word32_t sum_adapt;
    spx_word16_t leak_estimate;
-   
+
    spx_word16_t *e;      /* scratch */
    spx_word16_t *x;      /* Far-end input buffer (2N) */
@@ -172 +174 @@
    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];
+   spx_mem_t *notch_mem;
 
    /* NOTE: If you only use speex_echo_cancel() and want to save some memory, remove this */
@@ -183 +185 @@
 };
 
-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;
@@ -191 +193 @@
 #else
    den2 = radius*radius + .7*(1-radius)*(1-radius);
-#endif   
+#endif
    /*printf ("%d %d %d %d %d %d\n", num[0], num[1], num[2], den[0], den[1], den[2]);*/
    for (i=0;i<len;i++)
    {
-      spx_word16_t vin = in[i];
+      spx_word16_t vin = in[i*stride];
       spx_word32_t vout = mem[0] + SHL32(EXTEND32(vin),15);
 #ifdef FIXED_POINT
@@ -235 +237 @@
 }
 
+/** 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;i<N-1;i+=2,j++)
+   {
+      ps[j] +=  MULT16_16(X[i],X[i]) + MULT16_16(X[i+1],X[i+1]);
+   }
+   ps[j]+=MULT16_16(X[i],X[i]);
+}
+
 /** Compute cross-power spectrum of a half-complex (packed) vectors and add to acc */
 #ifdef FIXED_POINT
@@ -331 +345 @@
 }
 
-static inline void mdf_adjust_prop(const spx_word32_t *W, int N, int M, spx_word16_t *prop)
-{
-   int i, j;
+static inline void mdf_adjust_prop(const spx_word32_t *W, int N, int M, int P, spx_word16_t *prop)
+{
+   int i, j, p;
    spx_word16_t max_sum = 1;
    spx_word32_t prop_sum = 1;
@@ -339 +353 @@
    {
       spx_word32_t tmp = 1;
-      for (j=0;j<N;j++)
-         tmp += MULT16_16(EXTRACT16(SHR32(W[i*N+j],18)), EXTRACT16(SHR32(W[i*N+j],18)));
+      for (p=0;p<P;p++)
+         for (j=0;j<N;j++)
+            tmp += MULT16_16(EXTRACT16(SHR32(W[p*N*M + i*N+j],18)), EXTRACT16(SHR32(W[p*N*M + i*N+j],18)));
 #ifdef FIXED_POINT
       /* Just a security in case an overflow were to occur */
@@ -379 +394 @@
 
 /** Creates a new echo canceller state */
-SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length)
-{
-   int i,N,M;
+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));
 
+   st->K = nb_speakers;
+   st->C = nb_mic;
+   C=st->C;
+   K=st->K;
 #ifdef DUMP_ECHO_CANCEL_DATA
    if (rFile || pFile || oFile)
@@ -391 +415 @@
    oFile = fopen("aec_out.sw", "wb");
 #endif
-   
+
    st->frame_size = frame_size;
    st->window_size = 2*frame_size;
@@ -413 +437 @@
 
    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->input = (spx_word16_t*)speex_alloc(st->frame_size*sizeof(spx_word16_t));
-   st->y = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t));
-   st->last_y = (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->input = (spx_word16_t*)speex_alloc(C*st->frame_size*sizeof(spx_word16_t));
+   st->y = (spx_word16_t*)speex_alloc(C*N*sizeof(spx_word16_t));
+   st->last_y = (spx_word16_t*)speex_alloc(C*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));
@@ -425 +449 @@
    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->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));
 #ifdef TWO_PATH
-   st->foreground = (spx_word16_t*)speex_alloc(M*N*sizeof(spx_word16_t));
+   st->foreground = (spx_word16_t*)speex_alloc(M*N*C*K*sizeof(spx_word16_t));
 #endif
    st->PHI = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
@@ -451 +475 @@
    for (i=0;i<=st->frame_size;i++)
       st->power_1[i] = FLOAT_ONE;
-   for (i=0;i<N*M;i++)
+   for (i=0;i<N*M*K*C;i++)
       st->W[i] = 0;
    {
@@ -466 +490 @@
       for (i=M-1;i>=0;i--)
       {
-         st->prop[i] = DIV32(MULT16_16(QCONST16(.8,15), st->prop[i]),sum);
-      }
-   }
-   
-   st->memX=st->memD=st->memE=0;
+         st->prop[i] = DIV32(MULT16_16(QCONST16(.8f,15), st->prop[i]),sum);
+      }
+   }
+
+   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)
@@ -479 +505 @@
       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;
-   
+
 #ifdef TWO_PATH
    st->Davg1 = st->Davg2 = 0;
    st->Dvar1 = st->Dvar2 = FLOAT_ZERO;
 #endif
-   
-   st->play_buf = (spx_int16_t*)speex_alloc((PLAYBACK_DELAY+1)*st->frame_size*sizeof(spx_int16_t));
+
+   st->play_buf = (spx_int16_t*)speex_alloc(K*(PLAYBACK_DELAY+1)*st->frame_size*sizeof(spx_int16_t));
    st->play_buf_pos = PLAYBACK_DELAY*st->frame_size;
    st->play_buf_started = 0;
-   
+
    return st;
 }
 
 /** Resets echo canceller state */
-void speex_echo_state_reset(SpeexEchoState *st)
-{
-   int i, M, N;
+EXPORT void speex_echo_state_reset(SpeexEchoState *st)
+{
+   int i, M, N, C, K;
    st->cancel_count=0;
    st->screwed_up = 0;
    N = st->window_size;
    M = st->M;
+   C=st->C;
+   K=st->K;
    for (i=0;i<N*M;i++)
       st->W[i] = 0;
@@ -522 +550 @@
       st->last_y[i] = 0;
    }
-   for (i=0;i<N;i++)
+   for (i=0;i<N*C;i++)
    {
       st->E[i] = 0;
+   }
+   for (i=0;i<N*K;i++)
+   {
       st->x[i] = 0;
    }
-   st->notch_mem[0] = st->notch_mem[1] = 0;
-   st->memX=st->memD=st->memE=0;
+   for (i=0;i<2*C;i++)
+      st->notch_mem[i] = 0;
+   for (i=0;i<C;i++)
+      st->memD[i]=st->memE[i]=0;
+   for (i=0;i<K;i++)
+      st->memX[i]=0;
 
    st->saturated = 0;
@@ -546 +581 @@
 
 /** Destroys an echo canceller state */
-void speex_echo_state_destroy(SpeexEchoState *st)
+EXPORT void speex_echo_state_destroy(SpeexEchoState *st)
 {
    spx_fft_destroy(st->fft_table);
@@ -577 +612 @@
    speex_free(st->wtmp2);
 #endif
+   speex_free(st->memX);
+   speex_free(st->memD);
+   speex_free(st->memE);
+   speex_free(st->notch_mem);
+
    speex_free(st->play_buf);
    speex_free(st);
-   
+
 #ifdef DUMP_ECHO_CANCEL_DATA
    fclose(rFile);
@@ -588 +628 @@
 }
 
-void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out)
+EXPORT void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out)
 {
    int i;
@@ -611 +651 @@
 }
 
-void speex_echo_playback(SpeexEchoState *st, const spx_int16_t *play)
+EXPORT void speex_echo_playback(SpeexEchoState *st, const spx_int16_t *play)
 {
    /*speex_warning_int("playback with fill level ", st->play_buf_pos/st->frame_size);*/
@@ -638 +678 @@
 
 /** Performs echo cancellation on a frame (deprecated, last arg now ignored) */
-void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *in, const spx_int16_t *far_end, spx_int16_t *out, spx_int32_t *Yout)
+EXPORT void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *in, const spx_int16_t *far_end, spx_int16_t *out, spx_int32_t *Yout)
 {
    speex_echo_cancellation(st, in, far_end, out);
@@ -644 +684 @@
 
 /** Performs echo cancellation on a frame */
-void speex_echo_cancellation(SpeexEchoState *st, const spx_int16_t *in, const spx_int16_t *far_end, spx_int16_t *out)
-{
-   int i,j;
-   int N,M;
+EXPORT void speex_echo_cancellation(SpeexEchoState *st, const spx_int16_t *in, const spx_int16_t *far_end, spx_int16_t *out)
+{
+   int i,j, chan, speak;
+   int N,M, C, K;
    spx_word32_t Syy,See,Sxx,Sdd, Sff;
 #ifdef TWO_PATH
@@ -659 +699 @@
    spx_word16_t RER;
    spx_word32_t tmp32;
-   
+
    N = st->window_size;
    M = st->M;
+   C = st->C;
+   K = st->K;
+
    st->cancel_count++;
 #ifdef FIXED_POINT
@@ -671 +714 @@
 #endif
 
-   /* Apply a notch filter to make sure DC doesn't end up causing problems */
-   filter_dc_notch16(in, st->notch_radius, st->input, st->frame_size, st->notch_mem);
-   /* Copy input data to buffer and apply pre-emphasis */
-   for (i=0;i<st->frame_size;i++)
-   {
-      spx_word32_t tmp32;
-      tmp32 = SUB32(EXTEND32(far_end[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX)));
-#ifdef FIXED_POINT
-      /* If saturation occurs here, we need to freeze adaptation for M+1 frames (not just one) */
-      if (tmp32 > 32767)
-      {
-         tmp32 = 32767;
-         st->saturated = M+1;
-      }
-      if (tmp32 < -32767)
-      {
-         tmp32 = -32767;
-         st->saturated = M+1;
-      }      
-#endif
-      st->x[i+st->frame_size] = EXTRACT16(tmp32);
-      st->memX = far_end[i];
-      
-      tmp32 = SUB32(EXTEND32(st->input[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memD)));
-#ifdef FIXED_POINT
-      if (tmp32 > 32767)
-      {
-         tmp32 = 32767;
-         if (st->saturated == 0)
-            st->saturated = 1;
-      }      
-      if (tmp32 < -32767)
-      {
-         tmp32 = -32767;
-         if (st->saturated == 0)
-            st->saturated = 1;
-      }
-#endif
-      st->memD = st->input[i];
-      st->input[i] = tmp32;
-   }
-
-   /* Shift memory: this could be optimized eventually*/
-   for (j=M-1;j>=0;j--)
-   {
-      for (i=0;i<N;i++)
-         st->X[(j+1)*N+i] = st->X[j*N+i];
-   }
-
-   /* Convert x (far end) to frequency domain */
-   spx_fft(st->fft_table, st->x, &st->X[0]);
-   for (i=0;i<N;i++)
-      st->last_y[i] = st->x[i];
-   Sxx = mdf_inner_prod(st->x+st->frame_size, st->x+st->frame_size, st->frame_size);
-   for (i=0;i<st->frame_size;i++)
-      st->x[i] = st->x[i+st->frame_size];
-   /* From here on, the top part of x is used as scratch space */
-   
+   for (chan = 0; chan < C; chan++)
+   {
+      /* Apply a notch filter to make sure DC doesn't end up causing problems */
+      filter_dc_notch16(in+chan, st->notch_radius, st->input+chan*st->frame_size, st->frame_size, st->notch_mem+2*chan, C);
+      /* Copy input data to buffer and apply pre-emphasis */
+      /* Copy input data to buffer */
+      for (i=0;i<st->frame_size;i++)
+      {
+         spx_word32_t tmp32;
+         /* FIXME: This core has changed a bit, need to merge properly */
+         tmp32 = SUB32(EXTEND32(st->input[chan*st->frame_size+i]), EXTEND32(MULT16_16_P15(st->preemph, st->memD[chan])));
+#ifdef FIXED_POINT
+         if (tmp32 > 32767)
+         {
+            tmp32 = 32767;
+            if (st->saturated == 0)
+               st->saturated = 1;
+         }
+         if (tmp32 < -32767)
+         {
+            tmp32 = -32767;
+            if (st->saturated == 0)
+               st->saturated = 1;
+         }
+#endif
+         st->memD[chan] = st->input[chan*st->frame_size+i];
+         st->input[chan*st->frame_size+i] = EXTRACT16(tmp32);
+      }
+   }
+
+   for (speak = 0; speak < K; speak++)
+   {
+      for (i=0;i<st->frame_size;i++)
+      {
+         spx_word32_t tmp32;
+         st->x[speak*N+i] = st->x[speak*N+i+st->frame_size];
+         tmp32 = SUB32(EXTEND32(far_end[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)
+         {
+            tmp32 = 32767;
+            st->saturated = M+1;
+         }
+         if (tmp32 < -32767)
+         {
+            tmp32 = -32767;
+            st->saturated = M+1;
+         }
+#endif
+         st->x[speak*N+i+st->frame_size] = EXTRACT16(tmp32);
+         st->memX[speak] = far_end[i*K+speak];
+      }
+   }
+
+   for (speak = 0; speak < K; speak++)
+   {
+      /* Shift memory: this could be optimized eventually*/
+      for (j=M-1;j>=0;j--)
+      {
+         for (i=0;i<N;i++)
+            st->X[(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]);
+   }
+
+   Sxx = 0;
+   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);
+   }
+
+   Sff = 0;
+   for (chan = 0; chan < C; chan++)
+   {
 #ifdef TWO_PATH
-   /* Compute foreground filter */
-   spectral_mul_accum16(st->X, st->foreground, st->Y, N, M);   
-   spx_ifft(st->fft_table, st->Y, st->e);
-   for (i=0;i<st->frame_size;i++)
-      st->e[i] = SUB16(st->input[i], st->e[i+st->frame_size]);
-   Sff = mdf_inner_prod(st->e, st->e, st->frame_size);
-#endif
-   
+      /* Compute foreground filter */
+      spectral_mul_accum16(st->X, st->foreground+chan*N*K*M, st->Y+chan*N, N, M*K);
+      spx_ifft(st->fft_table, st->Y+chan*N, st->e+chan*N);
+      for (i=0;i<st->frame_size;i++)
+         st->e[chan*N+i] = SUB16(st->input[chan*st->frame_size+i], st->e[chan*N+i+st->frame_size]);
+      Sff += mdf_inner_prod(st->e+chan*N, st->e+chan*N, st->frame_size);
+#endif
+   }
+
    /* Adjust proportional adaption rate */
-   mdf_adjust_prop (st->W, N, M, st->prop);
+   /* FIXME: Adjust that for C, K*/
+   if (st->adapted)
+      mdf_adjust_prop (st->W, N, M, C*K, st->prop);
    /* Compute weight gradient */
    if (st->saturated == 0)
    {
-      for (j=M-1;j>=0;j--)
-      {
-         weighted_spectral_mul_conj(st->power_1, FLOAT_SHL(PSEUDOFLOAT(st->prop[j]),-15), &st->X[(j+1)*N], st->E, st->PHI, N);
-         for (i=0;i<N;i++)
-            st->W[j*N+i] = ADD32(st->W[j*N+i], st->PHI[i]);
-         
+      for (chan = 0; chan < C; chan++)
+      {
+         for (speak = 0; speak < K; speak++)
+         {
+            for (j=M-1;j>=0;j--)
+            {
+               weighted_spectral_mul_conj(st->power_1, FLOAT_SHL(PSEUDOFLOAT(st->prop[j]),-15), &st->X[(j+1)*N*K+speak*N], st->E+chan*N, st->PHI, N);
+               for (i=0;i<N;i++)
+                  st->W[chan*N*K*M + j*N*K + speak*N + i] += st->PHI[i];
+            }
+         }
       }
    } else {
       st->saturated--;
    }
-   
+
+   /* FIXME: MC conversion required */
    /* Update weight to prevent circular convolution (MDF / AUMDF) */
-   for (j=0;j<M;j++)
-   {
-      /* 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;i<N;i++)
-            st->wtmp2[i] = EXTRACT16(PSHR32(st->W[j*N+i],NORMALIZE_SCALEDOWN+16));
-         spx_ifft(st->fft_table, st->wtmp2, st->wtmp);
-         for (i=0;i<st->frame_size;i++)
+   for (chan = 0; chan < C; chan++)
+   {
+      for (speak = 0; speak < K; speak++)
+      {
+         for (j=0;j<M;j++)
          {
-            st->wtmp[i]=0;
+            /* 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;i<N;i++)
+                  st->wtmp2[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;i<st->frame_size;i++)
+               {
+                  st->wtmp[i]=0;
+               }
+               for (i=st->frame_size;i<N;i++)
+               {
+                  st->wtmp[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;i<N;i++)
+                  st->W[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[chan*N*K*M + j*N*K + speak*N], st->wtmp);
+               for (i=st->frame_size;i<N;i++)
+               {
+                  st->wtmp[i]=0;
+               }
+               spx_fft(st->fft_table, st->wtmp, &st->W[chan*N*K*M + j*N*K + speak*N]);
+#endif
+            }
          }
-         for (i=st->frame_size;i<N;i++)
-         {
-            st->wtmp[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;i<N;i++)
-            st->W[j*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;i<N;i++)
-         {
-            st->wtmp[i]=0;
-         }
-         spx_fft(st->fft_table, st->wtmp, &st->W[j*N]);
-#endif
-      }
-   }
-
-   /* Compute filter response Y */
-   spectral_mul_accum(st->X, st->W, st->Y, N, M);
-   spx_ifft(st->fft_table, st->Y, st->y);
-
+      }
+   }
+
+   /* 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;
+
+   Dbf = 0;
+   See = 0;
 #ifdef TWO_PATH
    /* Difference in response, this is used to estimate the variance of our residual power estimate */
-   for (i=0;i<st->frame_size;i++)
-      st->e[i] = SUB16(st->e[i+st->frame_size], st->y[i+st->frame_size]);
-   Dbf = 10+mdf_inner_prod(st->e, st->e, st->frame_size);
-#endif
-
-   for (i=0;i<st->frame_size;i++)
-      st->e[i] = SUB16(st->input[i], st->y[i+st->frame_size]);
-   See = mdf_inner_prod(st->e, st->e, st->frame_size);
+   for (chan = 0; chan < C; chan++)
+   {
+      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);
+      for (i=0;i<st->frame_size;i++)
+         st->e[chan*N+i] = SUB16(st->e[chan*N+i+st->frame_size], st->y[chan*N+i+st->frame_size]);
+      Dbf += 10+mdf_inner_prod(st->e+chan*N, st->e+chan*N, st->frame_size);
+      for (i=0;i<st->frame_size;i++)
+         st->e[chan*N+i] = SUB16(st->input[chan*st->frame_size+i], st->y[chan*N+i+st->frame_size]);
+      See += mdf_inner_prod(st->e+chan*N, st->e+chan*N, st->frame_size);
+   }
+#endif
+
 #ifndef TWO_PATH
    Sff = See;
@@ -808 +892 @@
 #ifdef TWO_PATH
    /* Logic for updating the foreground filter */
-   
+
    /* For two time windows, compute the mean of the energy difference, as well as the variance */
    st->Davg1 = ADD32(MULT16_32_Q15(QCONST16(.6f,15),st->Davg1), MULT16_32_Q15(QCONST16(.4f,15),SUB32(Sff,See)));
@@ -814 +898 @@
    st->Dvar1 = FLOAT_ADD(FLOAT_MULT(VAR1_SMOOTH, st->Dvar1), FLOAT_MUL32U(MULT16_32_Q15(QCONST16(.4f,15),Sff), MULT16_32_Q15(QCONST16(.4f,15),Dbf)));
    st->Dvar2 = FLOAT_ADD(FLOAT_MULT(VAR2_SMOOTH, st->Dvar2), FLOAT_MUL32U(MULT16_32_Q15(QCONST16(.15f,15),Sff), MULT16_32_Q15(QCONST16(.15f,15),Dbf)));
-   
+
    /* Equivalent float code:
    st->Davg1 = .6*st->Davg1 + .4*(Sff-See);
@@ -821 +905 @@
    st->Dvar2 = .7225*st->Dvar2 + .0225*Sff*Dbf;
    */
-   
+
    update_foreground = 0;
    /* Check if we have a statistically significant reduction in the residual echo */
@@ -831 +915 @@
    else if (FLOAT_GT(FLOAT_MUL32U(st->Davg2, ABS32(st->Davg2)), FLOAT_MULT(VAR2_UPDATE,(st->Dvar2))))
       update_foreground = 1;
-   
+
    /* Do we update? */
    if (update_foreground)
@@ -838 +922 @@
       st->Dvar1 = st->Dvar2 = FLOAT_ZERO;
       /* Copy background filter to foreground filter */
-      for (i=0;i<N*M;i++)
+      for (i=0;i<N*M*C*K;i++)
          st->foreground[i] = EXTRACT16(PSHR32(st->W[i],16));
       /* Apply a smooth transition so as to not introduce blocking artifacts */
-      for (i=0;i<st->frame_size;i++)
-         st->e[i+st->frame_size] = 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]);
+      for (chan = 0; chan < C; chan++)
+         for (i=0;i<st->frame_size;i++)
+            st->e[chan*N+i+st->frame_size] = 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]);
    } else {
       int reset_background=0;
@@ -855 +940 @@
       {
          /* Copy foreground filter to background filter */
-         for (i=0;i<N*M;i++)
+         for (i=0;i<N*M*C*K;i++)
             st->W[i] = SHL32(EXTEND32(st->foreground[i]),16);
          /* We also need to copy the output so as to get correct adaptation */
-         for (i=0;i<st->frame_size;i++)
-            st->y[i+st->frame_size] = st->e[i+st->frame_size];
-         for (i=0;i<st->frame_size;i++)
-            st->e[i] = SUB16(st->input[i], st->y[i+st->frame_size]);
+         for (chan = 0; chan < C; chan++)
+         {
+            for (i=0;i<st->frame_size;i++)
+               st->y[chan*N+i+st->frame_size] = st->e[chan*N+i+st->frame_size];
+            for (i=0;i<st->frame_size;i++)
+               st->e[chan*N+i] = SUB16(st->input[chan*st->frame_size+i], st->y[chan*N+i+st->frame_size]);
+         }
          See = Sff;
          st->Davg1 = st->Davg2 = 0;
@@ -869 +957 @@
 #endif
 
-   /* Compute error signal (for the output with de-emphasis) */ 
-   for (i=0;i<st->frame_size;i++)
-   {
-      spx_word32_t tmp_out;
+   Sey = Syy = Sdd = 0;
+   for (chan = 0; chan < C; chan++)
+   {
+      /* Compute error signal (for the output with de-emphasis) */
+      for (i=0;i<st->frame_size;i++)
+      {
+         spx_word32_t tmp_out;
 #ifdef TWO_PATH
-      tmp_out = SUB32(EXTEND32(st->input[i]), EXTEND32(st->e[i+st->frame_size]));
-#else
-      tmp_out = SUB32(EXTEND32(st->input[i]), EXTEND32(st->y[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)));
+         tmp_out = SUB32(EXTEND32(st->input[chan*st->frame_size+i]), EXTEND32(st->e[chan*N+i+st->frame_size]));
+#else
+         tmp_out = SUB32(EXTEND32(st->input[chan*st->frame_size+i]), EXTEND32(st->y[chan*N+i+st->frame_size]));
+#endif
+         tmp_out = ADD32(tmp_out, EXTEND32(MULT16_16_P15(st->preemph, st->memE[chan])));
       /* This is an arbitrary test for saturation in the microphone signal */
-      if (in[i] <= -32000 || in[i] >= 32000)
-      {
-         tmp_out = 0;
+         if (in[i*C+chan] <= -32000 || in[i*C+chan] >= 32000)
+         {
          if (st->saturated == 0)
             st->saturated = 1;
-      }
-      out[i] = (spx_int16_t)tmp_out;
-      st->memE = tmp_out;
-   }
-   
+         }
+         out[i*C+chan] = WORD2INT(tmp_out);
+         st->memE[chan] = tmp_out;
+      }
+
 #ifdef DUMP_ECHO_CANCEL_DATA
-   dump_audio(in, far_end, out, st->frame_size);
-#endif
-   
-   /* Compute error signal (filter update version) */ 
-   for (i=0;i<st->frame_size;i++)
-   {
-      st->e[i+st->frame_size] = st->e[i];
-      st->e[i] = 0;
-   }
-
-   /* Compute a bunch of correlations */
-   Sey = mdf_inner_prod(st->e+st->frame_size, st->y+st->frame_size, st->frame_size);
-   Syy = mdf_inner_prod(st->y+st->frame_size, st->y+st->frame_size, st->frame_size);
-   Sdd = mdf_inner_prod(st->input, st->input, st->frame_size);
-   
+      dump_audio(in, far_end, out, st->frame_size);
+#endif
+
+      /* Compute error signal (filter update version) */
+      for (i=0;i<st->frame_size;i++)
+      {
+         st->e[chan*N+i+st->frame_size] = st->e[chan*N+i];
+         st->e[chan*N+i] = 0;
+      }
+
+      /* Compute a bunch of correlations */
+      /* FIXME: bad merge */
+      Sey += mdf_inner_prod(st->e+chan*N+st->frame_size, st->y+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);
+      Sdd += mdf_inner_prod(st->input+chan*st->frame_size, st->input+chan*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;i<st->frame_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);
+
+   }
+
    /*printf ("%f %f %f %f\n", Sff, See, Syy, Sdd, st->update_cond);*/
-   
+
    /* Do some sanity check */
    if (!(Syy>=0 && Sxx>=0 && See >= 0)
@@ -922 +1020 @@
       /* Things have gone really bad */
       st->screwed_up += 50;
-      for (i=0;i<st->frame_size;i++)
+      for (i=0;i<st->frame_size*C;i++)
          out[i] = 0;
    } else if (SHR32(Sff, 2) > ADD32(Sdd, SHR32(MULT16_16(N, 10000),6)))
@@ -942 +1040 @@
    See = MAX32(See, SHR32(MULT16_16(N, 100),6));
 
-   /* Convert error to frequency domain */
-   spx_fft(st->fft_table, st->e, st->E);
-   for (i=0;i<st->frame_size;i++)
-      st->y[i] = 0;
-   spx_fft(st->fft_table, st->y, st->Y);
-
-   /* Compute power spectrum of far end (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);
+   }
+
+
    /* Smooth far end 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] = 100;
-      for (i=0;i<M;i++)
-      {
-         power_spectrum(&st->X[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 */
@@ -988 +1067 @@
 #endif
    }
-   
+
    Pyy = FLOAT_SQRT(Pyy);
    Pey = FLOAT_DIVU(Pey,Pyy);
@@ -1016 +1095 @@
       st->leak_estimate = SHL16(st->leak_estimate,1);
    /*printf ("%f\n", st->leak_estimate);*/
-   
+
    /* Compute Residual to Error Ratio */
 #ifdef FIXED_POINT
@@ -1072 +1151 @@
       spx_word16_t adapt_rate=0;
 
-      if (Sxx > SHR32(MULT16_16(N, 1000),6)) 
+      if (Sxx > SHR32(MULT16_16(N, 1000),6))
       {
          tmp32 = MULT16_32_Q15(QCONST16(.25f, 15), Sxx);
@@ -1092 +1171 @@
    }
 
-   /* Save residual echo so it can be used by the nonlinear processor */
-   if (st->adapted)
-   {
-      /* If the filter is adapted, take the filtered echo */
+   /* FIXME: MC conversion required */
       for (i=0;i<st->frame_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;i<st->frame_size;i++)
          st->last_y[st->frame_size+i] = in[i]-out[i];
@@ -1114 +1193 @@
    spx_word16_t leak2;
    int N;
-   
+
    N = st->window_size;
 
@@ -1120 +1199 @@
    for (i=0;i<N;i++)
       st->y[i] = MULT16_16_Q15(st->window[i],st->last_y[i]);
-      
+
    /* Compute power spectrum of the echo */
    spx_fft(st->fft_table, st->y, st->Y);
    power_spectrum(st->Y, residual_echo, N);
-      
+
 #ifdef FIXED_POINT
    if (st->leak_estimate > 16383)
@@ -1139 +1218 @@
    for (i=0;i<=st->frame_size;i++)
       residual_echo[i] = (spx_int32_t)MULT16_32_Q15(leak2,residual_echo[i]);
-   
-}
-
-int speex_echo_ctl(SpeexEchoState *st, int request, void *ptr)
+
+}
+
+EXPORT int speex_echo_ctl(SpeexEchoState *st, int request, void *ptr)
 {
    switch(request)
    {
-      
+
       case SPEEX_ECHO_GET_FRAME_SIZE:
          (*(int*)ptr) = st->frame_size;
@@ -1170 +1249 @@
          (*(int*)ptr) = st->sampling_rate;
          break;
+      case SPEEX_ECHO_GET_IMPULSE_RESPONSE_SIZE:
+         /*FIXME: Implement this for multiple channels */
+         *((spx_int32_t *)ptr) = st->M * st->frame_size;
+         break;
+      case SPEEX_ECHO_GET_IMPULSE_RESPONSE:
+      {
+         int M = st->M, N = st->window_size, n = st->frame_size, i, j;
+         spx_int32_t *filt = (spx_int32_t *) ptr;
+         for(j=0;j<M;j++)
+         {
+            /*FIXME: Implement this for multiple channels */
+#ifdef FIXED_POINT
+            for (i=0;i<N;i++)
+               st->wtmp2[i] = EXTRACT16(PSHR32(st->W[j*N+i],16+NORMALIZE_SCALEDOWN));
+            spx_ifft(st->fft_table, st->wtmp2, st->wtmp);
+#else
+            spx_ifft(st->fft_table, &st->W[j*N], st->wtmp);
+#endif
+            for(i=0;i<n;i++)
+               filt[j*n+i] = PSHR32(MULT16_16(32767,st->wtmp[i]), WEIGHT_SHIFT-NORMALIZE_SCALEDOWN);
+         }
+      }
+         break;
       default:
          speex_warning_int("Unknown speex_echo_ctl request: ", request);

pjproject/trunk/third_party/speex/libspeex/misc_bfin.h

@@ -2 +2 @@
 /**
    @file misc_bfin.h
-   @author Jean-Marc Valin 
+   @author Jean-Marc Valin
    @brief Various compatibility routines for Speex (Blackfin version)
 */
@@ -9 +9 @@
    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
@@ -33 +33 @@
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
+
+#include "bfin.h"
 
 #define OVERRIDE_SPEEX_MOVE
@@ -49 +51 @@
    : "=a" (src), "=a" (dest)
    : "a" ((n>>2)-1), "0" (src), "1" (dest)
-   : "R0", "I0", "L0", "memory"
+   : "R0", "I0", "L0", "memory" BFIN_HWLOOP0_REGS
          );
    return dest;

pjproject/trunk/third_party/speex/libspeex/os_support.h

@@ -1 +1 @@
 /* Copyright (C) 2007 Jean-Marc Valin
-      
+
    File: os_support.h
    This is the (tiny) OS abstraction layer. Aside from math.h, this is the
@@ -46 +46 @@
 #endif
 
-/** Speex wrapper for calloc. To do your own dynamic allocation, all you need to do is replace this function, speex_realloc and speex_free 
+/** Speex wrapper for calloc. To do your own dynamic allocation, all you need to do is replace this function, speex_realloc and speex_free
     NOTE: speex_alloc needs to CLEAR THE MEMORY */
 #ifndef OVERRIDE_SPEEX_ALLOC
 static inline void *speex_alloc (int size)
 {
-   /* WARNING: this is not equivalent to malloc(). If you want to use malloc() 
+   /* WARNING: this is not equivalent to malloc(). If you want to use malloc()
       or your own allocator, YOU NEED TO CLEAR THE MEMORY ALLOCATED. Otherwise
       you will experience strange bugs */
@@ -91 +91 @@
 #endif
 
-/** Copy n bytes of memory from src to dst. The 0* term provides compile-time type checking  */
+/** Copy n elements from src to dst. The 0* term provides compile-time type checking  */
 #ifndef OVERRIDE_SPEEX_COPY
 #define SPEEX_COPY(dst, src, n) (memcpy((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) ))
 #endif
 
-/** Copy n bytes of memory from src to dst, allowing overlapping regions. The 0* term 
+/** Copy n elements from src to dst, allowing overlapping regions. The 0* term
     provides compile-time type checking */
 #ifndef OVERRIDE_SPEEX_MOVE
@@ -102 +102 @@
 #endif
 
-/** Set n bytes of memory to value of c, starting at address s */
+/** For n elements worth of memory, set every byte to the value of c, starting at address dst */
 #ifndef OVERRIDE_SPEEX_MEMSET
 #define SPEEX_MEMSET(dst, c, n) (memset((dst), (c), (n)*sizeof(*(dst))))

pjproject/trunk/third_party/speex/libspeex/resample.c

@@ -1 +1 @@
 /* Copyright (C) 2007-2008 Jean-Marc Valin
    Copyright (C) 2008      Thorvald Natvig
-      
+
    File: resample.c
    Arbitrary resampling code
@@ -39 +39 @@
       - Good *perceptual* quality (and not best SNR)
 
-   Warning: This resampler is relatively new. Although I think I got rid of 
+   Warning: This resampler is relatively new. Although I think I got rid of
    all the major bugs and I don't expect the API to change anymore, there
    may be something I've missed. So use with caution.
@@ -45 +45 @@
    This algorithm is based on this original resampling algorithm:
    Smith, Julius O. Digital Audio Resampling Home Page
-   Center for Computer Research in Music and Acoustics (CCRMA), 
+   Center for Computer Research in Music and Acoustics (CCRMA),
    Stanford University, 2007.
-   Web published at http://www-ccrma.stanford.edu/~jos/resample/.
-
-   There is one main difference, though. This resampler uses cubic 
+   Web published at https://ccrma.stanford.edu/~jos/resample/.
+
+   There is one main difference, though. This resampler uses cubic
    interpolation instead of linear interpolation in the above paper. This
    makes the table much smaller and makes it possible to compute that table
-   on a per-stream basis. In turn, being able to tweak the table for each 
-   stream makes it possible to both reduce complexity on simple ratios 
-   (e.g. 2/3), and get rid of the rounding operations in the inner loop. 
+   on a per-stream basis. In turn, being able to tweak the table for each
+   stream makes it possible to both reduce complexity on simple ratios
+   (e.g. 2/3), and get rid of the rounding operations in the inner loop.
    The latter both reduces CPU time and makes the algorithm more SIMD-friendly.
 */
@@ -64 +64 @@
 #ifdef OUTSIDE_SPEEX
 #include <stdlib.h>
-static void *speex_alloc (int size) {return calloc(size,1);}
-static void *speex_realloc (void *ptr, int size) {return realloc(ptr, size);}
-static void speex_free (void *ptr) {free(ptr);}
+static void *speex_alloc(int size) {return calloc(size,1);}
+static void *speex_realloc(void *ptr, int size) {return realloc(ptr, size);}
+static void speex_free(void *ptr) {free(ptr);}
+#ifndef EXPORT
+#define EXPORT
+#endif
 #include "speex_resampler.h"
 #include "arch.h"
 #else /* OUTSIDE_SPEEX */
-               
+
 #include "speex/speex_resampler.h"
 #include "arch.h"
@@ -76 +79 @@
 #endif /* OUTSIDE_SPEEX */
 
-#include "stack_alloc.h"
 #include <math.h>
+#include <limits.h>
 
 #ifndef M_PI
-#define M_PI 3.14159263
-#endif
-
-#ifdef FIXED_POINT
-#define WORD2INT(x) ((x) < -32767 ? -32768 : ((x) > 32766 ? 32767 : (x)))  
-#else
-#define WORD2INT(x) ((x) < -32767.5f ? -32768 : ((x) > 32766.5f ? 32767 : floor(.5+(x))))  
-#endif
-               
+#define M_PI 3.14159265358979323846
+#endif
+
 #define IMAX(a,b) ((a) > (b) ? (a) : (b))
 #define IMIN(a,b) ((a) < (b) ? (a) : (b))
@@ -96 +93 @@
 #endif
 
-#ifdef _USE_SSE
+#ifndef UINT32_MAX
+#define UINT32_MAX 4294967295U
+#endif
+
+#ifdef USE_SSE
 #include "resample_sse.h"
+#endif
+
+#ifdef USE_NEON
+#include "resample_neon.h"
 #endif
 
@@ -114 +119 @@
    spx_uint32_t num_rate;
    spx_uint32_t den_rate;
-   
+
    int    quality;
    spx_uint32_t nb_channels;
@@ -126 +131 @@
    int          initialised;
    int          started;
-   
+
    /* These are per-channel */
    spx_int32_t  *last_sample;
    spx_uint32_t *samp_frac_num;
    spx_uint32_t *magic_samples;
-   
+
    spx_word16_t *mem;
    spx_word16_t *sinc_table;
    spx_uint32_t sinc_table_length;
    resampler_basic_func resampler_ptr;
-         
+
    int    in_stride;
    int    out_stride;
 } ;
 
-static double kaiser12_table[68] = {
+static const double kaiser12_table[68] = {
    0.99859849, 1.00000000, 0.99859849, 0.99440475, 0.98745105, 0.97779076,
    0.96549770, 0.95066529, 0.93340547, 0.91384741, 0.89213598, 0.86843014,
@@ -155 +160 @@
    0.00001000, 0.00000000};
 /*
-static double kaiser12_table[36] = {
+static const double kaiser12_table[36] = {
    0.99440475, 1.00000000, 0.99440475, 0.97779076, 0.95066529, 0.91384741,
    0.86843014, 0.81573067, 0.75723148, 0.69451601, 0.62920216, 0.56287762,
@@ -163 +168 @@
    0.00153438, 0.00069463, 0.00025272, 0.0000527734, 0.00000500, 0.00000000};
 */
-static double kaiser10_table[36] = {
+static const double kaiser10_table[36] = {
    0.99537781, 1.00000000, 0.99537781, 0.98162644, 0.95908712, 0.92831446,
    0.89005583, 0.84522401, 0.79486424, 0.74011713, 0.68217934, 0.62226347,
@@ -171 +176 @@
    0.00488951, 0.00257636, 0.00115101, 0.00035515, 0.00000000, 0.00000000};
 
-static double kaiser8_table[36] = {
+static const double kaiser8_table[36] = {
    0.99635258, 1.00000000, 0.99635258, 0.98548012, 0.96759014, 0.94302200,
    0.91223751, 0.87580811, 0.83439927, 0.78875245, 0.73966538, 0.68797126,
@@ -178 +183 @@
    0.10562887, 0.08273982, 0.06335451, 0.04724088, 0.03412321, 0.02369490,
    0.01563093, 0.00959968, 0.00527363, 0.00233883, 0.00050000, 0.00000000};
-   
-static double kaiser6_table[36] = {
+
+static const double kaiser6_table[36] = {
    0.99733006, 1.00000000, 0.99733006, 0.98935595, 0.97618418, 0.95799003,
    0.93501423, 0.90755855, 0.87598009, 0.84068475, 0.80211977, 0.76076565,
@@ -188 +193 @@
 
 struct FuncDef {
-   double *table;
+   const double *table;
    int oversample;
 };
-      
-static struct FuncDef _KAISER12 = {kaiser12_table, 64};
-#define KAISER12 (&_KAISER12)
-/*static struct FuncDef _KAISER12 = {kaiser12_table, 32};
-#define KAISER12 (&_KAISER12)*/
-static struct FuncDef _KAISER10 = {kaiser10_table, 32};
-#define KAISER10 (&_KAISER10)
-static struct FuncDef _KAISER8 = {kaiser8_table, 32};
-#define KAISER8 (&_KAISER8)
-static struct FuncDef _KAISER6 = {kaiser6_table, 32};
-#define KAISER6 (&_KAISER6)
+
+static const struct FuncDef kaiser12_funcdef = {kaiser12_table, 64};
+#define KAISER12 (&kaiser12_funcdef)
+static const struct FuncDef kaiser10_funcdef = {kaiser10_table, 32};
+#define KAISER10 (&kaiser10_funcdef)
+static const struct FuncDef kaiser8_funcdef = {kaiser8_table, 32};
+#define KAISER8 (&kaiser8_funcdef)
+static const struct FuncDef kaiser6_funcdef = {kaiser6_table, 32};
+#define KAISER6 (&kaiser6_funcdef)
 
 struct QualityMapping {
@@ -208 +211 @@
    float downsample_bandwidth;
    float upsample_bandwidth;
-   struct FuncDef *window_func;
+   const struct FuncDef *window_func;
 };
 
 
 /* This table maps conversion quality to internal parameters. There are two
-   reasons that explain why the up-sampling bandwidth is larger than the 
+   reasons that explain why the up-sampling bandwidth is larger than the
    down-sampling bandwidth:
    1) When up-sampling, we can assume that the spectrum is already attenuated
@@ -235 +238 @@
 };
 /*8,24,40,56,80,104,128,160,200,256,320*/
-static double compute_func(float x, struct FuncDef *func)
+static double compute_func(float x, const struct FuncDef *func)
 {
    float y, frac;
    double interp[4];
-   int ind; 
+   int ind;
    y = x*func->oversample;
    ind = (int)floor(y);
@@ -250 +253 @@
    /* Just to make sure we don't have rounding problems */
    interp[1] = 1.f-interp[3]-interp[2]-interp[0];
-   
+
    /*sum = frac*accum[1] + (1-frac)*accum[2];*/
    return interp[0]*func->table[ind] + interp[1]*func->table[ind+1] + interp[2]*func->table[ind+2] + interp[3]*func->table[ind+3];
@@ -270 +273 @@
 #ifdef FIXED_POINT
 /* The slow way of computing a sinc for the table. Should improve that some day */
-static spx_word16_t sinc(float cutoff, float x, int N, struct FuncDef *window_func)
+static spx_word16_t sinc(float cutoff, float x, int N, const struct FuncDef *window_func)
 {
    /*fprintf (stderr, "%f ", x);*/
@@ -283 +286 @@
 #else
 /* The slow way of computing a sinc for the table. Should improve that some day */
-static spx_word16_t sinc(float cutoff, float x, int N, struct FuncDef *window_func)
+static spx_word16_t sinc(float cutoff, float x, int N, const struct FuncDef *window_func)
 {
    /*fprintf (stderr, "%f ", x);*/
@@ -338 +341 @@
    const spx_uint32_t den_rate = st->den_rate;
    spx_word32_t sum;
-   int j;
 
    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
    {
-      const spx_word16_t *sinc = & sinc_table[samp_frac_num*N];
+      const spx_word16_t *sinct = & sinc_table[samp_frac_num*N];
       const spx_word16_t *iptr = & in[last_sample];
 
 #ifndef OVERRIDE_INNER_PRODUCT_SINGLE
-      float accum[4] = {0,0,0,0};
-
+      int j;
+      sum = 0;
+      for(j=0;j<N;j++) sum += MULT16_16(sinct[j], iptr[j]);
+
+/*    This code is slower on most DSPs which have only 2 accumulators.
+      Plus this this forces truncation to 32 bits and you lose the HW guard bits.
+      I think we can trust the compiler and let it vectorize and/or unroll itself.
+      spx_word32_t accum[4] = {0,0,0,0};
       for(j=0;j<N;j+=4) {
-        accum[0] += sinc[j]*iptr[j];
-        accum[1] += sinc[j+1]*iptr[j+1];
-        accum[2] += sinc[j+2]*iptr[j+2];
-        accum[3] += sinc[j+3]*iptr[j+3];
+        accum[0] += MULT16_16(sinct[j], iptr[j]);
+        accum[1] += MULT16_16(sinct[j+1], iptr[j+1]);
+        accum[2] += MULT16_16(sinct[j+2], iptr[j+2]);
+        accum[3] += MULT16_16(sinct[j+3], iptr[j+3]);
       }
       sum = accum[0] + accum[1] + accum[2] + accum[3];
-#else
-      sum = inner_product_single(sinc, iptr, N);
-#endif
-
-      out[out_stride * out_sample++] = PSHR32(sum, 15);
+*/
+      sum = SATURATE32PSHR(sum, 15, 32767);
+#else
+      sum = inner_product_single(sinct, iptr, N);
+#endif
+
+      out[out_stride * out_sample++] = sum;
       last_sample += int_advance;
       samp_frac_num += frac_advance;
@@ -389 +399 @@
    const spx_uint32_t den_rate = st->den_rate;
    double sum;
-   int j;
 
    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
    {
-      const spx_word16_t *sinc = & sinc_table[samp_frac_num*N];
+      const spx_word16_t *sinct = & sinc_table[samp_frac_num*N];
       const spx_word16_t *iptr = & in[last_sample];
 
 #ifndef OVERRIDE_INNER_PRODUCT_DOUBLE
+      int j;
       double accum[4] = {0,0,0,0};
 
       for(j=0;j<N;j+=4) {
-        accum[0] += sinc[j]*iptr[j];
-        accum[1] += sinc[j+1]*iptr[j+1];
-        accum[2] += sinc[j+2]*iptr[j+2];
-        accum[3] += sinc[j+3]*iptr[j+3];
+        accum[0] += sinct[j]*iptr[j];
+        accum[1] += sinct[j+1]*iptr[j+1];
+        accum[2] += sinct[j+2]*iptr[j+2];
+        accum[3] += sinct[j+3]*iptr[j+3];
       }
       sum = accum[0] + accum[1] + accum[2] + accum[3];
 #else
-      sum = inner_product_double(sinc, iptr, N);
+      sum = inner_product_double(sinct, iptr, N);
 #endif
 
@@ -436 +446 @@
    const int frac_advance = st->frac_advance;
    const spx_uint32_t den_rate = st->den_rate;
-   int j;
    spx_word32_t sum;
 
@@ -453 +462 @@
 
 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_SINGLE
+      int j;
       spx_word32_t accum[4] = {0,0,0,0};
 
@@ -464 +474 @@
 
       cubic_coef(frac, interp);
-      sum = MULT16_32_Q15(interp[0],accum[0]) + MULT16_32_Q15(interp[1],accum[1]) + MULT16_32_Q15(interp[2],accum[2]) + MULT16_32_Q15(interp[3],accum[3]);
+      sum = MULT16_32_Q15(interp[0],SHR32(accum[0], 1)) + MULT16_32_Q15(interp[1],SHR32(accum[1], 1)) + MULT16_32_Q15(interp[2],SHR32(accum[2], 1)) + MULT16_32_Q15(interp[3],SHR32(accum[3], 1));
+      sum = SATURATE32PSHR(sum, 15, 32767);
 #else
       cubic_coef(frac, interp);
       sum = interpolate_product_single(iptr, st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample, interp);
 #endif
-      
-      out[out_stride * out_sample++] = PSHR32(sum,15);
+
+      out[out_stride * out_sample++] = sum;
       last_sample += int_advance;
       samp_frac_num += frac_advance;
@@ -498 +509 @@
    const int frac_advance = st->frac_advance;
    const spx_uint32_t den_rate = st->den_rate;
-   int j;
    spx_word32_t sum;
 
@@ -515 +525 @@
 
 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE
+      int j;
       double accum[4] = {0,0,0,0};
 
@@ -531 +542 @@
       sum = interpolate_product_double(iptr, st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample, interp);
 #endif
-      
+
       out[out_stride * out_sample++] = PSHR32(sum,15);
       last_sample += int_advance;
@@ -548 +559 @@
 #endif
 
-static void update_filter(SpeexResamplerState *st)
-{
-   spx_uint32_t old_length;
-   
-   old_length = st->filt_len;
+/* This resampler is used to produce zero output in situations where memory
+   for the filter could not be allocated.  The expected numbers of input and
+   output samples are still processed so that callers failing to check error
+   codes are not surprised, possibly getting into infinite loops. */
+static int resampler_basic_zero(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_word16_t *in, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
+{
+   int out_sample = 0;
+   int last_sample = st->last_sample[channel_index];
+   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
+   const int out_stride = st->out_stride;
+   const int int_advance = st->int_advance;
+   const int frac_advance = st->frac_advance;
+   const spx_uint32_t den_rate = st->den_rate;
+
+   (void)in;
+   while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
+   {
+      out[out_stride * out_sample++] = 0;
+      last_sample += int_advance;
+      samp_frac_num += frac_advance;
+      if (samp_frac_num >= den_rate)
+      {
+         samp_frac_num -= den_rate;
+         last_sample++;
+      }
+   }
+
+   st->last_sample[channel_index] = last_sample;
+   st->samp_frac_num[channel_index] = samp_frac_num;
+   return out_sample;
+}
+
+static int multiply_frac(spx_uint32_t *result, spx_uint32_t value, spx_uint32_t num, spx_uint32_t den)
+{
+   spx_uint32_t major = value / den;
+   spx_uint32_t remain = value % den;
+   /* TODO: Could use 64 bits operation to check for overflow. But only guaranteed in C99+ */
+   if (remain > UINT32_MAX / num || major > UINT32_MAX / num
+       || major * num > UINT32_MAX - remain * num / den)
+      return RESAMPLER_ERR_OVERFLOW;
+   *result = remain * num / den + major * num;
+   return RESAMPLER_ERR_SUCCESS;
+}
+
+static int update_filter(SpeexResamplerState *st)
+{
+   spx_uint32_t old_length = st->filt_len;
+   spx_uint32_t old_alloc_size = st->mem_alloc_size;
+   int use_direct;
+   spx_uint32_t min_sinc_table_length;
+   spx_uint32_t min_alloc_size;
+
+   st->int_advance = st->num_rate/st->den_rate;
+   st->frac_advance = st->num_rate%st->den_rate;
    st->oversample = quality_map[st->quality].oversample;
    st->filt_len = quality_map[st->quality].base_length;
-   
+
    if (st->num_rate > st->den_rate)
    {
       /* down-sampling */
       st->cutoff = quality_map[st->quality].downsample_bandwidth * st->den_rate / st->num_rate;
-      /* FIXME: divide the numerator and denominator by a certain amount if they're too large */
-      st->filt_len = st->filt_len*st->num_rate / st->den_rate;
-      /* Round down to make sure we have a multiple of 4 */
-      st->filt_len &= (~0x3);
+      if (multiply_frac(&st->filt_len,st->filt_len,st->num_rate,st->den_rate) != RESAMPLER_ERR_SUCCESS)
+         goto fail;
+      /* Round up to make sure we have a multiple of 8 for SSE */
+      st->filt_len = ((st->filt_len-1)&(~0x7))+8;
       if (2*st->den_rate < st->num_rate)
          st->oversample >>= 1;
@@ -578 +638 @@
       st->cutoff = quality_map[st->quality].upsample_bandwidth;
    }
-   
+
+#ifdef RESAMPLE_FULL_SINC_TABLE
+   use_direct = 1;
+   if (INT_MAX/sizeof(spx_word16_t)/st->den_rate < st->filt_len)
+      goto fail;
+#else
    /* Choose the resampling type that requires the least amount of memory */
-   if (st->den_rate <= st->oversample)
+   use_direct = st->filt_len*st->den_rate <= st->filt_len*st->oversample+8
+                && INT_MAX/sizeof(spx_word16_t)/st->den_rate >= st->filt_len;
+#endif
+   if (use_direct)
+   {
+      min_sinc_table_length = st->filt_len*st->den_rate;
+   } else {
+      if ((INT_MAX/sizeof(spx_word16_t)-8)/st->oversample < st->filt_len)
+         goto fail;
+
+      min_sinc_table_length = st->filt_len*st->oversample+8;
+   }
+   if (st->sinc_table_length < min_sinc_table_length)
+   {
+      spx_word16_t *sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,min_sinc_table_length*sizeof(spx_word16_t));
+      if (!sinc_table)
+         goto fail;
+
+      st->sinc_table = sinc_table;
+      st->sinc_table_length = min_sinc_table_length;
+   }
+   if (use_direct)
    {
       spx_uint32_t i;
-      if (!st->sinc_table)
-         st->sinc_table = (spx_word16_t *)speex_alloc(st->filt_len*st->den_rate*sizeof(spx_word16_t));
-      else if (st->sinc_table_length < st->filt_len*st->den_rate)
-      {
-         st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,st->filt_len*st->den_rate*sizeof(spx_word16_t));
-         st->sinc_table_length = st->filt_len*st->den_rate;
-      }
       for (i=0;i<st->den_rate;i++)
       {
@@ -609 +688 @@
    } else {
       spx_int32_t i;
-      if (!st->sinc_table)
-         st->sinc_table = (spx_word16_t *)speex_alloc((st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
-      else if (st->sinc_table_length < st->filt_len*st->oversample+8)
-      {
-         st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,(st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
-         st->sinc_table_length = st->filt_len*st->oversample+8;
-      }
       for (i=-4;i<(spx_int32_t)(st->oversample*st->filt_len+4);i++)
          st->sinc_table[i+4] = sinc(st->cutoff,(i/(float)st->oversample - st->filt_len/2), st->filt_len, quality_map[st->quality].window_func);
@@ -628 +700 @@
       /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff);*/
    }
-   st->int_advance = st->num_rate/st->den_rate;
-   st->frac_advance = st->num_rate%st->den_rate;
-
-   
+
    /* Here's the place where we update the filter memory to take into account
       the change in filter length. It's probably the messiest part of the code
       due to handling of lots of corner cases. */
-   if (!st->mem)
+
+   /* Adding buffer_size to filt_len won't overflow here because filt_len
+      could be multiplied by sizeof(spx_word16_t) above. */
+   min_alloc_size = st->filt_len-1 + st->buffer_size;
+   if (min_alloc_size > st->mem_alloc_size)
+   {
+      spx_word16_t *mem;
+      if (INT_MAX/sizeof(spx_word16_t)/st->nb_channels < min_alloc_size)
+          goto fail;
+      else if (!(mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*min_alloc_size * sizeof(*mem))))
+          goto fail;
+
+      st->mem = mem;
+      st->mem_alloc_size = min_alloc_size;
+   }
+   if (!st->started)
    {
       spx_uint32_t i;
-      st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
-      st->mem = (spx_word16_t*)speex_alloc(st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
-      for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
-         st->mem[i] = 0;
-      /*speex_warning("init filter");*/
-   } else if (!st->started)
-   {
-      spx_uint32_t i;
-      st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
-      st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
          st->mem[i] = 0;
@@ -653 +727 @@
    } else if (st->filt_len > old_length)
    {
-      spx_int32_t i;
+      spx_uint32_t i;
       /* Increase the filter length */
       /*speex_warning("increase filter size");*/
-      int old_alloc_size = st->mem_alloc_size;
-      if ((st->filt_len-1 + st->buffer_size) > st->mem_alloc_size)
+      for (i=st->nb_channels;i--;)
       {
-         st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
-         st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
-      }
-      for (i=st->nb_channels-1;i>=0;i--)
-      {
-         spx_int32_t j;
+         spx_uint32_t j;
          spx_uint32_t olen = old_length;
          /*if (st->magic_samples[i])*/
          {
             /* Try and remove the magic samples as if nothing had happened */
-            
+
             /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */
             olen = old_length + 2*st->magic_samples[i];
-            for (j=old_length-2+st->magic_samples[i];j>=0;j--)
+            for (j=old_length-1+st->magic_samples[i];j--;)
                st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]] = st->mem[i*old_alloc_size+j];
             for (j=0;j<st->magic_samples[i];j++)
@@ -713 +781 @@
       }
    }
-
+   return RESAMPLER_ERR_SUCCESS;
+
+fail:
+   st->resampler_ptr = resampler_basic_zero;
+   /* st->mem may still contain consumed input samples for the filter.
+      Restore filt_len so that filt_len - 1 still points to the position after
+      the last of these samples. */
+   st->filt_len = old_length;
+   return RESAMPLER_ERR_ALLOC_FAILED;
 }
 
@@ -723 +799 @@
 EXPORT SpeexResamplerState *speex_resampler_init_frac(spx_uint32_t nb_channels, spx_uint32_t ratio_num, spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
 {
-   spx_uint32_t i;
    SpeexResamplerState *st;
-   if (quality > 10 || quality < 0)
+   int filter_err;
+
+   if (nb_channels == 0 || ratio_num == 0 || ratio_den == 0 || quality > 10 || quality < 0)
    {
       if (err)
@@ -732 +809 @@
    }
    st = (SpeexResamplerState *)speex_alloc(sizeof(SpeexResamplerState));
+   if (!st)
+   {
+      if (err)
+         *err = RESAMPLER_ERR_ALLOC_FAILED;
+      return NULL;
+   }
    st->initialised = 0;
    st->started = 0;
@@ -744 +827 @@
    st->mem = 0;
    st->resampler_ptr = 0;
-         
+
    st->cutoff = 1.f;
    st->nb_channels = nb_channels;
    st->in_stride = 1;
    st->out_stride = 1;
-   
-#ifdef FIXED_POINT
+
    st->buffer_size = 160;
-#else
-   st->buffer_size = 160;
-#endif
-   
+
    /* Per channel data */
-   st->last_sample = (spx_int32_t*)speex_alloc(nb_channels*sizeof(int));
-   st->magic_samples = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(int));
-   st->samp_frac_num = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(int));
-   for (i=0;i<nb_channels;i++)
-   {
-      st->last_sample[i] = 0;
-      st->magic_samples[i] = 0;
-      st->samp_frac_num[i] = 0;
-   }
+   if (!(st->last_sample = (spx_int32_t*)speex_alloc(nb_channels*sizeof(spx_int32_t))))
+      goto fail;
+   if (!(st->magic_samples = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(spx_uint32_t))))
+      goto fail;
+   if (!(st->samp_frac_num = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(spx_uint32_t))))
+      goto fail;
 
    speex_resampler_set_quality(st, quality);
    speex_resampler_set_rate_frac(st, ratio_num, ratio_den, in_rate, out_rate);
 
-   
-   update_filter(st);
-   
-   st->initialised = 1;
+   filter_err = update_filter(st);
+   if (filter_err == RESAMPLER_ERR_SUCCESS)
+   {
+      st->initialised = 1;
+   } else {
+      speex_resampler_destroy(st);
+      st = NULL;
+   }
    if (err)
-      *err = RESAMPLER_ERR_SUCCESS;
+      *err = filter_err;
 
    return st;
+
+fail:
+   if (err)
+      *err = RESAMPLER_ERR_ALLOC_FAILED;
+   speex_resampler_destroy(st);
+   return NULL;
 }
 
@@ -797 +883 @@
    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
    spx_uint32_t ilen;
-   
+
    st->started = 1;
-   
+
    /* Call the right resampler through the function ptr */
    out_sample = st->resampler_ptr(st, channel_index, mem, in_len, out, out_len);
-   
+
    if (st->last_sample[channel_index] < (spx_int32_t)*in_len)
       *in_len = st->last_sample[channel_index];
    *out_len = out_sample;
    st->last_sample[channel_index] -= *in_len;
-   
+
    ilen = *in_len;
 
@@ -820 +906 @@
    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
    const int N = st->filt_len;
-   
+
    speex_resampler_process_native(st, channel_index, &tmp_in_len, *out, &out_len);
 
    st->magic_samples[channel_index] -= tmp_in_len;
-   
+
    /* If we couldn't process all "magic" input samples, save the rest for next time */
    if (st->magic_samples[channel_index])
@@ -850 +936 @@
    const int istride = st->in_stride;
 
-   if (st->magic_samples[channel_index]) 
+   if (st->magic_samples[channel_index])
       olen -= speex_resampler_magic(st, channel_index, &out, olen);
    if (! st->magic_samples[channel_index]) {
@@ -856 +942 @@
         spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
         spx_uint32_t ochunk = olen;
- 
+
         if (in) {
            for(j=0;j<ichunk;++j)
@@ -874 +960 @@
    *in_len -= ilen;
    *out_len -= olen;
-   return RESAMPLER_ERR_SUCCESS;
+   return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
 }
 
@@ -892 +978 @@
 #ifdef VAR_ARRAYS
    const unsigned int ylen = (olen < FIXED_STACK_ALLOC) ? olen : FIXED_STACK_ALLOC;
-   VARDECL(spx_word16_t *ystack);
-   ALLOC(ystack, ylen, spx_word16_t);
+   spx_word16_t ystack[ylen];
 #else
    const unsigned int ylen = FIXED_STACK_ALLOC;
@@ -900 +985 @@
 
    st->out_stride = 1;
-   
+
    while (ilen && olen) {
      spx_word16_t *y = ystack;
@@ -937 +1022 @@
         out[j*ostride_save] = WORD2INT(ystack[j]);
 #endif
-     
+
      ilen -= ichunk;
      olen -= ochunk;
@@ -948 +1033 @@
    *out_len -= olen;
 
-   return RESAMPLER_ERR_SUCCESS;
+   return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
 }
 
@@ -955 +1040 @@
    spx_uint32_t i;
    int istride_save, ostride_save;
-   spx_uint32_t bak_len = *out_len;
+   spx_uint32_t bak_out_len = *out_len;
+   spx_uint32_t bak_in_len = *in_len;
    istride_save = st->in_stride;
    ostride_save = st->out_stride;
@@ -961 +1047 @@
    for (i=0;i<st->nb_channels;i++)
    {
-      *out_len = bak_len;
+      *out_len = bak_out_len;
+      *in_len = bak_in_len;
       if (in != NULL)
          speex_resampler_process_float(st, i, in+i, in_len, out+i, out_len);
@@ -969 +1056 @@
    st->in_stride = istride_save;
    st->out_stride = ostride_save;
-   return RESAMPLER_ERR_SUCCESS;
-}
-               
+   return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
+}
+
 EXPORT int speex_resampler_process_interleaved_int(SpeexResamplerState *st, const spx_int16_t *in, spx_uint32_t *in_len, spx_int16_t *out, spx_uint32_t *out_len)
 {
    spx_uint32_t i;
    int istride_save, ostride_save;
-   spx_uint32_t bak_len = *out_len;
+   spx_uint32_t bak_out_len = *out_len;
+   spx_uint32_t bak_in_len = *in_len;
    istride_save = st->in_stride;
    ostride_save = st->out_stride;
@@ -982 +1070 @@
    for (i=0;i<st->nb_channels;i++)
    {
-      *out_len = bak_len;
+      *out_len = bak_out_len;
+      *in_len = bak_in_len;
       if (in != NULL)
          speex_resampler_process_int(st, i, in+i, in_len, out+i, out_len);
@@ -990 +1079 @@
    st->in_stride = istride_save;
    st->out_stride = ostride_save;
-   return RESAMPLER_ERR_SUCCESS;
+   return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
 }
 
@@ -1002 +1091 @@
    *in_rate = st->in_rate;
    *out_rate = st->out_rate;
+}
+
+static inline spx_uint32_t compute_gcd(spx_uint32_t a, spx_uint32_t b)
+{
+   while (b != 0)
+   {
+      spx_uint32_t temp = a;
+
+      a = b;
+      b = temp % b;
+   }
+   return a;
 }
 
@@ -1009 +1110 @@
    spx_uint32_t old_den;
    spx_uint32_t i;
+
+   if (ratio_num == 0 || ratio_den == 0)
+      return RESAMPLER_ERR_INVALID_ARG;
+
    if (st->in_rate == in_rate && st->out_rate == out_rate && st->num_rate == ratio_num && st->den_rate == ratio_den)
       return RESAMPLER_ERR_SUCCESS;
-   
+
    old_den = st->den_rate;
    st->in_rate = in_rate;
@@ -1017 +1122 @@
    st->num_rate = ratio_num;
    st->den_rate = ratio_den;
-   /* FIXME: This is terribly inefficient, but who cares (at least for now)? */
-   for (fact=2;fact<=IMIN(st->num_rate, st->den_rate);fact++)
-   {
-      while ((st->num_rate % fact == 0) && (st->den_rate % fact == 0))
-      {
-         st->num_rate /= fact;
-         st->den_rate /= fact;
-      }
-   }
-      
+
+   fact = compute_gcd(st->num_rate, st->den_rate);
+
+   st->num_rate /= fact;
+   st->den_rate /= fact;
+
    if (old_den > 0)
    {
       for (i=0;i<st->nb_channels;i++)
       {
-         st->samp_frac_num[i]=st->samp_frac_num[i]*st->den_rate/old_den;
+         if (multiply_frac(&st->samp_frac_num[i],st->samp_frac_num[i],st->den_rate,old_den) != RESAMPLER_ERR_SUCCESS)
+            return RESAMPLER_ERR_OVERFLOW;
          /* Safety net */
          if (st->samp_frac_num[i] >= st->den_rate)
@@ -1037 +1139 @@
       }
    }
-   
+
    if (st->initialised)
-      update_filter(st);
+      return update_filter(st);
    return RESAMPLER_ERR_SUCCESS;
 }
@@ -1057 +1159 @@
    st->quality = quality;
    if (st->initialised)
-      update_filter(st);
+      return update_filter(st);
    return RESAMPLER_ERR_SUCCESS;
 }
@@ -1107 +1209 @@
 {
    spx_uint32_t i;
+   for (i=0;i<st->nb_channels;i++)
+   {
+      st->last_sample[i] = 0;
+      st->magic_samples[i] = 0;
+      st->samp_frac_num[i] = 0;
+   }
    for (i=0;i<st->nb_channels*(st->filt_len-1);i++)
       st->mem[i] = 0;

pjproject/trunk/third_party/speex/libspeex/resample_sse.h

@@ -10 +10 @@
    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
@@ -72 +72 @@
 }
 
-#ifdef _USE_SSE2
+#ifdef USE_SSE2
 #include <emmintrin.h>
 #define OVERRIDE_INNER_PRODUCT_DOUBLE
@@ -92 +92 @@
       sum = _mm_add_pd(sum, _mm_cvtps_pd(_mm_movehl_ps(t, t)));
    }
-   sum = _mm_add_sd(sum, (__m128d) _mm_movehl_ps((__m128) sum, (__m128) sum));
+   sum = _mm_add_sd(sum, _mm_unpackhi_pd(sum, sum));
    _mm_store_sd(&ret, sum);
    return ret;
@@ -121 +121 @@
   sum2 = _mm_mul_pd(f2, sum2);
   sum = _mm_add_pd(sum1, sum2);
-  sum = _mm_add_sd(sum, (__m128d) _mm_movehl_ps((__m128) sum, (__m128) sum));
+  sum = _mm_add_sd(sum, _mm_unpackhi_pd(sum, sum));
   _mm_store_sd(&ret, sum);
   return ret;

pjproject/trunk/third_party/speex/libspeex/scal.c

@@ -52 +52 @@
 #include <math.h>
 #include <stdlib.h>
+
+#ifndef M_PI
+#define M_PI           3.14159265358979323846  /* pi */
+#endif
 
 #define ALLPASS_ORDER 20
@@ -172 +176 @@
 
       x = buff+st->frame_size;
-      beta = 1.-.3*amount*amount;
       if (amount>1)
          beta = 1-sqrt(.4*amount);