From 67d4e56dc8a365871c3dca4f04fcf8b9c9f47ee6 Mon Sep 17 00:00:00 2001 From: Benny Prijono Date: Thu, 23 Nov 2006 10:19:46 +0000 Subject: Updated Speex to their latest SVN (1.2-beta). AEC seems to work much better now and take less CPU, so I increased default tail length in PJSUA to 800ms. git-svn-id: http://svn.pjsip.org/repos/pjproject/trunk@823 74dad513-b988-da41-8d7b-12977e46ad98 --- pjmedia/src/pjmedia-codec/speex/VERSION.TXT | 9 + pjmedia/src/pjmedia-codec/speex/_kiss_fft_guts.h | 6 + pjmedia/src/pjmedia-codec/speex/arch.h | 6 + pjmedia/src/pjmedia-codec/speex/cb_search.c | 28 +- pjmedia/src/pjmedia-codec/speex/cb_search.h | 6 +- pjmedia/src/pjmedia-codec/speex/cb_search_bfin.h | 5 +- pjmedia/src/pjmedia-codec/speex/fftwrap.c | 9 +- pjmedia/src/pjmedia-codec/speex/filterbank.c | 200 ++++ pjmedia/src/pjmedia-codec/speex/filterbank.h | 66 + pjmedia/src/pjmedia-codec/speex/filters.c | 94 +- pjmedia/src/pjmedia-codec/speex/filters.h | 9 + pjmedia/src/pjmedia-codec/speex/fixed_debug.h | 47 +- pjmedia/src/pjmedia-codec/speex/fixed_generic.h | 8 +- pjmedia/src/pjmedia-codec/speex/jitter.c | 4 +- pjmedia/src/pjmedia-codec/speex/kiss_fft.c | 26 +- pjmedia/src/pjmedia-codec/speex/kiss_fftr.c | 27 +- pjmedia/src/pjmedia-codec/speex/lbr_48k_tables.c | 678 +++++++++++ pjmedia/src/pjmedia-codec/speex/lpc.c | 201 ---- pjmedia/src/pjmedia-codec/speex/lpc_spx.c | 8 +- pjmedia/src/pjmedia-codec/speex/lsp.c | 6 +- pjmedia/src/pjmedia-codec/speex/lsp_bfin.h | 89 ++ pjmedia/src/pjmedia-codec/speex/ltp.c | 96 +- pjmedia/src/pjmedia-codec/speex/math_approx.c | 202 +++- pjmedia/src/pjmedia-codec/speex/math_approx.h | 14 +- pjmedia/src/pjmedia-codec/speex/mdf.c | 524 +++++--- pjmedia/src/pjmedia-codec/speex/medfilter.c | 97 ++ pjmedia/src/pjmedia-codec/speex/medfilter.h | 51 + pjmedia/src/pjmedia-codec/speex/misc.h | 4 +- pjmedia/src/pjmedia-codec/speex/modes.c | 16 +- pjmedia/src/pjmedia-codec/speex/modes.h | 4 +- pjmedia/src/pjmedia-codec/speex/nb_celp.c | 241 ++-- pjmedia/src/pjmedia-codec/speex/nb_celp.h | 10 +- pjmedia/src/pjmedia-codec/speex/preprocess_spx.c | 1258 +++++++++++--------- pjmedia/src/pjmedia-codec/speex/pseudofloat.h | 131 +- pjmedia/src/pjmedia-codec/speex/quant_lsp_bfin.h | 165 +++ pjmedia/src/pjmedia-codec/speex/sb_celp.c | 336 +++--- pjmedia/src/pjmedia-codec/speex/sb_celp.h | 2 +- pjmedia/src/pjmedia-codec/speex/speex.c | 8 +- pjmedia/src/pjmedia-codec/speex/speex.h | 17 +- pjmedia/src/pjmedia-codec/speex/speex_bits.h | 24 +- pjmedia/src/pjmedia-codec/speex/speex_callbacks.c | 12 +- pjmedia/src/pjmedia-codec/speex/speex_callbacks.h | 10 +- pjmedia/src/pjmedia-codec/speex/speex_echo.h | 57 +- pjmedia/src/pjmedia-codec/speex/speex_header.h | 7 +- pjmedia/src/pjmedia-codec/speex/speex_jitter.h | 68 +- pjmedia/src/pjmedia-codec/speex/speex_preprocess.h | 141 ++- pjmedia/src/pjmedia-codec/speex/speex_stereo.h | 6 +- pjmedia/src/pjmedia-codec/speex/vbr.c | 32 +- pjmedia/src/pjmedia-codec/speex/vorbis_psy.c | 508 ++++++++ pjmedia/src/pjmedia-codec/speex/vorbis_psy.h | 97 ++ pjmedia/src/pjmedia-codec/speex/window.c | 50 +- 51 files changed, 4104 insertions(+), 1616 deletions(-) create mode 100644 pjmedia/src/pjmedia-codec/speex/VERSION.TXT create mode 100644 pjmedia/src/pjmedia-codec/speex/filterbank.c create mode 100644 pjmedia/src/pjmedia-codec/speex/filterbank.h create mode 100644 pjmedia/src/pjmedia-codec/speex/lbr_48k_tables.c delete mode 100644 pjmedia/src/pjmedia-codec/speex/lpc.c create mode 100644 pjmedia/src/pjmedia-codec/speex/lsp_bfin.h create mode 100644 pjmedia/src/pjmedia-codec/speex/medfilter.c create mode 100644 pjmedia/src/pjmedia-codec/speex/medfilter.h create mode 100644 pjmedia/src/pjmedia-codec/speex/quant_lsp_bfin.h create mode 100644 pjmedia/src/pjmedia-codec/speex/vorbis_psy.c create mode 100644 pjmedia/src/pjmedia-codec/speex/vorbis_psy.h (limited to 'pjmedia/src/pjmedia-codec') diff --git a/pjmedia/src/pjmedia-codec/speex/VERSION.TXT b/pjmedia/src/pjmedia-codec/speex/VERSION.TXT new file mode 100644 index 00000000..743c1a21 --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/VERSION.TXT @@ -0,0 +1,9 @@ +This Speex taken from SVN + url="http://svn.xiph.org/trunk/speex", + committed-rev="12136" + name="" + committed-date="2006-11-22T02:21:19.831305Z" + last-author="jm" + kind="dir" + uuid="0101bb08-14d6-0310-b084-bc0e0c8e3800" + revision="12140" diff --git a/pjmedia/src/pjmedia-codec/speex/_kiss_fft_guts.h b/pjmedia/src/pjmedia-codec/speex/_kiss_fft_guts.h index 72acee18..43a3ba57 100644 --- a/pjmedia/src/pjmedia-codec/speex/_kiss_fft_guts.h +++ b/pjmedia/src/pjmedia-codec/speex/_kiss_fft_guts.h @@ -20,6 +20,7 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND and defines typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ #include "kiss_fft.h" +#include "math_approx.h" #define MAXFACTORS 32 /* e.g. an fft of length 128 has 4 factors @@ -140,6 +141,11 @@ struct kiss_fft_state{ (x)->r = KISS_FFT_COS(phase);\ (x)->i = KISS_FFT_SIN(phase);\ }while(0) +#define kf_cexp2(x,phase) \ + do{ \ + (x)->r = spx_cos_norm((phase));\ + (x)->i = spx_cos_norm((phase)-32768);\ +}while(0) /* a debugging function */ diff --git a/pjmedia/src/pjmedia-codec/speex/arch.h b/pjmedia/src/pjmedia-codec/speex/arch.h index 05004373..2bc5061b 100644 --- a/pjmedia/src/pjmedia-codec/speex/arch.h +++ b/pjmedia/src/pjmedia-codec/speex/arch.h @@ -39,8 +39,10 @@ #define ABS(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute integer value. */ #define ABS16(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute 16-bit value. */ +#define MIN16(a,b) ((a) < (b) ? (a) : (b)) /**< Maximum 16-bit value. */ #define MAX16(a,b) ((a) > (b) ? (a) : (b)) /**< Maximum 16-bit value. */ #define ABS32(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute 32-bit value. */ +#define MIN32(a,b) ((a) < (b) ? (a) : (b)) /**< Maximum 32-bit value. */ #define MAX32(a,b) ((a) > (b) ? (a) : (b)) /**< Maximum 32-bit value. */ #ifdef FIXED_POINT @@ -68,6 +70,7 @@ typedef spx_word32_t spx_sig_t; #define VERY_SMALL 0 #define VERY_LARGE32 ((spx_word32_t)2147483647) #define VERY_LARGE16 ((spx_word16_t)32767) +#define Q15_ONE ((spx_word16_t)32767) #ifdef FIXED_DEBUG @@ -113,6 +116,7 @@ typedef float spx_word32_t; #define VERY_SMALL 1e-15f #define VERY_LARGE32 1e15f #define VERY_LARGE16 1e15f +#define Q15_ONE ((spx_word16_t)1.f) #define QCONST16(x,bits) (x) #define QCONST32(x,bits) (x) @@ -127,6 +131,7 @@ typedef float spx_word32_t; #define SHL32(a,shift) (a) #define PSHR16(a,shift) (a) #define PSHR32(a,shift) (a) +#define VSHR32(a,shift) (a) #define SATURATE16(x,a) (x) #define SATURATE32(x,a) (x) @@ -147,6 +152,7 @@ typedef float spx_word32_t; #define MULT16_32_Q13(a,b) ((a)*(b)) #define MULT16_32_Q14(a,b) ((a)*(b)) #define MULT16_32_Q15(a,b) ((a)*(b)) +#define MULT16_32_P15(a,b) ((a)*(b)) #define MAC16_32_Q11(c,a,b) ((c)+(a)*(b)) #define MAC16_32_Q15(c,a,b) ((c)+(a)*(b)) diff --git a/pjmedia/src/pjmedia-codec/speex/cb_search.c b/pjmedia/src/pjmedia-codec/speex/cb_search.c index 0cf91057..5c688260 100644 --- a/pjmedia/src/pjmedia-codec/speex/cb_search.c +++ b/pjmedia/src/pjmedia-codec/speex/cb_search.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin +/* Copyright (C) 2002-2006 Jean-Marc Valin File: cb_search.c Redistribution and use in source and binary forms, with or without @@ -85,7 +85,7 @@ static void compute_weighted_codebook(const signed char *shape_cb, const spx_wor #endif #ifndef OVERRIDE_TARGET_UPDATE -PJ_INLINE(void) target_update(spx_word16_t *t, spx_word16_t g, spx_word16_t *r, int len) +static inline void target_update(spx_word16_t *t, spx_word16_t g, spx_word16_t *r, int len) { int n; for (n=0;n10) - N=10; - if (N<1) - N=1; params = (const split_cb_params *) par; subvect_size = params->subvect_size; @@ -152,7 +145,7 @@ int update_target ALLOC(t, nsf, spx_word16_t); ALLOC(e, nsf, spx_sig_t); - /* FIXME: make that adaptive? */ + /* FIXME: Do we still need to copy the target? */ for (i=0;i10) N=10; + /* Complexity isn't as important for the codebooks as it is for the pitch */ + N=(2*N)/3; if (N<1) N=1; - if (N==1) { - split_cb_search_shape_sign_N1(target,ak,awk1,awk2,par,p,nsf,exc,r,bits,stack,complexity,update_target); + split_cb_search_shape_sign_N1(target,ak,awk1,awk2,par,p,nsf,exc,r,bits,stack,update_target); return; } ALLOC(ot2, N, spx_word16_t*); @@ -347,7 +341,6 @@ int update_target oind[i]=itmp+(2*i+1)*nb_subvect; } - /* FIXME: make that adaptive? */ for (i=0;in; + speex_warning("FFT should not be done in-place"); for (i=0;i<((struct drft_lookup *)table)->n;i++) out[i] = scale*in[i]; } else { @@ -120,7 +120,6 @@ 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; @@ -145,8 +144,8 @@ struct kiss_config { void *spx_fft_init(int size) { struct kiss_config *table; - table = speex_alloc(sizeof(struct kiss_config)); - table->freq_data = speex_alloc(sizeof(kiss_fft_cpx)*((size>>1)+1)); + table = (struct kiss_config*)speex_alloc(sizeof(struct kiss_config)); + table->freq_data = (kiss_fft_cpx*)speex_alloc(sizeof(kiss_fft_cpx)*((size>>1)+1)); table->forward = kiss_fftr_alloc(size,0,NULL,NULL); table->backward = kiss_fftr_alloc(size,1,NULL,NULL); table->N = size; diff --git a/pjmedia/src/pjmedia-codec/speex/filterbank.c b/pjmedia/src/pjmedia-codec/speex/filterbank.c new file mode 100644 index 00000000..6da6ba49 --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/filterbank.c @@ -0,0 +1,200 @@ +/* Copyright (C) 2006 Jean-Marc Valin */ +/** + @file filterbank.c + @brief Converting between psd and filterbank + */ +/* + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 2. 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. + + 3. The name of the author may not be used to endorse or promote products + derived from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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. +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "filterbank.h" +#include "misc.h" +#include +#include "math_approx.h" + +#ifdef FIXED_POINT + +#define toBARK(n) (MULT16_16(26829,spx_atan(SHR32(MULT16_16(97,n),2))) + MULT16_16(4588,spx_atan(MULT16_32_Q15(20,MULT16_16(n,n)))) + MULT16_16(3355,n)) + +#else +#define toBARK(n) (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n)) +#endif + +#define toMEL(n) (2595.f*log10(1.f+(n)/700.f)) + +FilterBank *filterbank_new(int banks, spx_word32_t sampling, int len, int type) +{ + FilterBank *bank; + spx_word32_t df; + spx_word32_t max_mel, mel_interval; + int i; + int id1; + int id2; + df = DIV32(SHL32(sampling,15),MULT16_16(2,len)); + max_mel = toBARK(EXTRACT16(MULT16_16_Q15(QCONST16(.5f,15),sampling))); + mel_interval = PDIV32(max_mel,banks-1); + + bank = (FilterBank*)speex_alloc(sizeof(FilterBank)); + bank->nb_banks = banks; + bank->len = len; + bank->bank_left = (int*)speex_alloc(len*sizeof(int)); + bank->bank_right = (int*)speex_alloc(len*sizeof(int)); + bank->filter_left = (spx_word16_t*)speex_alloc(len*sizeof(spx_word16_t)); + bank->filter_right = (spx_word16_t*)speex_alloc(len*sizeof(spx_word16_t)); + /* Think I can safely disable normalisation that for fixed-point (and probably float as well) */ +#ifndef FIXED_POINT + bank->scaling = (float*)speex_alloc(banks*sizeof(float)); +#endif + for (i=0;i max_mel) + break; +#ifdef FIXED_POINT + id1 = DIV32(mel,mel_interval); +#else + id1 = (int)(floor(mel/mel_interval)); +#endif + if (id1>banks-2) + { + id1 = banks-2; + val = Q15_ONE; + } else { + val = DIV32_16(mel - id1*mel_interval,EXTRACT16(PSHR32(mel_interval,15))); + } + id2 = id1+1; + bank->bank_left[i] = id1; + bank->filter_left[i] = SUB16(Q15_ONE,val); + bank->bank_right[i] = id2; + bank->filter_right[i] = val; + } + + /* Think I can safely disable normalisation for fixed-point (and probably float as well) */ +#ifndef FIXED_POINT + for (i=0;inb_banks;i++) + bank->scaling[i] = 0; + for (i=0;ilen;i++) + { + int id = bank->bank_left[i]; + bank->scaling[id] += bank->filter_left[i]; + id = bank->bank_right[i]; + bank->scaling[id] += bank->filter_right[i]; + } + for (i=0;inb_banks;i++) + bank->scaling[i] = Q15_ONE/(bank->scaling[i]); +#endif + return bank; +} + +void filterbank_destroy(FilterBank *bank) +{ + speex_free(bank->bank_left); + speex_free(bank->bank_right); + speex_free(bank->filter_left); + speex_free(bank->filter_right); +#ifndef FIXED_POINT + speex_free(bank->scaling); +#endif + speex_free(bank); +} + +void filterbank_compute_bank32(FilterBank *bank, spx_word32_t *ps, spx_word32_t *mel) +{ + int i; + for (i=0;inb_banks;i++) + mel[i] = 0; + + for (i=0;ilen;i++) + { + int id; + id = bank->bank_left[i]; + mel[id] += MULT16_32_P15(bank->filter_left[i],ps[i]); + id = bank->bank_right[i]; + mel[id] += MULT16_32_P15(bank->filter_right[i],ps[i]); + } + /* Think I can safely disable normalisation that for fixed-point (and probably float as well) */ +#ifndef FIXED_POINT + /*for (i=0;inb_banks;i++) + mel[i] = MULT16_32_P15(Q15(bank->scaling[i]),mel[i]); + */ +#endif +} + +void filterbank_compute_psd16(FilterBank *bank, spx_word16_t *mel, spx_word16_t *ps) +{ + int i; + for (i=0;ilen;i++) + { + spx_word32_t tmp; + int id1, id2; + id1 = bank->bank_left[i]; + id2 = bank->bank_right[i]; + tmp = MULT16_16(mel[id1],bank->filter_left[i]); + tmp += MULT16_16(mel[id2],bank->filter_right[i]); + ps[i] = EXTRACT16(PSHR32(tmp,15)); + } +} + + +#ifndef FIXED_POINT +void filterbank_compute_bank(FilterBank *bank, float *ps, float *mel) +{ + int i; + for (i=0;inb_banks;i++) + mel[i] = 0; + + for (i=0;ilen;i++) + { + int id = bank->bank_left[i]; + mel[id] += bank->filter_left[i]*ps[i]; + id = bank->bank_right[i]; + mel[id] += bank->filter_right[i]*ps[i]; + } + for (i=0;inb_banks;i++) + mel[i] *= bank->scaling[i]; +} + +void filterbank_compute_psd(FilterBank *bank, float *mel, float *ps) +{ + int i; + for (i=0;ilen;i++) + { + int id = bank->bank_left[i]; + ps[i] = mel[id]*bank->filter_left[i]; + id = bank->bank_right[i]; + ps[i] += mel[id]*bank->filter_right[i]; + } +} +#endif diff --git a/pjmedia/src/pjmedia-codec/speex/filterbank.h b/pjmedia/src/pjmedia-codec/speex/filterbank.h new file mode 100644 index 00000000..5ded6b93 --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/filterbank.h @@ -0,0 +1,66 @@ +/* Copyright (C) 2006 Jean-Marc Valin */ +/** + @file filterbank.h + @brief Converting between psd and filterbank + */ +/* + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 2. 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. + + 3. The name of the author may not be used to endorse or promote products + derived from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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 FILTERBANK_H +#define FILTERBANK_H + +#include "misc.h" + +typedef struct { + int *bank_left; + int *bank_right; + spx_word16_t *filter_left; + spx_word16_t *filter_right; +#ifndef FIXED_POINT + float *scaling; +#endif + int nb_banks; + int len; +} FilterBank; + + +FilterBank *filterbank_new(int banks, spx_word32_t sampling, int len, int type); + +void filterbank_destroy(FilterBank *bank); + +void filterbank_compute_bank32(FilterBank *bank, spx_word32_t *ps, spx_word32_t *mel); + +void filterbank_compute_psd16(FilterBank *bank, spx_word16_t *mel, spx_word16_t *psd); + +#ifndef FIXED_POINT +void filterbank_compute_bank(FilterBank *bank, float *psd, float *mel); +void filterbank_compute_psd(FilterBank *bank, float *mel, float *psd); +#endif + + +#endif diff --git a/pjmedia/src/pjmedia-codec/speex/filters.c b/pjmedia/src/pjmedia-codec/speex/filters.c index 73cb3912..19bb1e38 100644 --- a/pjmedia/src/pjmedia-codec/speex/filters.c +++ b/pjmedia/src/pjmedia-codec/speex/filters.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin +/* Copyright (C) 2002-2006 Jean-Marc Valin File: filters.c Various analysis/synthesis filters @@ -62,6 +62,32 @@ void bw_lpc(spx_word16_t gamma, const spx_coef_t *lpc_in, spx_coef_t *lpc_out, i } } +void highpass(const spx_word16_t *x, spx_word16_t *y, int len, int filtID, spx_mem_t *mem) +{ + int i; +#ifdef FIXED_POINT + const spx_word16_t Pcoef[5][3] = {{16384, -31313, 14991}, {16384, -31569, 15249}, {16384, -31677, 15328}, {16384, -32313, 15947}, {16384, -22446, 6537}}; + const spx_word16_t Zcoef[5][3] = {{15672, -31344, 15672}, {15802, -31601, 15802}, {15847, -31694, 15847}, {16162, -32322, 16162}, {14418, -28836, 14418}}; +#else + const spx_word16_t Pcoef[5][3] = {{1.00000f, -1.91120f, 0.91498f}, {1.00000f, -1.92683f, 0.93071f}, {1.00000f, -1.93338f, 0.93553f}, {1.00000f, -1.97226f, 0.97332f}, {1.00000f, -1.37000f, 0.39900f}}; + const spx_word16_t Zcoef[5][3] = {{0.95654f, -1.91309f, 0.95654f}, {0.96446f, -1.92879f, 0.96446f}, {0.96723f, -1.93445f, 0.96723f}, {0.98645f, -1.97277f, 0.98645f}, {0.88000f, -1.76000f, 0.88000f}}; +#endif + const spx_word16_t *den, *num; + if (filtID>4) + filtID=4; + + den = Pcoef[filtID]; num = Zcoef[filtID]; + /*return;*/ + for (i=0;i16383) + { + scaledown = 1; + break; + } + } + if (scaledown) + { + for (i=0;i>shift; if (!VERIFY_INT(res)) + { fprintf (stderr, "SHR32: output is not int: %d\n", (int)res); + } spx_mips++; return res; } @@ -143,18 +145,21 @@ static inline int SHL32(long long a, int shift) long long res; if (!VERIFY_INT(a) || !VERIFY_SHORT(shift)) { - fprintf (stderr, "SHR32: inputs are not int: %d %d\n", (int)a, shift); + fprintf (stderr, "SHL32: inputs are not int: %d %d\n", (int)a, shift); } res = a<>1)),shift)) +#define PSHR32(a,shift) (SHR32((a)+((1<<((shift))>>1)),shift)) +#define VSHR32(a, shift) (((shift)>0) ? SHR32(a, shift) : SHL32(a, -(shift))) -#define PSHR16(a,shift) (SHR16(ADD16(a,(1<<((shift)-1))),shift)) -#define PSHR32(a,shift) (SHR32(ADD32(a,(1<<((shift)-1))),shift)) #define SATURATE16(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x))) #define SATURATE32(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x))) @@ -170,7 +175,9 @@ static inline short ADD16(int a, int b) } res = a+b; if (!VERIFY_SHORT(res)) - fprintf (stderr, "ADD16: output is not short: %d+%d=%d\n", a,b,res); + { + fprintf (stderr, "ADD16: output is not short: %d+%d=%d\n", a,b,res); + } spx_mips++; return res; } @@ -220,8 +227,7 @@ static inline int SUB32(long long a, long long b) #define ADD64(a,b) (MIPS_INC(a)+(b)) -#define PSHR(a,shift) (SHR((a)+(1<<((shift)-1)),shift)) - +#define PSHR(a,shift) (SHR((a)+((1<<((shift))>>1)),shift)) /* result fits in 16 bits */ static inline short MULT16_16_16(int a, int b) { @@ -264,6 +270,8 @@ static inline int MULT16_32_QX(int a, long long b, int Q) { fprintf (stderr, "MULT16_32_Q%d: inputs are not short+int: %d %d\n", Q, (int)a, (int)b); } + if (ABS32(b)>=(1<<(15+Q))) + fprintf (stderr, "MULT16_32_Q%d: second operand too large: %d %d\n", Q, (int)a, (int)b); res = (((long long)a)*(long long)b) >> Q; if (!VERIFY_INT(res)) fprintf (stderr, "MULT16_32_Q%d: output is not int: %d*%d=%d\n", Q, (int)a, (int)b,(int)res); @@ -271,6 +279,22 @@ static inline int MULT16_32_QX(int a, long long b, int Q) return res; } +static inline int MULT16_32_PX(int a, long long b, int Q) +{ + long long res; + if (!VERIFY_SHORT(a) || !VERIFY_INT(b)) + { + fprintf (stderr, "MULT16_32_P%d: inputs are not short+int: %d %d\n", Q, (int)a, (int)b); + } + if (ABS32(b)>=(1<<(15+Q))) + fprintf (stderr, "MULT16_32_Q%d: second operand too large: %d %d\n", Q, (int)a, (int)b); + res = ((((long long)a)*(long long)b) + ((1<>1))>> Q; + if (!VERIFY_INT(res)) + fprintf (stderr, "MULT16_32_P%d: output is not int: %d*%d=%d\n", Q, (int)a, (int)b,(int)res); + spx_mips+=5; + return res; +} + #define MULT16_32_Q11(a,b) MULT16_32_QX(a,b,11) #define MAC16_32_Q11(c,a,b) ADD32((c),MULT16_32_Q11((a),(b))) @@ -278,6 +302,7 @@ static inline int MULT16_32_QX(int a, long long b, int Q) #define MULT16_32_Q13(a,b) MULT16_32_QX(a,b,13) #define MULT16_32_Q14(a,b) MULT16_32_QX(a,b,14) #define MULT16_32_Q15(a,b) MULT16_32_QX(a,b,15) +#define MULT16_32_P15(a,b) MULT16_32_PX(a,b,15) #define MAC16_32_Q15(c,a,b) ADD32((c),MULT16_32_Q15((a),(b))) static inline int SATURATE(int a, int b) @@ -341,7 +366,9 @@ static inline short MULT16_16_Q15(int a, int b) res = ((long long)a)*b; res >>= 15; if (!VERIFY_SHORT(res)) + { fprintf (stderr, "MULT16_16_Q15: output is not short: %d\n", (int)res); + } spx_mips+=3; return res; } diff --git a/pjmedia/src/pjmedia-codec/speex/fixed_generic.h b/pjmedia/src/pjmedia-codec/speex/fixed_generic.h index 375050c3..2948177c 100644 --- a/pjmedia/src/pjmedia-codec/speex/fixed_generic.h +++ b/pjmedia/src/pjmedia-codec/speex/fixed_generic.h @@ -46,14 +46,15 @@ #define SHL16(a,shift) ((a) << (shift)) #define SHR32(a,shift) ((a) >> (shift)) #define SHL32(a,shift) ((a) << (shift)) -#define PSHR16(a,shift) (SHR16((a)+(1<<((shift)-1)),shift)) -#define PSHR32(a,shift) (SHR32((a)+(1<<((shift)-1)),shift)) +#define PSHR16(a,shift) (SHR16((a)+((1<<((shift))>>1)),shift)) +#define PSHR32(a,shift) (SHR32((a)+((1<<((shift))>>1)),shift)) +#define VSHR32(a, shift) (((shift)>0) ? SHR32(a, shift) : SHL32(a, -(shift))) #define SATURATE16(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x))) #define SATURATE32(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x))) #define SHR(a,shift) ((a) >> (shift)) #define SHL(a,shift) ((spx_word32_t)(a) << (shift)) -#define PSHR(a,shift) (SHR((a)+(1<<((shift)-1)),shift)) +#define PSHR(a,shift) (SHR((a)+((1<<((shift))>>1)),shift)) #define SATURATE(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x))) @@ -77,6 +78,7 @@ #define MULT16_32_Q11(a,b) ADD32(MULT16_16((a),SHR((b),11)), SHR(MULT16_16((a),((b)&0x000007ff)),11)) #define MAC16_32_Q11(c,a,b) ADD32(c,ADD32(MULT16_16((a),SHR((b),11)), SHR(MULT16_16((a),((b)&0x000007ff)),11))) +#define MULT16_32_P15(a,b) ADD32(MULT16_16((a),SHR((b),15)), PSHR(MULT16_16((a),((b)&0x00007fff)),15)) #define MULT16_32_Q15(a,b) ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15)) #define MAC16_32_Q15(c,a,b) ADD32(c,ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15))) diff --git a/pjmedia/src/pjmedia-codec/speex/jitter.c b/pjmedia/src/pjmedia-codec/speex/jitter.c index a4c07516..6d5f2ada 100644 --- a/pjmedia/src/pjmedia-codec/speex/jitter.c +++ b/pjmedia/src/pjmedia-codec/speex/jitter.c @@ -78,7 +78,7 @@ struct JitterBuffer_ { /** Initialise jitter buffer */ JitterBuffer *jitter_buffer_init(int tick) { - JitterBuffer *jitter = speex_alloc(sizeof(JitterBuffer)); + JitterBuffer *jitter = (JitterBuffer*)speex_alloc(sizeof(JitterBuffer)); if (jitter) { int i; @@ -180,7 +180,7 @@ void jitter_buffer_put(JitterBuffer *jitter, const JitterBufferPacket *packet) } /* Copy packet in buffer */ - jitter->buf[i]=speex_alloc(packet->len); + jitter->buf[i]=(char*)speex_alloc(packet->len); for (j=0;jlen;j++) jitter->buf[i][j]=packet->data[j]; jitter->timestamp[i]=packet->timestamp; diff --git a/pjmedia/src/pjmedia-codec/speex/kiss_fft.c b/pjmedia/src/pjmedia-codec/speex/kiss_fft.c index a0b3724b..c4511f8a 100644 --- a/pjmedia/src/pjmedia-codec/speex/kiss_fft.c +++ b/pjmedia/src/pjmedia-codec/speex/kiss_fft.c @@ -338,8 +338,6 @@ static void kf_factor(int n,int * facbuf) { int p=4; - double floor_sqrt; - floor_sqrt = floor( sqrt((double)n) ); /*factor out powers of 4, powers of 2, then any remaining primes */ do { @@ -349,7 +347,7 @@ void kf_factor(int n,int * facbuf) case 2: p = 3; break; default: p += 2; break; } - if (p > floor_sqrt) + if (p>32000 || (spx_int32_t)p*(spx_int32_t)p > n) p = n; /* no more factors, skip to end */ } n /= p; @@ -357,7 +355,6 @@ void kf_factor(int n,int * facbuf) *facbuf++ = n; } while (n > 1); } - /* * * User-callable function to allocate all necessary storage space for the fft. @@ -382,15 +379,22 @@ kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem int i; st->nfft=nfft; st->inverse = inverse_fft; - +#ifdef FIXED_POINT for (i=0;iinverse) - phase *= -1; - kf_cexp(st->twiddles+i, phase ); + spx_word32_t phase = i; + if (!st->inverse) + phase = -phase; + kf_cexp2(st->twiddles+i, DIV32(SHL32(phase,17),nfft)); } - +#else + for (i=0;iinverse) + phase *= -1; + kf_cexp(st->twiddles+i, phase ); + } +#endif kf_factor(nfft,st->factors); } return st; diff --git a/pjmedia/src/pjmedia-codec/speex/kiss_fftr.c b/pjmedia/src/pjmedia-codec/speex/kiss_fftr.c index b90b7254..95c45733 100644 --- a/pjmedia/src/pjmedia-codec/speex/kiss_fftr.c +++ b/pjmedia/src/pjmedia-codec/speex/kiss_fftr.c @@ -58,13 +58,22 @@ kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenme st->super_twiddles = st->tmpbuf + nfft; kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); - for (i = 0; i < nfft; ++i) { - double phase = - -3.14159265358979323846264338327 * ((double) i / nfft + .5); - if (inverse_fft) - phase *= -1; - kf_cexp (st->super_twiddles+i,phase); +#ifdef FIXED_POINT + for (i=0;i>1); + if (!inverse_fft) + phase = -phase; + kf_cexp2(st->super_twiddles+i, DIV32(SHL32(phase,16),nfft)); } +#else + for (i=0;isuper_twiddles+i, phase ); + } +#endif return st; } @@ -75,8 +84,7 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; if ( st->substate->inverse) { - speex_warning("kiss fft usage error: improper alloc\n"); - exit(1); + speex_error("kiss fft usage error: improper alloc\n"); } ncfft = st->substate->nfft; @@ -130,8 +138,7 @@ void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *t int k, ncfft; if (st->substate->inverse == 0) { - speex_warning ("kiss fft usage error: improper alloc\n"); - exit (1); + speex_error ("kiss fft usage error: improper alloc\n"); } ncfft = st->substate->nfft; diff --git a/pjmedia/src/pjmedia-codec/speex/lbr_48k_tables.c b/pjmedia/src/pjmedia-codec/speex/lbr_48k_tables.c new file mode 100644 index 00000000..b25ec15d --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/lbr_48k_tables.c @@ -0,0 +1,678 @@ +/* Copyright (C) 2002 Jean-Marc Valin + File: lbr_48k_tables.c + + 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. +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + + +const int dummy_epic_48k_variable=0; +#ifdef EPIC_48K + +const signed char gain_cdbk_ulbr[192] = { +-31, -48, -30, +-19, -10, -18, +-33, -22, -45, +-5, -56, -43, +-30, -56, -3, +-59, -17, -52, +-41, -60, -58, +-64, -47, -22, +-30, -31, -31, +-29, -14, -31, +-22, -37, -58, +-31, -44, 13, +-37, 0, 1, +-46, -55, -35, +-56, -14, -53, +-8, 1, -36, +-29, -15, -27, +-29, -39, -28, +-43, -5, 3, +-51, -27, -54, +10, -46, -36, +3, -3, -42, +-27, 16, -22, +-34, -52, 13, +-31, -21, -28, +-34, -45, -40, +-20, -48, 4, +-40, -27, 16, +-6, 11, -44, +-35, 12, -5, +19, -33, -37, +-29, 18, -32, +-29, -23, -19, +16, -47, -28, +-34, -30, 17, +-20, 2, -26, +-38, -40, -36, +15, -14, -40, +-39, 14, -9, +-15, 25, -39, +-26, 19, -32, +-39, 17, -14, +10, -36, -26, +14, -13, -40, +-29, -21, -12, +-8, 19, -39, +-36, -18, 15, +-32, -38, -38, +-19, 4, -23, +-38, -7, 11, +9, -10, -39, +-37, 24, -19, +-34, -5, -8, +-20, 23, -41, +-4, 17, -31, +-17, -26, -26, +-24, 28, -36, +-7, 15, -39, +-42, 16, -11, +-29, 14, -6, +-36, 28, -27, +-21, 5, -26, +11, -9, -39, +-38, -7, 13, +}; + + +const signed char exc_12_32_table[384] = { +34, 55, 9, 55, 4, 44, -2, 25, 4, -6, 13, -22, +20, 26, -13, -56, -37, 18, 5, 28, 4, 10, 6, -7, +37, -24, -31, 22, 12, -6, -4, -7, 2, 0, -3, -2, +-16, -13, -1, 9, -2, 4, 6, 5, -3, 3, 8, -1, +-1, -6, -2, -1, 8, 24, 19, 33, -73, -53, 6, -18, +14, 7, 11, 8, -33, -94, -5, 7, 0, 44, 1, 19, +-9, -7, -34, -16, 8, 2, 5, 0, 3, 1, -2, 3, +-22, 6, -2, 12, 16, 30, 39, 25, 25, 2, 10, -2, +-1, -40, -6, -51, -5, -48, -9, -33, -14, -1, -24, 15, +104, 39, 12, -9, -20, -12, -30, -10, -31, -7, -30, -8, +-71, -53, -4, -11, 9, -10, 7, -10, 10, -1, 11, 8, +24, 14, 6, -3, 10, 8, 8, 11, -6, 11, 0, -2, +-6, -2, 1, -1, -3, 8, -41, 27, 57, -7, 11, -16, +-61, 50, 10, -10, 4, -13, 14, -7, 1, 5, -4, 4, +0, 2, -1, -2, -1, 1, 1, 0, -1, -1, -2, -3, +-3, -15, 69, 60, 10, -10, -10, -29, -21, -7, -16, 2, +24, -32, 24, -18, -14, -2, -11, 11, -6, 10, 1, 3, +24, -10, 14, 18, -13, 17, -16, 4, -3, -21, -3, -11, +-19, 12, -14, 26, 20, -9, 24, -15, 18, 1, -32, -2, +-1, 8, -3, 4, 11, -47, 7, 46, -4, -10, -10, -2, +-24, 29, -33, 6, -20, -3, 0, -12, 5, -30, 8, -13, +28, 9, 5, -11, 0, -14, -13, -22, -12, -8, -4, 1, +-6, 28, 45, -18, -31, -5, 1, 2, 1, 5, 0, -3, +-19, -10, 10, 27, 8, -16, -28, -9, 2, -5, 8, -1, +100, -49, 4, -43, 25, -7, 1, 9, -13, 13, -18, 13, +-1, -1, 0, 2, -2, -8, 9, -46, -7, 70, 23, 7, +-103, 20, 8, 42, -5, 21, -4, 4, 1, -8, 16, -8, +3, 3, 8, 4, 7, -3, -3, -4, 9, 6, 2, 13, +6, 3, -15, 11, -43, 31, 40, -13, 12, -21, -2, -3, +-10, -9, 16, -35, 31, -3, -12, 8, -34, 7, 12, 22, +-3, -4, -7, -12, 24, 53, -19, -43, 4, -3, -4, 6, +-18, -30, -58, -17, -11, 17, 23, 34, 30, 28, 28, 15, +}; + + +const signed char cdbk_lsp_vlbr[5120]={ +23, 34, 108, 100, 102, 82, 69, 48, 52, 25, +0, -37, -55, -78, -111, -79, 58, 57, 45, 32, +27, -9, -12, -14, -41, -29, -17, -41, 44, 35, +-24, -68, -72, 61, 100, 73, 100, 80, 70, 37, +12, -5, 22, 11, -10, -40, -33, -17, 19, 12, +-20, -57, -94, -92, 56, 71, 48, 31, 22, -5, +41, 28, 6, -6, -12, -39, -18, -16, -30, -23, +65, 54, 41, 28, 23, 9, 26, 18, 22, 6, +17, -16, -33, -54, -87, -79, 8, -8, 44, 35, +-20, -62, -78, 22, 78, 47, 44, 33, 26, 14, +8, 1, 45, 47, 72, 68, 55, 31, 36, 17, +-27, -68, -86, -65, -10, 23, 8, -22, -31, 25, +-4, -38, -55, -68, -96, -118, -39, 30, 28, 31, +-21, -66, -47, 99, 91, 68, 78, 56, 64, 36, +33, 22, 13, -13, -36, -22, 44, 37, 54, 33, +-31, -76, -106, -100, -5, 21, 7, -17, 13, 48, +-26, -65, -84, -84, -46, 67, 97, 66, 58, 31, +-20, -52, -32, -20, 3, 16, 27, 40, 54, 29, +-6, -35, -56, -64, -8, -31, -36, 21, 26, -3, +32, 23, 1, -23, -19, -44, -45, -7, 10, -10, +-24, -55, 2, 67, 72, 85, 90, 74, 77, 45, +-21, -58, -45, -49, 16, 34, 13, -15, -16, 16, +8, -31, -34, -61, -83, 10, 24, 8, 56, 25, +-8, -49, -74, -95, -123, -77, 6, 40, 46, 42, +-21, -60, -59, -34, -12, 27, 8, -19, -48, -17, +-25, -66, -78, -73, -81, -16, 14, 0, -2, 33, +78, 79, 69, 49, 44, 32, 50, 44, 46, 22, +24, 9, -4, -18, -37, -56, 22, 34, 22, 11, +-19, -59, -85, -41, 46, 72, 60, 33, 29, -3, +-21, -66, -70, 65, 92, 57, 61, 41, 40, 23, +-4, -41, -60, -72, -102, -106, 4, 56, 57, 31, +-5, -48, -62, -91, -109, 1, 76, 54, 72, 39, +-21, -61, -86, -46, -34, -39, 42, 25, 15, 12, +5, -16, -36, -56, 5, 18, 11, 13, 52, 23, +12, -6, 30, 40, 59, 40, 27, 8, 19, 6, +25, 8, -9, -19, -25, -53, -40, -38, -46, -4, +-17, -59, -83, 2, 58, 29, 18, -2, -17, -5, +-35, -80, -111, -117, -41, -9, 14, 23, 36, 56, +48, 67, 93, 71, 77, 91, 110, 95, 83, 47, +-25, -62, -97, -93, 76, 96, 73, 52, 61, 28, +-9, -55, -46, 49, 33, 8, 1, -25, 28, 23, +-10, -47, -60, -45, -62, -58, 56, 57, 48, 28, +34, 12, -9, 0, 34, 4, 6, 11, 3, -18, +21, -9, -16, -13, -39, -41, 14, -8, 33, 28, +-7, -49, -61, 15, 34, 3, 2, -13, -28, -17, +-14, -50, -46, -65, -76, -13, -10, -29, -30, 22, +-28, -68, -97, -98, -8, 38, 36, 26, 25, 15, +6, -33, -9, 1, -28, -11, -19, -24, 61, 36, +-15, -60, -19, 81, 58, 52, 42, 28, 66, 36, +-15, -52, -71, -15, 11, -13, 38, 28, 11, -4, +34, 1, -9, -27, -57, -19, 36, 6, 14, 0, +-1, -43, -14, 16, -12, -5, -14, -29, -33, -32, +-13, -57, -75, -100, -111, 1, 2, 13, 48, 33, +12, 16, 100, 85, 69, 49, 40, 29, 46, 23, +-4, -26, -41, -44, -7, -26, -39, -27, 18, 0, +-4, -47, -51, 17, 7, -19, 13, -10, -16, 9, +-24, -63, -93, -53, 25, 14, 73, 51, 35, 8, +-34, -77, -106, -83, -51, -47, 2, 12, 41, 53, +-13, -47, -67, -44, 42, 20, 24, 33, 21, -3, +11, -15, -29, -51, -79, -88, 22, 56, 43, 20, +11, -22, -37, -1, 61, 40, 28, 24, 22, -6, +-3, -33, -50, -66, -93, -100, -16, -16, 3, 41, +-18, -58, -82, -5, 95, 78, 56, 39, 30, 1, +-6, -47, -28, -26, -36, 49, 55, 51, 71, 35, +-6, -50, -42, -4, -32, -1, -1, -18, 67, 40, +-23, -63, -56, -48, -32, 0, -14, -43, -46, 25, +-17, -61, -63, 13, -1, 28, 23, 10, 67, 36, +45, 92, 124, 111, 108, 86, 77, 56, 57, 28, +50, 35, 13, 3, -2, -32, 3, 14, 6, -8, +12, -17, -24, -42, -67, -23, 67, 49, 64, 38, +-21, -60, -90, -45, 32, 6, 7, -3, -15, 9, +-16, -62, -73, 50, 46, 18, 7, -13, 63, 39, +19, -16, -19, 20, 5, -15, 16, -9, 5, 8, +-11, -46, -42, -39, -55, -68, -62, -27, -18, 23, +-23, -61, -67, -71, -29, 44, 32, 10, -15, -12, +-6, -45, -43, -40, -67, -22, 42, 19, 61, 38, +9, -13, -38, -37, 40, 30, 15, 9, 11, -16, +0, -18, -29, -34, -17, -44, -50, -3, 47, 15, +-3, -46, -26, 20, -10, 16, 20, -2, 43, 18, +-23, -46, 46, 91, 99, 100, 99, 79, 72, 42, +-1, -44, -33, -36, -56, 22, 17, 4, 71, 37, +0, -38, -49, 0, -1, -30, -21, -35, -44, -6, +-32, -74, -101, -98, -14, -21, -23, 7, 26, 45, +8, -28, -44, -63, -96, -84, 34, 21, 13, 23, +10, -24, -38, -17, -29, -53, -16, -41, -14, 23, +-19, -61, -76, -12, 97, 99, 79, 60, 59, 25, +5, -11, -26, -54, -26, -8, -13, 3, 25, 4, +17, 16, 4, -24, -1, 42, 60, 63, 70, 37, +10, -27, -22, -43, -62, -7, -16, 10, 75, 40, +-26, -64, -96, -106, -3, 73, 73, 46, 55, 29, +-2, -45, -43, 45, 46, 17, 37, 10, 30, 32, +21, -2, -18, -28, -47, -63, -38, -56, -3, 27, +-17, -48, -9, 9, 3, 28, 50, 58, 73, 44, +8, -14, -32, -56, -81, -106, -35, 41, 53, 26, +1, -38, -51, -66, -100, -61, 32, 17, 66, 42, +17, 5, -6, -21, -26, -52, -36, 23, 56, 22, +6, -20, -31, -22, -19, -48, 16, 38, 22, -2, +-25, -67, -93, -51, 82, 62, 71, 69, 63, 35, +-12, -51, -71, -60, -76, -91, -14, 41, 35, 20, +-16, -31, 22, 32, 55, 80, 98, 91, 85, 49, +-21, -63, -92, -6, 57, 27, 36, 11, 60, 39, +-7, -45, -67, -81, -114, -110, 0, 24, 23, 45, +-18, -55, -61, -56, -60, -64, -32, -2, 58, 36, +-3, -33, -51, -24, 19, -10, -19, 4, 14, -15, +-19, -59, -81, -12, 7, -12, 36, 16, 48, 36, +17, -13, -32, -49, -78, -95, -1, -3, -10, 25, +15, 5, 41, 59, 108, 101, 103, 81, 70, 35, +-14, -52, -37, 15, 93, 83, 66, 50, 47, 15, +-3, -31, -49, -52, -9, -31, -10, 37, 62, 27, +-15, -56, -82, -17, 75, 56, 36, 22, 7, -16, +-24, -63, -93, -84, 25, 94, 98, 65, 60, 31, +-2, -45, -39, -61, -61, 48, 35, 32, 54, 18, +-19, -51, -45, -57, -28, -8, 10, 14, 38, 26, +-2, -46, -38, 45, 26, 22, 48, 21, 63, 40, +-22, -61, -73, -75, -67, -31, 13, 18, 51, 34, +-12, -2, -1, -17, -3, -27, -3, 6, -1, -15, +-16, -59, -78, 10, 36, 9, 4, -18, 33, 22, +-25, -62, -97, -107, 39, 87, 69, 46, 42, 12, +11, -7, -30, -36, 19, 2, -10, -7, -4, -24, +11, -8, 25, 28, 28, 5, -4, -10, 5, -2, +-10, -48, -37, -17, -38, -9, -2, -19, -30, -22, +-23, -61, -79, -81, -2, 15, -4, 17, 20, 2, +25, -14, -3, -10, -38, 1, 14, -14, -9, -27, +2, -18, -38, -36, -11, -39, -36, -28, -36, -11, +32, 59, 127, 124, 127, 108, 91, 68, 64, 34, +3, -32, -37, -10, -25, -46, 12, 1, -17, -24, +-29, -69, -102, -100, 2, -7, 11, 14, 1, 31, +30, -6, -4, -16, -44, -5, 3, -9, 66, 40, +-9, -45, -52, -5, 37, 19, 26, 6, 51, 32, +-31, -73, -96, -45, -25, -37, -15, -16, 32, 39, +3, -15, 18, 21, 28, 33, 58, 58, 69, 42, +-31, -73, -99, -99, -48, 14, 21, 5, 2, 39, +7, -35, -20, 29, 2, -8, -8, -28, 38, 26, +-5, -39, -64, -36, 15, -15, -11, -21, -23, 5, +-8, -51, -56, 15, -1, -14, -8, -31, 36, 22, +-8, -53, -68, -98, -101, 42, 49, 38, 41, 12, +10, -27, -22, 4, -23, -21, 30, -1, 22, 26, +-13, -56, -42, 31, 9, -1, -10, -2, 22, -4, +15, 8, 56, 57, 45, 55, 57, 46, 72, 44, +-7, -53, -26, 53, 21, 17, 0, 0, 74, 41, +3, -18, -2, 0, 19, 17, 42, 36, 47, 26, +24, -7, -23, -34, -62, -60, 6, -22, 18, 25, +-11, -42, -46, -61, -83, -99, -67, -11, 28, 39, +30, -3, -10, -1, -24, -30, -1, -28, 15, 18, +19, -15, -10, -6, -35, -26, 33, 10, 56, 39, +-13, -53, -82, -42, 53, 37, 18, 10, -3, -21, +-21, -60, -89, -46, 89, 94, 71, 46, 42, 9, +-2, -34, -44, -46, -64, -84, -1, 37, 16, 0, +-17, -51, -65, -64, -7, -17, -29, -11, 52, 27, +22, -15, -16, -39, -55, 26, 36, 21, 62, 28, +-2, -26, -38, -49, -55, -80, -75, 8, 20, 9, +-6, -47, -61, -82, -103, -17, -15, -25, 53, 40, +-8, -47, -66, -18, 56, 43, 25, 29, 39, 3, +-27, -66, -86, -69, -50, -59, -34, -1, 19, 42, +3, -20, 2, 21, 72, 57, 52, 36, 31, 7, +-12, -49, -61, -13, -1, -33, 5, 37, 26, 2, +-27, -69, -92, -62, 2, 43, 88, 67, 64, 36, +0, -40, -4, -6, -20, 43, 33, 25, 50, 20, +14, -20, -30, -44, -73, -37, -24, -47, 26, 20, +31, 53, 111, 118, 127, 126, 121, 99, 85, 46, +-14, -45, -51, -39, 24, 5, -6, 17, 46, 14, +-4, -43, -45, -70, -63, 8, 14, 58, 78, 39, +-8, -47, -66, -84, -114, -55, 10, -8, 32, 40, +28, 22, 42, 26, 8, -21, -16, -6, 22, 10, +24, 10, 34, 31, 35, 31, 46, 39, 59, 36, +-4, -43, -62, -10, 20, -14, 2, 14, -6, -19, +-21, -62, -89, -22, 62, 33, 30, 16, 15, 15, +0, -22, -31, -45, -58, -80, -66, 13, 68, 34, +-16, -45, -6, 7, -6, -17, -14, -15, 2, 2, +10, -5, -22, -38, -40, -70, -60, -15, -23, 0, +22, -11, -22, -39, -67, -25, 30, 5, 58, 37, +-21, -47, 6, 43, 37, 45, 65, 66, 73, 43, +2, -25, -40, -53, -72, -94, -35, 24, 9, 8, +-3, 0, -3, -9, 4, -23, -10, 20, 43, 14, +-2, -41, -60, -9, 57, 32, 17, 16, 6, -19, +1, -31, -36, -36, -54, -68, -77, -75, 21, 37, +-19, -32, 79, 90, 92, 81, 67, 47, 52, 28, +-6, -36, -57, -62, 27, 40, 21, 11, 9, -19, +-10, -47, -49, -59, -74, -18, -14, -30, 25, 18, +-23, -69, -82, 60, 66, 40, 75, 54, 65, 38, +-19, -57, -92, -68, 66, 58, 34, 18, 1, -16, +-29, -68, -99, -88, -37, -38, 13, 8, 5, 40, +-22, -63, -75, 14, 15, 7, 75, 58, 59, 34, +-23, -62, -82, -39, -31, -53, -27, 5, -3, 20, +13, -26, -20, 22, 2, -3, 35, 13, 54, 39, +32, 5, -13, -22, -45, -58, -1, -20, -19, 7, +30, 46, 70, 55, 89, 88, 91, 67, 56, 28, +-13, -50, -63, -25, -28, -50, -23, -32, -34, 19, +-13, -54, -65, -9, -20, -37, 29, 6, 11, 25, +0, -40, -55, -78, -107, -25, 47, 20, 34, 16, +-20, -58, -96, -103, 38, 43, 27, 30, 15, -1, +-16, -49, -52, -66, -80, -57, -44, -39, 6, 38, +0, -38, 5, 13, -8, 23, 24, 1, 7, -9, +-18, -56, -64, -7, 38, 13, 11, 32, 28, 0, +14, -3, -20, -17, 4, -26, -34, -8, 19, -9, +-23, -60, -83, -38, -8, -32, 11, 19, -1, -5, +-5, -47, -12, 56, 38, 22, 18, -8, -5, -8, +18, -4, -24, -16, 27, 2, -6, 5, 25, -5, +13, 0, -19, -35, -23, -45, -59, -30, 19, 3, +19, -12, -23, 1, -7, -35, -14, -32, -23, 4, +-23, -64, -67, -22, -27, -5, -5, -20, 20, 5, +20, 11, 83, 92, 85, 89, 69, 53, 80, 48, +15, -2, -21, -29, -18, -48, -52, -12, -11, -21, +-6, -38, -55, -68, -9, 33, 22, 19, 25, -1, +-8, -46, -49, -67, -64, 16, 8, -6, 32, 15, +3, -25, -46, -46, 39, 50, 34, 21, 46, 14, +8, -33, -37, -68, -82, 31, 34, 13, 19, -6, +33, 0, 5, -7, -32, 2, 22, -3, 35, 17, +-23, -62, -91, -64, 6, 3, 36, 26, 7, -3, +-12, -54, -60, 26, 46, 16, 30, 22, 8, -4, +-23, -61, -40, 31, 58, 73, 88, 77, 74, 41, +-2, -42, -49, 13, 5, -15, 22, -4, 26, 27, +-13, -54, -39, 18, 2, -8, -12, 34, 56, 23, +-20, -31, 27, 23, 24, 28, 39, 33, 47, 27, +36, 17, -4, -20, -30, -61, -8, 20, 0, -15, +-10, -51, -72, -82, -111, -73, 34, 25, 19, 38, +-10, -45, -63, -55, -46, -75, -45, 34, 34, 12, +6, -18, 29, 26, 7, -9, 0, 5, 38, 22, +-7, -52, -16, 69, 43, 26, 23, 2, 51, 34, +-12, -51, -59, -78, -88, 15, 20, 0, -14, 12, +-3, -36, -59, -45, 60, 49, 28, 20, 16, -13, +-28, -70, -90, 9, 67, 48, 90, 77, 70, 38, +-10, -39, -58, -54, 15, -12, 3, 35, 27, -3, +12, -1, 28, 29, 55, 53, 80, 65, 51, 23, +-17, -61, -39, 74, 56, 43, 75, 51, 58, 36, +-30, -71, -93, -43, -29, -26, 4, -19, -14, 37, +3, -13, -31, -38, 11, -5, -22, -11, 43, 14, +-25, -65, -80, -79, -71, 3, 37, 32, 20, 9, +-20, -60, -77, -26, 18, 43, 44, 24, 22, -3, +-4, -42, -22, -19, -45, -32, -35, -39, -46, 1, +-25, -59, -27, -10, -7, -4, 7, 13, 25, 12, +8, -25, -32, -47, -74, -32, 27, 6, 25, 7, +41, 40, 62, 64, 64, 50, 54, 42, 49, 25, +-21, -63, -88, -21, 16, -3, -4, -26, 57, 38, +8, -25, -34, 2, -8, -28, 2, -22, 12, 23, +-19, -49, 10, 71, 84, 71, 66, 48, 42, 22, +-20, -58, -89, -57, 62, 44, 33, 36, 25, -1, +-22, -55, -27, 1, 43, 37, 46, 50, 51, 26, +1, -38, -46, 22, 34, 4, 20, -2, 3, 9, +-4, -42, -49, -75, -89, -24, -25, 19, 71, 39, +5, -28, -45, -43, -63, -75, -17, -38, 14, 30, +-4, -36, -62, -59, -29, -43, -4, -16, 11, 23, +-19, -57, -82, -39, 26, 2, -2, 20, 11, -10, +-28, -68, -92, -70, 9, -1, -15, -30, 11, 31, +1, -22, -41, -49, -30, -58, -48, 8, 4, -9, +38, 41, 108, 115, 96, 98, 103, 84, 86, 51, +15, 1, 58, 46, 26, 6, 16, 18, 41, 24, +4, -34, -14, -27, -42, 20, 18, 2, 23, 1, +-22, -59, -83, -70, -22, -42, -26, 29, 29, 15, +-14, -34, 11, -1, -21, -35, -3, 1, 29, 16, +-16, -57, -78, -7, 17, -13, 8, -13, -6, 22, +-22, -32, -21, -20, 20, -4, 10, 13, 12, -4, +8, -30, -30, -46, -71, -4, 3, -11, 4, -11, +16, 5, -15, -21, 3, -23, -25, -19, -28, -32, +-28, -68, -98, -101, -34, 19, 71, 52, 49, 30, +-18, -57, -82, -56, -56, -66, 15, 12, 1, 29, +-21, -62, -76, -27, -33, -38, 18, 30, 54, 32, +3, -36, -10, -17, -34, -3, -8, 32, 63, 27, +1, -30, -44, -20, -13, -49, -25, 3, -14, -18, +-26, -68, -80, -46, -28, 17, 42, 37, 58, 34, +30, 26, 57, 55, 49, 25, 16, 3, 24, 11, +35, 35, 67, 57, 60, 82, 114, 103, 93, 55, +18, -8, -23, -32, -53, -68, 15, 11, -6, -7, +-2, -43, -29, 0, -28, -5, -5, -15, 25, -1, +4, -13, -35, -45, 14, -6, -2, 19, 16, -9, +-30, -72, -93, -93, -73, -10, -15, 6, 30, 45, +-23, -58, -50, -55, -74, -60, -23, 0, 6, 21, +-4, -40, -63, -24, 7, -19, 4, -18, 27, 28, +-12, 1, 88, 76, 74, 88, 93, 90, 80, 44, +-13, -59, -43, 52, 27, 21, 15, 12, 42, 11, +22, -14, -17, -33, -57, -4, 5, -18, 40, 18, +-3, -23, -43, -44, 8, -16, -14, -4, -20, -29, +35, 45, 75, 82, 111, 117, 125, 105, 89, 51, +-3, -38, -57, -30, 79, 71, 48, 33, 33, 0, +-17, -62, -57, 66, 67, 35, 29, 5, 22, 17, +3, -31, -34, -21, -44, -49, 42, 23, 24, 26, +-23, -62, -74, -49, -30, -30, -37, -50, 9, 35, +-17, -57, -71, -81, -45, 61, 58, 37, 31, 9, +3, -7, 28, 14, -2, 0, 40, 41, 58, 33, +-11, -51, -74, -17, 40, 12, 8, 13, -4, -22, +-16, -46, -31, -35, -49, -49, -26, -9, -7, -7, +17, -9, -24, -41, -68, -73, 38, 33, 19, 16, +-15, -50, -47, -16, -24, -21, 59, 56, 53, 30, +-14, -54, -57, 2, -17, -33, -34, -21, 4, -4, +-23, -62, -93, -72, 48, 31, 21, 6, 3, 17, +-18, -63, -79, 44, 68, 36, 45, 20, 57, 37, +-29, -72, -99, -111, -86, -31, 7, 25, 39, 55, +-14, -49, -53, -63, -80, -31, 24, 13, 1, -1, +-9, -45, -55, -27, -31, -63, -23, 25, 13, -5, +-20, -61, -80, 7, 44, 16, 54, 40, 32, 17, +24, 7, -8, -43, -62, -54, -11, 7, 35, 27, +-12, -55, -59, -48, -69, -4, -1, -12, 68, 39, +-12, -31, 52, 63, 53, 34, 29, 22, 36, 19, +-26, -66, -97, -79, 50, 41, 40, 48, 54, 28, +-2, -37, -41, -2, -11, -30, 29, 16, 4, -2, +40, 49, 56, 37, 39, 40, 64, 59, 67, 39, +11, -5, 20, 14, 25, 16, 25, 22, 37, 17, +-3, -43, -46, -10, -35, -38, -35, -39, 67, 43, +-7, -47, -33, -39, -60, -12, -18, 11, 43, 11, +-25, -65, -91, -76, -91, -81, 0, 13, 34, 50, +-9, -50, -52, 17, 0, -4, 43, 18, 63, 42, +-8, -15, 15, 41, 56, 35, 51, 45, 51, 29, +0, -14, -24, -36, -43, -70, -39, 27, 33, 5, +-25, -62, -81, -66, -12, -26, -16, -4, -13, 21, +-29, -68, -60, -24, -3, 11, 18, 19, 30, 20, +1, -35, -42, -30, -57, -51, 13, -17, 3, 22, +-8, -27, -12, -2, -7, -21, 36, 41, 34, 12, +-17, -56, -62, -72, -73, -17, -26, 9, 16, 13, +11, -21, -37, -3, 16, -17, 1, -3, -18, -19, +15, -20, -19, -22, -49, -30, -7, -29, 3, -2, +17, -4, 11, 6, 51, 40, 36, 34, 48, 22, +-19, -55, -29, 37, 68, 49, 45, 33, 42, 23, +7, -30, -22, 3, -22, -36, -36, -54, 20, 22, +20, 2, -15, -39, -59, -85, -10, 37, 21, 2, +-15, -54, -77, -54, 74, 70, 48, 32, 51, 20, +-25, -64, -70, -75, -52, 17, 6, -20, -30, 26, +-13, -55, -15, 39, 16, 42, 30, 33, 62, 28, +-21, -56, -30, -35, 6, 13, -4, -29, 33, 27, +-17, -55, -75, -31, 3, -28, -26, 16, 18, -4, +-13, -44, -60, -52, -9, -36, -38, 1, -9, -9, +-12, -50, -77, -70, 43, 47, 28, 13, 43, 16, +-13, -57, -80, -104, -113, -30, 43, 45, 52, 39, +3, -28, -42, -37, -58, -67, 23, 4, 38, 33, +-21, -64, -74, -22, 43, 83, 81, 56, 62, 32, +34, 26, 23, 9, 14, 17, 26, 16, 37, 19, +-5, -48, -49, -75, -65, 9, -6, 41, 45, 15, +32, 30, 63, 83, 90, 91, 100, 84, 85, 52, +-19, -54, -68, -71, 11, 30, 13, 1, 63, 37, +44, 76, 99, 87, 117, 113, 103, 77, 64, 32, +-16, -45, -12, 30, 27, 15, 57, 49, 42, 22, +9, -15, -31, -28, -36, -61, -13, -18, -33, -5, +-12, -55, -8, 82, 64, 47, 42, 21, 27, 9, +-15, -56, -74, -12, -19, -28, 7, -16, 53, 38, +-7, -45, -64, -55, -74, -80, 35, 45, 24, 15, +-25, -48, -34, -42, -6, -27, -9, 9, 13, -7, +-25, -65, -84, -35, 30, 14, 24, 39, 48, 28, +-22, -62, -86, -51, 64, 104, 94, 61, 62, 31, +16, -15, -25, -28, -55, -56, -10, -23, 52, 36, +10, -10, -22, -46, -71, -92, -45, -13, 16, 26, +-22, -65, -84, 17, 85, 55, 66, 48, 55, 31, +-1, -38, -33, -25, -49, -15, 18, -10, 41, 30, +-3, -24, -47, -60, -30, -46, -17, -13, -27, 1, +-7, -41, -61, -54, -50, -78, -28, 5, -2, 20, +17, -12, -27, -4, 8, -27, -5, 23, 14, -12, +-36, -81, -111, -75, -17, -9, 9, 9, 39, 49, +-13, -59, -54, 68, 51, 32, 35, 14, 64, 38, +-11, -44, -69, -57, 11, -11, -5, 3, -16, -15, +34, 36, 41, 16, -8, -24, 11, 23, 48, 28, +-17, -42, 8, 17, 45, 69, 71, 55, 49, 23, +3, -30, -46, -64, -95, -109, 2, 39, 19, 19, +25, -13, -7, 16, -11, -5, 8, -11, 52, 33, +-8, -37, -57, -60, 13, 7, -14, -4, 20, -7, +7, -31, 9, 44, 20, 22, 29, 10, 52, 31, +3, -22, -36, -53, -80, -77, -35, -41, 54, 41, +-21, -59, -87, -83, 12, 69, 57, 36, 32, 2, +6, -14, -34, -42, -4, -32, -27, 10, 4, -20, +-11, -56, -59, 25, 8, -5, -9, -26, 68, 43, +22, 13, 40, 39, 73, 81, 95, 88, 82, 45, +-18, -62, -79, 28, 60, 28, 29, 3, 23, 25, +6, -31, -39, -55, -85, -32, 7, -17, 48, 30, +7, -24, -42, -8, 39, 10, 7, 11, 1, -20, +-1, -36, -26, -30, -48, 2, 46, 26, 35, 14, +-17, -57, -56, -10, -4, 26, 22, 6, -4, -16, +-18, -55, -59, -67, -86, -50, 3, 29, 29, 16, +-25, -61, -33, 2, 26, 25, 23, 10, 24, 16, +26, 15, -8, -20, 6, -21, -16, 3, 0, -22, +13, -17, -26, -12, -31, -48, 15, -3, 1, 14, +4, -27, -33, -21, -29, -53, -52, -64, -42, 22, +-11, -41, -42, -40, -42, -64, 0, 48, 50, 19, +-13, -47, -42, -56, -50, 10, 3, -9, -30, -16, +-4, -47, -12, 14, -13, 21, 13, 6, 73, 40, +-15, -50, -63, -41, -31, -55, -60, -13, 28, 15, +-6, -50, -56, -81, -95, 30, 29, 21, 71, 35, +-14, -58, -65, 37, 40, 9, 18, -10, -4, 20, +31, 33, 79, 106, 119, 103, 100, 77, 64, 33, +14, 7, 56, 81, 97, 85, 85, 61, 45, 20, +-24, -66, -74, -51, -17, 16, 5, -21, 22, 26, +-1, -25, -38, -24, -2, -35, -26, 21, 34, 1, +20, 15, 75, 59, 39, 26, 48, 43, 50, 29, +26, 1, -15, -4, -9, -38, 9, 2, -9, -8, +20, 14, 19, 13, 2, -16, 24, 25, 19, 3, +-8, -41, -58, -78, -109, -106, -27, 9, 53, 46, +17, 5, -13, -25, -24, -54, -20, 2, -21, -25, +-11, -50, -48, 19, 10, -7, 46, 26, 24, 17, +-28, -68, -69, -50, -49, -33, -7, 10, 20, 21, +4, -36, -21, 38, 19, 1, 25, -3, 20, 22, +-18, -58, -87, -44, 82, 73, 49, 31, 19, -7, +-21, -61, -78, -44, -58, -66, -9, -23, 10, 40, +-25, -66, -76, -33, -33, 6, 16, -6, -21, 5, +-27, -69, -77, 5, -2, -7, 6, 9, 24, 6, +1, -41, -41, -58, -79, 25, 69, 43, 68, 36, +-28, -67, -85, -71, -34, 0, 14, 12, -2, 4, +-2, -33, -55, -37, 37, 25, 6, 6, 29, -4, +-8, -25, -11, -9, 40, 23, 17, 25, 37, 12, +-21, -27, 52, 60, 47, 58, 76, 70, 69, 38, +23, 15, 22, -8, -32, -50, -10, 3, 31, 21, +-10, -44, -67, -61, -29, -54, 8, 39, 21, 4, +31, 18, 30, 36, 46, 28, 50, 42, 35, 13, +-21, -57, -24, -16, -15, 14, 3, -25, -17, 20, +23, 6, 24, 35, 90, 72, 64, 55, 55, 23, +7, -28, -42, -19, -36, -50, -15, -40, 29, 28, +-21, -59, -66, -59, -6, 85, 83, 53, 54, 26, +2, -20, -42, -52, -27, -49, 5, 9, -9, -8, +-18, -56, -84, -72, 24, 51, 32, 20, 13, -10, +-13, -53, -64, -39, -63, -40, 24, 0, 34, 29, +0, -31, -45, -63, -90, -53, 3, -18, -9, 24, +-13, -57, -71, 27, 19, -3, 25, -3, 45, 35, +12, -9, -26, -40, -51, -78, -24, 11, -9, -6, +-12, -49, -45, -33, -50, -48, -46, -52, 2, 25, +-14, -56, -84, -108, -122, -50, 4, 22, 42, 53, +-6, -44, -54, -28, -41, -61, 19, 25, 6, 1, +-32, -75, -95, -38, -1, -15, 7, 14, 23, 23, +11, 12, 37, 30, 38, 51, 80, 80, 82, 47, +-19, -56, -69, -82, -98, -64, -29, 2, 28, 42, +-18, -49, 3, 34, 41, 39, 32, 18, 21, 7, +-8, 23, 39, 30, 30, 27, 41, 36, 44, 23, +-16, -49, -69, -46, 1, -27, 41, 48, 35, 15, +6, -32, -36, 5, -17, -30, 2, -16, 51, 35, +-23, -64, -91, -21, 71, 44, 52, 44, 40, 21, +-22, -55, -40, -20, 62, 52, 38, 29, 27, 5, +-27, -69, -75, -6, -8, 3, -2, -30, -42, 13, +2, -29, -42, -10, 29, 2, 15, 30, 26, 0, +-27, -69, -85, -75, -54, -8, -14, -31, 16, 42, +-2, -44, -54, -75, -101, -6, 4, -3, 32, 13, +7, -3, -22, -43, -14, -27, -24, -11, -14, -12, +-18, -57, -85, -66, 47, 86, 75, 45, 42, 10, +18, 3, 39, 76, 80, 48, 48, 41, 44, 21, +-13, -51, -73, -22, -15, -41, 3, -4, -22, 3, +-10, -50, -63, -9, -23, -42, -6, -30, 11, 28, +15, -11, -29, -39, -57, -71, -23, -40, -33, 20, +-2, -40, -48, -5, -21, -25, 31, 7, 53, 35, +-19, -63, -72, 39, 34, 16, 37, 17, 54, 31, +-18, -49, -28, -27, -40, -29, 4, 4, 19, 11, +-1, -43, -41, 24, 15, -12, -4, -31, -4, 23, +-11, -46, -79, -74, 34, 21, 9, 15, 3, -13, +-16, -51, -56, -55, 8, 62, 50, 30, 43, 14, +-1, -25, -27, -32, -46, -62, -66, -36, 53, 32, +-12, -49, -77, -50, 2, -21, 23, 6, 14, 21, +-5, -47, -58, -77, -105, -14, 30, 9, 73, 44, +-24, -52, -4, 10, 13, 13, 27, 27, 42, 24, +-10, -38, -27, -15, -24, -52, -53, 1, 14, -6, +-17, -45, -13, 2, 19, 53, 83, 79, 76, 43, +-21, -62, -86, -48, 40, 17, 14, -9, 40, 30, +-1, -32, -51, -33, -3, -35, 2, 17, -2, -12, +-21, -60, -85, -70, 33, 73, 58, 37, 67, 36, +-1, -45, -39, 37, 17, 3, 18, -8, 53, 35, +-8, -47, -65, -61, -87, -93, 9, 9, 2, 33, +-13, -55, -56, -19, -29, 14, 17, 6, 55, 28, +5, 1, -7, -23, -26, -56, -39, 14, 11, -13, +-28, -69, -89, -74, -83, -45, 1, 0, 16, 45, +-3, 5, 91, 104, 119, 111, 97, 76, 72, 39, +19, 5, -12, -34, -41, -72, -67, 14, 21, 1, +-17, -49, 15, 27, 13, 6, 2, 12, 27, 10, +-10, -30, 18, 36, 93, 87, 87, 69, 58, 26, +37, 32, 64, 54, 53, 67, 78, 73, 80, 48, +48, 46, 38, 9, -1, 8, 47, 44, 58, 33, +-23, -61, -88, -76, 27, 14, 9, 36, 36, 10, +-24, -66, -88, -97, -55, 20, 19, 16, 52, 31, +4, -16, -1, -8, 2, 0, 15, 13, 29, 14, +-9, -44, -66, -72, 8, 32, 37, 38, 46, 16, +-22, -65, -62, 30, 22, 51, 57, 45, 68, 38, +-4, -42, -53, -33, -54, -53, -4, -18, 62, 40, +-5, -37, -61, -41, 35, 15, -2, -3, -12, -28, +-18, -65, -55, 90, 80, 49, 44, 21, 59, 36, +5, -16, -31, -32, -37, -62, -19, -11, 9, 16, +-22, -60, -67, -51, -61, -35, -5, -18, -27, 24, +-18, -55, -39, -42, -40, 20, 25, 6, 6, 7, +0, -25, -42, -50, -59, -84, -31, -13, -22, 17, +-32, -73, -100, -89, -21, -10, 18, 38, 31, 23, +-15, -54, -57, -13, -18, -41, -32, 17, 50, 21, +-16, -57, -71, -10, -8, -26, -38, -47, 42, 25, +-17, -58, -82, -7, 33, 3, 30, 11, 13, 24, +-23, -61, -97, -83, 82, 81, 57, 39, 31, 2, +26, 32, 104, 86, 62, 55, 77, 70, 74, 43, +-8, -29, -33, -52, -74, -73, -17, 14, 39, 25, +-21, -60, -68, -22, 43, 37, 51, 54, 64, 36, +-5, -30, -50, -52, 22, 7, 10, 13, 0, -20, +-15, -53, -61, -34, -50, -25, 15, 6, 0, -14, +-10, -51, -60, 20, 77, 50, 34, 22, 8, -8, +-8, -45, -52, -59, -76, -35, -43, -49, 47, 40, +41, 44, 53, 40, 41, 30, 38, 31, 46, 24, +19, 3, -15, -26, -16, -45, -13, 24, 18, -5, +-3, -39, -54, -35, -49, -70, 2, -6, -19, 15, +-13, -53, -36, 6, -3, 45, 50, 30, 33, 10, +15, -25, -16, -4, -32, 7, 23, 6, 67, 36, +-21, -58, -87, -71, 10, -12, -16, 19, 10, 0, +9, -15, -6, 25, 31, 7, 30, 28, 20, 2, +2, -28, -42, -50, -67, -75, 12, 12, -4, 8, +-17, -63, -58, 56, 49, 36, 60, 38, 37, 14, +-13, -48, -38, -29, -44, -25, -20, -33, 38, 20, +2, -37, -42, -1, -18, -42, 3, -17, -19, 12, +-20, -60, -73, -7, 12, -13, -22, -38, 4, 10, +-8, -41, -63, -65, 42, 63, 45, 31, 31, 0, +-4, -46, -38, -9, -35, 15, 50, 27, 67, 39, +3, -7, -13, -34, -51, -47, 9, 39, 54, 29, +-29, -71, -89, -52, -39, -6, 3, 2, 38, 28, +-14, -42, -13, -15, -19, 0, 20, 12, 37, 25, +11, -26, -24, -40, -65, -11, -13, -27, 65, 37, +0, -31, -46, -17, -21, -45, 11, -2, 25, 24, +51, 91, 102, 87, 85, 63, 57, 42, 48, 23, +-4, -34, -56, -70, -10, -14, 4, 18, 3, -8, +-23, -59, -39, -44, -42, -14, -10, -23, -4, 17, +-2, 4, 35, 63, 69, 75, 82, 63, 78, 48, +-17, -55, -60, -3, -4, -19, 4, -4, -18, -29, +2, -20, -36, -50, -32, -59, -16, 30, 16, -6, +-12, -47, 24, 68, 45, 46, 41, 32, 65, 37, +-4, -40, -54, -67, -96, -66, -9, -25, 42, 38, +-15, 13, 58, 58, 84, 104, 119, 104, 89, 51, +-15, -24, -9, -24, -27, -50, -7, 28, 29, 6, +-7, -33, -35, -49, -65, -53, -37, -10, 33, 14, +31, 19, 46, 72, 67, 45, 83, 68, 63, 41, +-14, -53, -59, -17, 55, 79, 64, 39, 43, 10, +}; + +const signed char cdbk_lsp2_vlbr[160]={ +-20, -30, -24, 17, 7, -13, -21, 61, 56, 16, +12, 1, 10, 77, 32, 3, 7, 3, -25, -31, +-4, 2, -36, -83, 18, 5, -5, 5, 11, 23, +-2, -1, -11, -12, -20, -28, 68, 50, -17, -20, +5, 2, 1, 20, 17, 4, -52, -66, 36, 24, +-4, -10, 7, -15, -32, 80, 37, 8, -13, -29, +33, 37, 28, 15, 8, 14, 35, 18, 50, 36, +-4, -1, 4, -7, 3, 3, -11, -58, -75, 13, +13, 21, 24, -11, -12, -38, -72, 33, 15, -12, +-44, -17, 83, 21, 2, 7, 0, 4, 0, -1, +-25, -42, -51, 33, 20, 15, 30, -13, 9, 32, +6, 2, -8, 7, -38, -77, 6, -13, -7, 32, +48, 57, 32, -12, -10, -4, 2, -15, -29, -29, +2, 10, -9, -16, 79, 44, 7, 12, -5, -18, +-23, -29, -35, -3, -3, -18, -34, -3, -39, -50, +-5, -10, -8, -37, -76, 11, -4, -19, 30, 16, +}; + +#endif diff --git a/pjmedia/src/pjmedia-codec/speex/lpc.c b/pjmedia/src/pjmedia-codec/speex/lpc.c deleted file mode 100644 index fd5d3821..00000000 --- a/pjmedia/src/pjmedia-codec/speex/lpc.c +++ /dev/null @@ -1,201 +0,0 @@ -/* - Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann, - Technische Universitaet Berlin - - Any use of this software is permitted provided that this notice is not - removed and that neither the authors nor the Technische Universitaet Berlin - are deemed to have made any representations as to the suitability of this - software for any purpose nor are held responsible for any defects of - this software. THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE. - - As a matter of courtesy, the authors request to be informed about uses - this software has found, about bugs in this software, and about any - improvements that may be of general interest. - - Berlin, 28.11.1994 - Jutta Degener - Carsten Bormann - - - Code modified by Jean-Marc Valin - - Speex License: - - 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. -*/ - -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -#include "lpc.h" - -#ifdef BFIN_ASM -#include "lpc_bfin.h" -#endif - -/* LPC analysis - * - * The next two functions calculate linear prediction coefficients - * and/or the related reflection coefficients from the first P_MAX+1 - * values of the autocorrelation function. - */ - -/* Invented by N. Levinson in 1947, modified by J. Durbin in 1959. - */ - -/* returns minimum mean square error */ -spx_word32_t _spx_lpc( -spx_coef_t *lpc, /* out: [0...p-1] LPC coefficients */ -const spx_word16_t *ac, /* in: [0...p] autocorrelation values */ -int p -) -{ - int i, j; - spx_word16_t r; - spx_word16_t error = ac[0]; - - if (ac[0] == 0) - { - for (i = 0; i < p; i++) - lpc[i] = 0; - return 0; - } - - for (i = 0; i < p; i++) { - - /* Sum up this iteration's reflection coefficient */ - spx_word32_t rr = NEG32(SHL32(EXTEND32(ac[i + 1]),13)); - for (j = 0; j < i; j++) - rr = SUB32(rr,MULT16_16(lpc[j],ac[i - j])); -#ifdef FIXED_POINT - r = DIV32_16(rr+PSHR32(error,1),ADD16(error,8)); -#else - r = rr/(error+.003*ac[0]); -#endif - /* Update LPC coefficients and total error */ - lpc[i] = r; - for (j = 0; j < i>>1; j++) - { - spx_word16_t tmp = lpc[j]; - lpc[j] = MAC16_16_P13(lpc[j],r,lpc[i-1-j]); - lpc[i-1-j] = MAC16_16_P13(lpc[i-1-j],r,tmp); - } - if (i & 1) - lpc[j] = MAC16_16_P13(lpc[j],lpc[j],r); - - error = SUB16(error,MULT16_16_Q13(r,MULT16_16_Q13(error,r))); - } - return error; -} - - -#ifdef FIXED_POINT - -/* Compute the autocorrelation - * ,--, - * ac(i) = > x(n) * x(n-i) for all n - * `--' - * for lags between 0 and lag-1, and x == 0 outside 0...n-1 - */ - -#ifndef OVERRIDE_SPEEX_AUTOCORR -void _spx_autocorr( -const spx_word16_t *x, /* in: [0...n-1] samples x */ -spx_word16_t *ac, /* out: [0...lag-1] ac values */ -int lag, -int n -) -{ - spx_word32_t d; - int i, j; - spx_word32_t ac0=1; - int shift, ac_shift; - - for (j=0;j x(n) * x(n-i) for all n - * `--' - * for lags between 0 and lag-1, and x == 0 outside 0...n-1 - */ -void _spx_autocorr( -const spx_word16_t *x, /* in: [0...n-1] samples x */ -float *ac, /* out: [0...lag-1] ac values */ -int lag, -int n -) -{ - float d; - int i; - while (lag--) - { - for (i = lag, d = 0; i < n; i++) - d += x[i] * x[i-lag]; - ac[lag] = d; - } - ac[0] += 10; -} - -#endif - - diff --git a/pjmedia/src/pjmedia-codec/speex/lpc_spx.c b/pjmedia/src/pjmedia-codec/speex/lpc_spx.c index c465faea..fd5d3821 100644 --- a/pjmedia/src/pjmedia-codec/speex/lpc_spx.c +++ b/pjmedia/src/pjmedia-codec/speex/lpc_spx.c @@ -94,7 +94,7 @@ int p for (j = 0; j < i; j++) rr = SUB32(rr,MULT16_16(lpc[j],ac[i - j])); #ifdef FIXED_POINT - r = DIV32_16(rr,ADD16(error,16)); + r = DIV32_16(rr+PSHR32(error,1),ADD16(error,8)); #else r = rr/(error+.003*ac[0]); #endif @@ -103,11 +103,11 @@ int p for (j = 0; j < i>>1; j++) { spx_word16_t tmp = lpc[j]; - lpc[j] = MAC16_16_Q13(lpc[j],r,lpc[i-1-j]); - lpc[i-1-j] = MAC16_16_Q13(lpc[i-1-j],r,tmp); + lpc[j] = MAC16_16_P13(lpc[j],r,lpc[i-1-j]); + lpc[i-1-j] = MAC16_16_P13(lpc[i-1-j],r,tmp); } if (i & 1) - lpc[j] = MAC16_16_Q13(lpc[j],lpc[j],r); + lpc[j] = MAC16_16_P13(lpc[j],lpc[j],r); error = SUB16(error,MULT16_16_Q13(r,MULT16_16_Q13(error,r))); } diff --git a/pjmedia/src/pjmedia-codec/speex/lsp.c b/pjmedia/src/pjmedia-codec/speex/lsp.c index 6e7ea311..a73d8835 100644 --- a/pjmedia/src/pjmedia-codec/speex/lsp.c +++ b/pjmedia/src/pjmedia-codec/speex/lsp.c @@ -4,8 +4,8 @@ Original copyright AUTHOR......: David Rowe DATE CREATED: 24/2/93 -Heavily modified by Jean-Marc Valin (fixed-point, optimizations, - additional functions, ...) +Heavily modified by Jean-Marc Valin (c) 2002-2006 (fixed-point, + optimizations, additional functions, ...) This file contains functions for converting Linear Prediction Coefficients (LPC) to Line Spectral Pair (LSP) and back. Note that the @@ -509,7 +509,7 @@ void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack) /* hard limit ak's to +/- 32767 */ - if (a < -32767) a = 32767; + if (a < -32767) a = -32767; if (a > 32767) a = 32767; ak[j-1] = (short)a; diff --git a/pjmedia/src/pjmedia-codec/speex/lsp_bfin.h b/pjmedia/src/pjmedia-codec/speex/lsp_bfin.h new file mode 100644 index 00000000..20e50528 --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/lsp_bfin.h @@ -0,0 +1,89 @@ +/* Copyright (C) 2006 David Rowe */ +/** + @file lsp_bfin.h + @author David Rowe + @brief LSP routines optimised for the Blackfin +*/ +/* + 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 OVERRIDE_CHEB_POLY_EVA +#ifdef OVERRIDE_CHEB_POLY_EVA +static inline spx_word32_t cheb_poly_eva( + spx_word16_t *coef, /* P or Q coefs in Q13 format */ + spx_word16_t x, /* cos of freq (-1.0 to 1.0) in Q14 format */ + int m, /* LPC order/2 */ + char *stack +) +{ + spx_word32_t sum; + + __asm__ __volatile__ + ( + "P0 = %2;\n\t" /* P0: coef[m], coef[m-1],..., coef[0] */ + "R4 = 8192;\n\t" /* R4: rounding constant */ + "R2 = %1;\n\t" /* R2: x */ + + "R5 = -16383;\n\t" + "R2 = MAX(R2,R5);\n\t" + "R5 = 16383;\n\t" + "R2 = MIN(R2,R5);\n\t" + + "R3 = W[P0--] (X);\n\t" /* R3: sum */ + "R5 = W[P0--] (X);\n\t" + "R5 = R5.L * R2.L (IS);\n\t" + "R5 = R5 + R4;\n\t" + "R5 >>>= 14;\n\t" + "R3 = R3 + R5;\n\t" + + "R0 = R2;\n\t" /* R0: b0 */ + "R1 = 16384;\n\t" /* R1: b1 */ + "LOOP cpe%= LC0 = %3;\n\t" + "LOOP_BEGIN cpe%=;\n\t" + "P1 = R0;\n\t" + "R0 = R2.L * R0.L (IS) || R5 = W[P0--] (X);\n\t" + "R0 >>>= 13;\n\t" + "R0 = R0 - R1;\n\t" + "R1 = P1;\n\t" + "R5 = R5.L * R0.L (IS);\n\t" + "R5 = R5 + R4;\n\t" + "R5 >>>= 14;\n\t" + "R3 = R3 + R5;\n\t" + "LOOP_END cpe%=;\n\t" + "%0 = R3;\n\t" + : "=&d" (sum) + : "a" (x), "a" (&coef[m]), "a" (m-1) + : "R0", "R1", "R3", "R2", "R4", "R5", "P0", "P1" + ); + return sum; +} +#endif + + + diff --git a/pjmedia/src/pjmedia-codec/speex/ltp.c b/pjmedia/src/pjmedia-codec/speex/ltp.c index 009923de..1801bd25 100644 --- a/pjmedia/src/pjmedia-codec/speex/ltp.c +++ b/pjmedia/src/pjmedia-codec/speex/ltp.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin +/* Copyright (C) 2002-2006 Jean-Marc Valin File: ltp.c Long-Term Prediction functions @@ -152,7 +152,7 @@ void pitch_xcorr(const spx_word16_t *_x, const spx_word16_t *_y, spx_word32_t *c #endif #ifndef OVERRIDE_COMPUTE_PITCH_ERROR -PJ_INLINE(spx_word32_t) compute_pitch_error(spx_word16_t *C, spx_word16_t *g, spx_word16_t pitch_control) +static inline spx_word32_t compute_pitch_error(spx_word16_t *C, spx_word16_t *g, spx_word16_t pitch_control) { spx_word32_t sum = 0; sum = ADD32(sum,MULT16_16(MULT16_16_16(g[0],pitch_control),C[0])); @@ -177,7 +177,10 @@ void open_loop_nbest_pitch(spx_word16_t *sw, int start, int end, int len, int *p spx_word32_t e0; VARDECL(spx_word32_t *corr); VARDECL(spx_word32_t *energy); - +#ifdef FIXED_POINT + int scaledown = 0; +#endif + ALLOC(best_score, N, spx_word32_t); ALLOC(best_ener, N, spx_word32_t); ALLOC(corr, end-start+1, spx_word32_t); @@ -189,7 +192,24 @@ void open_loop_nbest_pitch(spx_word16_t *sw, int start, int end, int len, int *p best_ener[i]=0; pitch[i]=start; } - +#ifdef FIXED_POINT + for (i=-end;i16383) + { + scaledown=1; + break; + } + } + /* If the weighted input is close to saturation, then we scale it down */ + if (scaledown) + { + for (i=-end;i=0;i--) { spx_word16_t e0=exc2[-pitch-1+i]; +#ifdef FIXED_POINT + /* Scale excitation down if needed (avoiding overflow) */ + if (scaledown) + e0 = SHR16(e0,1); +#endif x[i][0]=MULT16_16_Q14(r[0], e0); for (j=0;j16383) + { + scaledown=1; + break; + } + } + for (i=-end;i16383) + { + scaledown=1; + break; + } + } +#endif if (N>end-start+1) N=end-start+1; if (end != start) @@ -562,7 +631,7 @@ spx_word32_t *cumul_gain for (j=0;j16777216) + int r=0; + if (x>=(spx_int32_t)65536) { - x>>=10; - k+=5; + x >>= 16; + r += 16; } - if (x>1048576) + if (x>=256) { - x>>=6; - k+=3; + x >>= 8; + r += 8; } - if (x>262144) + if (x>=16) { - x>>=4; - k+=2; + x >>= 4; + r += 4; } - if (x>32768) + if (x>=4) { - x>>=2; - k+=1; + x >>= 2; + r += 2; } - if (x>16384) + if (x>=2) { - x>>=2; - k+=1; + r += 1; } -#else - while (x>16384) + return r; +} + +spx_int16_t spx_ilog4(spx_uint32_t x) +{ + int r=0; + if (x>=(spx_int32_t)65536) { - x>>=2; - k++; - } -#endif - while (x<4096) + x >>= 16; + r += 8; + } + if (x>=256) { - x<<=2; - k--; + x >>= 8; + r += 4; } + if (x>=16) + { + x >>= 4; + r += 2; + } + if (x>=4) + { + r += 1; + } + return r; +} + +#ifdef FIXED_POINT + +/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25723*x^3 (for .25 < x < 1) */ +/*#define C0 3634 +#define C1 21173 +#define C2 -12627 +#define C3 4215*/ + +/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25659*x^3 (for .25 < x < 1) */ +#define C0 3634 +#define C1 21173 +#define C2 -12627 +#define C3 4204 + +spx_word16_t spx_sqrt(spx_word32_t x) +{ + int k; + spx_word32_t rt; + k = spx_ilog4(x)-6; + x = VSHR32(x, (k<<1)); rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3))))))); - if (k>0) - rt <<= k; - else - rt >>= -k; - rt >>=7; + rt = VSHR32(rt,7-k); return rt; } @@ -149,6 +167,101 @@ spx_word16_t spx_cos(spx_word16_t x) } } +#define L1 32767 +#define L2 -7651 +#define L3 8277 +#define L4 -626 + +static inline spx_word16_t _spx_cos_pi_2(spx_word16_t x) +{ + 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)))))))); +} + +spx_word16_t spx_cos_norm(spx_word32_t x) +{ + x = x&0x0001ffff; + if (x>SHL32(EXTEND32(1), 16)) + x = SUB32(SHL32(EXTEND32(1), 17),x); + if (x&0x00007fff) + { + if (x14) + return 0x7fffffff; + else if (integer < -15) + return 0; + frac = SHL16(x-SHL16(integer,11),3); + frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac)))))); + return VSHR32(EXTEND32(frac), -integer-2); +} + +/* Input in Q11 format, output in Q16 */ +spx_word32_t spx_exp(spx_word16_t x) +{ + if (x>21290) + return 0x7fffffff; + else if (x<-21290) + return 0; + else + return spx_exp2(MULT16_16_P14(23637,x)); +} +#define M1 32767 +#define M2 -21 +#define M3 -11943 +#define M4 4936 + +static inline spx_word16_t spx_atan01(spx_word16_t x) +{ + return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x))))))); +} + +/* Input in Q15, output in Q14 */ +spx_word16_t spx_atan(spx_word32_t x) +{ + if (x <= 32767) + { + return SHR16(spx_atan01(x),1); + } else { + int e = spx_ilog2(x); + if (e>=29) + return 25736; + x = DIV32_16(SHL32(EXTEND32(32767),29-e), EXTRACT16(SHR32(x, e-14))); + return SUB16(25736, SHR16(spx_atan01(x),1)); + } +} #else #ifndef M_PI @@ -175,5 +288,4 @@ spx_word16_t spx_cos(spx_word16_t x) } } - #endif diff --git a/pjmedia/src/pjmedia-codec/speex/math_approx.h b/pjmedia/src/pjmedia-codec/speex/math_approx.h index 377bf1ac..49cfda6e 100644 --- a/pjmedia/src/pjmedia-codec/speex/math_approx.h +++ b/pjmedia/src/pjmedia-codec/speex/math_approx.h @@ -38,13 +38,25 @@ #include "misc.h" spx_word16_t spx_cos(spx_word16_t x); - +spx_int16_t spx_ilog2(spx_uint32_t x); +spx_int16_t spx_ilog4(spx_uint32_t x); #ifdef FIXED_POINT spx_word16_t spx_sqrt(spx_word32_t x); spx_word16_t spx_acos(spx_word16_t x); +spx_word32_t spx_exp(spx_word16_t x); +spx_word16_t spx_cos_norm(spx_word32_t x); + +/* Input in Q15, output in Q14 */ +spx_word16_t spx_atan(spx_word32_t x); + #else + #define spx_sqrt sqrt #define spx_acos acos +#define spx_exp exp +#define spx_cos_norm(x) (cos((.5f*M_PI)*(x))) +#define spx_atan atan + #endif #endif diff --git a/pjmedia/src/pjmedia-codec/speex/mdf.c b/pjmedia/src/pjmedia-codec/speex/mdf.c index 94e08f7f..3cea2343 100644 --- a/pjmedia/src/pjmedia-codec/speex/mdf.c +++ b/pjmedia/src/pjmedia-codec/speex/mdf.c @@ -41,8 +41,9 @@ double-talk is achieved using a variable learning rate as described in: Valin, J.-M., On Adjusting the Learning Rate in Frequency Domain Echo - Cancellation With Double-Talk. Submitted to IEEE Transactions on Speech - and Audio Processing, 2006. + Cancellation With Double-Talk. To appear in IEEE Transactions on Audio, + Speech and Language Processing, 2006. + 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 @@ -78,9 +79,6 @@ #define M_PI 3.14159265358979323846 #endif -#define min(a,b) ((a)<(b) ? (a) : (b)) -#define max(a,b) ((a)>(b) ? (a) : (b)) - #ifdef FIXED_POINT #define WEIGHT_SHIFT 11 #define NORMALIZE_SCALEDOWN 5 @@ -89,15 +87,24 @@ #define WEIGHT_SHIFT 0 #endif +/* If enabled, the transition between blocks is smooth, so there isn't any blocking +aftifact when adapting. The cost is an extra FFT and a matrix-vector multiply */ +#define SMOOTH_BLOCKS + #ifdef FIXED_POINT -static const spx_float_t MIN_LEAK = {16777, -24}; +static const spx_float_t MIN_LEAK = {16777, -19}; #define TOP16(x) ((x)>>16) #else -static const spx_float_t MIN_LEAK = .001f; +static const spx_float_t MIN_LEAK = .032f; #define TOP16(x) (x) #endif +#define PLAYBACK_DELAY 2 + +void speex_echo_get_residual(SpeexEchoState *st, spx_word32_t *Yout, int len); + + /** Speex echo cancellation state. */ struct SpeexEchoState_ { int frame_size; /**< Number of samples processed each time */ @@ -105,36 +112,40 @@ struct SpeexEchoState_ { int M; int cancel_count; int adapted; + int saturated; + int screwed_up; spx_int32_t sampling_rate; spx_word16_t spec_average; spx_word16_t beta0; spx_word16_t beta_max; spx_word32_t sum_adapt; - spx_word16_t *e; + spx_word16_t leak_estimate; + + spx_word16_t *e; /* scratch */ spx_word16_t *x; spx_word16_t *X; - spx_word16_t *d; - spx_word16_t *y; + spx_word16_t *input; /* scratch */ + spx_word16_t *y; /* scratch */ spx_word16_t *last_y; - spx_word32_t *Yps; - spx_word16_t *Y; + spx_word16_t *Y; /* scratch */ spx_word16_t *E; - spx_word32_t *PHI; + spx_word32_t *PHI; /* scratch */ spx_word32_t *W; spx_word32_t *power; spx_float_t *power_1; - spx_word16_t *wtmp; + spx_word16_t *wtmp; /* scratch */ #ifdef FIXED_POINT - spx_word16_t *wtmp2; + spx_word16_t *wtmp2; /* scratch */ #endif - spx_word32_t *Rf; - spx_word32_t *Yf; - spx_word32_t *Xf; + spx_word32_t *Rf; /* scratch */ + spx_word32_t *Yf; /* scratch */ + spx_word32_t *Xf; /* scratch */ spx_word32_t *Eh; spx_word32_t *Yh; spx_float_t Pey; spx_float_t Pyy; spx_word16_t *window; + spx_word16_t *prop; void *fft_table; spx_word16_t memX, memD, memE; spx_word16_t preemph; @@ -144,9 +155,10 @@ struct SpeexEchoState_ { /* NOTE: If you only use speex_echo_cancel() and want to save some memory, remove this */ spx_int16_t *play_buf; int play_buf_pos; + int play_buf_started; }; -PJ_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 i; spx_word16_t den2; @@ -170,7 +182,7 @@ PJ_INLINE(void) filter_dc_notch16(const spx_int16_t *in, spx_word16_t radius, sp } } -PJ_INLINE(spx_word32_t ) mdf_inner_prod(const spx_word16_t *x, const spx_word16_t *y, int len) +static inline spx_word32_t mdf_inner_prod(const spx_word16_t *x, const spx_word16_t *y, int len) { spx_word32_t sum=0; len >>= 1; @@ -186,7 +198,7 @@ PJ_INLINE(spx_word32_t ) mdf_inner_prod(const spx_word16_t *x, const spx_word16_ } /** Compute power spectrum of a half-complex (packed) vector */ -PJ_INLINE(void) power_spectrum(const spx_word16_t *X, spx_word32_t *ps, int N) +static inline void power_spectrum(const spx_word16_t *X, spx_word32_t *ps, int N) { int i, j; ps[0]=MULT16_16(X[0],X[0]); @@ -227,7 +239,7 @@ static inline void spectral_mul_accum(const spx_word16_t *X, const spx_word32_t acc[N-1] = PSHR32(tmp1,WEIGHT_SHIFT); } #else -PJ_INLINE(void) spectral_mul_accum(const spx_word16_t *X, const spx_word32_t *Y, spx_word16_t *acc, int N, int M) +static inline void spectral_mul_accum(const spx_word16_t *X, const spx_word32_t *Y, spx_word16_t *acc, int N, int M) { int i,j; for (i=0;iM = (filter_length+st->frame_size-1)/frame_size; st->cancel_count=0; st->sum_adapt = 0; - /* FIXME: Make that an init option (new API call?) */ + st->saturated = 0; + st->screwed_up = 0; + /* This is the default sampling rate */ st->sampling_rate = 8000; st->spec_average = DIV32_16(SHL32(EXTEND32(st->frame_size), 15), st->sampling_rate); #ifdef FIXED_POINT @@ -283,14 +301,14 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) st->beta0 = (2.0f*st->frame_size)/st->sampling_rate; st->beta_max = (.5f*st->frame_size)/st->sampling_rate; #endif + st->leak_estimate = 0; 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->d = (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->Yps = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); st->last_y = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); st->Yf = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t)); st->Rf = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t)); @@ -298,14 +316,15 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) st->Yh = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t)); st->Eh = (spx_word32_t*)speex_alloc((st->frame_size+1)*sizeof(spx_word32_t)); - st->X = (spx_word16_t*)speex_alloc(M*N*sizeof(spx_word16_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->PHI = (spx_word32_t*)speex_alloc(M*N*sizeof(spx_word32_t)); + st->PHI = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); st->power = (spx_word32_t*)speex_alloc((frame_size+1)*sizeof(spx_word32_t)); st->power_1 = (spx_float_t*)speex_alloc((frame_size+1)*sizeof(spx_float_t)); st->window = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); + st->prop = (spx_word16_t*)speex_alloc(M*sizeof(spx_word16_t)); st->wtmp = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); #ifdef FIXED_POINT st->wtmp2 = (spx_word16_t*)speex_alloc(N*sizeof(spx_word16_t)); @@ -318,10 +337,27 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) for (i=0;iwindow[i] = .5-.5*cos(2*M_PI*i/N); #endif + for (i=0;i<=st->frame_size;i++) + st->power_1[i] = FLOAT_ONE; for (i=0;iW[i] = 0; { - st->W[i] = st->PHI[i] = 0; + spx_word32_t sum = 0; + /* Ratio of ~10 between adaptation rate of first and last block */ + spx_word16_t decay = SHR32(spx_exp(NEG16(DIV32_16(QCONST16(2.4,11),M))),1); + st->prop[0] = QCONST16(.7, 15); + sum = EXTEND32(st->prop[0]); + for (i=1;iprop[i] = MULT16_16_Q15(st->prop[i-1], decay); + sum = ADD32(sum, EXTEND32(st->prop[i])); + } + 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->preemph = QCONST16(.9,15); if (st->sampling_rate<12000) @@ -335,9 +371,10 @@ SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length) st->adapted = 0; st->Pey = st->Pyy = FLOAT_ONE; - st->play_buf = (spx_int16_t*)speex_alloc(2*st->frame_size*sizeof(spx_int16_t)); - st->play_buf_pos = 0; - + st->play_buf = (spx_int16_t*)speex_alloc((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; } @@ -346,19 +383,40 @@ void speex_echo_state_reset(SpeexEchoState *st) { int i, M, N; st->cancel_count=0; + st->screwed_up = 0; N = st->window_size; M = st->M; for (i=0;iW[i] = 0; + for (i=0;iX[i] = 0; - } for (i=0;i<=st->frame_size;i++) + { st->power[i] = 0; - + st->power_1[i] = FLOAT_ONE; + st->Eh[i] = 0; + st->Yh[i] = 0; + } + for (i=0;iframe_size;i++) + { + st->last_y[i] = 0; + } + for (i=0;iE[i] = 0; + st->x[i] = 0; + } + st->notch_mem[0] = st->notch_mem[1] = 0; + st->memX=st->memD=st->memE=0; + + st->saturated = 0; st->adapted = 0; st->sum_adapt = 0; st->Pey = st->Pyy = FLOAT_ONE; + for (i=0;i<3*st->frame_size;i++) + st->play_buf[i] = 0; + st->play_buf_pos = PLAYBACK_DELAY*st->frame_size; + st->play_buf_started = 0; } @@ -369,10 +427,9 @@ void speex_echo_state_destroy(SpeexEchoState *st) speex_free(st->e); speex_free(st->x); - speex_free(st->d); + speex_free(st->input); speex_free(st->y); speex_free(st->last_y); - speex_free(st->Yps); speex_free(st->Yf); speex_free(st->Rf); speex_free(st->Xf); @@ -387,6 +444,7 @@ void speex_echo_state_destroy(SpeexEchoState *st) speex_free(st->power); speex_free(st->power_1); speex_free(st->window); + speex_free(st->prop); speex_free(st->wtmp); #ifdef FIXED_POINT speex_free(st->wtmp2); @@ -395,17 +453,19 @@ void speex_echo_state_destroy(SpeexEchoState *st) speex_free(st); } -void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out, spx_int32_t *Yout) +void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out) { int i; + /*speex_warning_int("capture with fill level ", st->play_buf_pos/st->frame_size);*/ + st->play_buf_started = 1; if (st->play_buf_pos>=st->frame_size) { - speex_echo_cancel(st, rec, st->play_buf, out, Yout); + speex_echo_cancellation(st, rec, st->play_buf, out); st->play_buf_pos -= st->frame_size; - for (i=0;iframe_size;i++) + for (i=0;iplay_buf_pos;i++) st->play_buf[i] = st->play_buf[i+st->frame_size]; } else { - speex_warning("no playback frame available"); + speex_warning("No playback frame available (your application is buggy and/or got xruns)"); if (st->play_buf_pos!=0) { speex_warning("internal playback buffer corruption?"); @@ -418,31 +478,48 @@ void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t void speex_echo_playback(SpeexEchoState *st, const spx_int16_t *play) { - if (st->play_buf_pos<=st->frame_size) + /*speex_warning_int("playback with fill level ", st->play_buf_pos/st->frame_size);*/ + if (!st->play_buf_started) + { + speex_warning("discarded first playback frame"); + return; + } + if (st->play_buf_pos<=PLAYBACK_DELAY*st->frame_size) { int i; for (i=0;iframe_size;i++) st->play_buf[st->play_buf_pos+i] = play[i]; st->play_buf_pos += st->frame_size; + if (st->play_buf_pos <= (PLAYBACK_DELAY-1)*st->frame_size) + { + speex_warning("Auto-filling the buffer (your application is buggy and/or got xruns)"); + for (i=0;iframe_size;i++) + st->play_buf[st->play_buf_pos+i] = play[i]; + st->play_buf_pos += st->frame_size; + } } else { - speex_warning("had to discard a playback frame"); + speex_warning("Had to discard a playback frame (your application is buggy and/or got xruns)"); } } /** Performs echo cancellation on a frame */ -void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int16_t *echo, spx_int16_t *out, spx_int32_t *Yout) +void 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); +} + +/** Performs echo cancellation on a frame (deprecated, last arg now ignored) */ +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; - spx_word32_t Syy,See; - spx_word16_t leak_estimate; + spx_word32_t Syy,See,Sxx,Sdd; + spx_word32_t Sey; spx_word16_t ss, ss_1; spx_float_t Pey = FLOAT_ONE, Pyy=FLOAT_ONE; spx_float_t alpha, alpha_1; spx_word16_t RER; spx_word32_t tmp32; - spx_word16_t M_1; - int saturated=0; N = st->window_size; M = st->M; @@ -450,82 +527,130 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int #ifdef FIXED_POINT ss=DIV32_16(11469,M); ss_1 = SUB16(32767,ss); - M_1 = DIV32_16(32767,M); #else ss=.35/M; ss_1 = 1-ss; - M_1 = 1.f/M; #endif - filter_dc_notch16(ref, st->notch_radius, st->d, st->frame_size, st->notch_mem); - /* Copy input data to buffer */ + 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;iframe_size;i++) { - spx_word16_t tmp; spx_word32_t tmp32; st->x[i] = st->x[i+st->frame_size]; - tmp32 = SUB32(EXTEND32(echo[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX))); + tmp32 = SUB32(EXTEND32(far_end[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memX))); #ifdef FIXED_POINT - /*FIXME: If saturation occurs here, we need to freeze adaptation for M frames (not just one) */ + /* If saturation occurs here, we need to freeze adaptation for M+1 frames (not just one) */ if (tmp32 > 32767) { tmp32 = 32767; - saturated = 1; - } + st->saturated = M+1; + } if (tmp32 < -32767) { tmp32 = -32767; - saturated = 1; + st->saturated = M+1; } #endif st->x[i+st->frame_size] = EXTRACT16(tmp32); - st->memX = echo[i]; + st->memX = far_end[i]; - tmp = st->d[i]; - st->d[i] = st->d[i+st->frame_size]; - tmp32 = SUB32(EXTEND32(tmp), EXTEND32(MULT16_16_P15(st->preemph, st->memD))); + tmp32 = SUB32(EXTEND32(st->input[i]), EXTEND32(MULT16_16_P15(st->preemph, st->memD))); #ifdef FIXED_POINT if (tmp32 > 32767) { tmp32 = 32767; - saturated = 1; + if (st->saturated == 0) + st->saturated = 1; } if (tmp32 < -32767) { tmp32 = -32767; - saturated = 1; + if (st->saturated == 0) + st->saturated = 1; } #endif - st->d[i+st->frame_size] = tmp32; - st->memD = tmp; + st->memD = st->input[i]; + st->input[i] = tmp32; } /* Shift memory: this could be optimized eventually*/ - for (i=0;iX[i]=st->X[i+N]; + for (j=M-1;j>=0;j--) + { + for (i=0;iX[(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]); + +#ifdef SMOOTH_BLOCKS + spectral_mul_accum(st->X, st->W, st->Y, N, M); + spx_ifft(st->fft_table, st->Y, st->e); +#endif - /* Convert x (echo input) to frequency domain */ - spx_fft(st->fft_table, st->x, &st->X[(M-1)*N]); + /* 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;iW[j*N+i] = ADD32(st->W[j*N+i], st->PHI[i]); + + } + } else { + st->saturated--; + } + /* Update weight to prevent circular convolution (MDF / AUMDF) */ + for (j=0;jcancel_count%(M-1) == j-1) + { +#ifdef FIXED_POINT + for (i=0;iwtmp2[i] = EXTRACT16(PSHR32(st->W[j*N+i],NORMALIZE_SCALEDOWN+16)); + spx_ifft(st->fft_table, st->wtmp2, st->wtmp); + for (i=0;iframe_size;i++) + { + st->wtmp[i]=0; + } + for (i=st->frame_size;iwtmp[i]=SHL16(st->wtmp[i],NORMALIZE_SCALEUP); + } + spx_fft(st->fft_table, st->wtmp, st->wtmp2); + /* The "-1" in the shift is a sort of kludge that trades less efficient update speed for decrease noise */ + for (i=0;iW[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;iwtmp[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); -#if 1 - spectral_mul_accum(st->X, st->PHI, st->Y, N, M); - spx_ifft(st->fft_table, st->Y, st->e); -#endif - + /* Compute error signal (for the output with de-emphasis) */ for (i=0;iframe_size;i++) { spx_word32_t tmp_out; -#if 1 +#ifdef SMOOTH_BLOCKS spx_word16_t y = MULT16_16_Q15(st->window[i+st->frame_size],st->e[i+st->frame_size]) + MULT16_16_Q15(st->window[i],st->y[i+st->frame_size]); - tmp_out = SUB32(EXTEND32(st->d[i+st->frame_size]), EXTEND32(y)); + tmp_out = SUB32(EXTEND32(st->input[i]), EXTEND32(y)); #else - tmp_out = SUB32(EXTEND32(st->d[i+st->frame_size]), EXTEND32(st->y[i+st->frame_size])); + tmp_out = SUB32(EXTEND32(st->input[i]), EXTEND32(st->y[i+st->frame_size])); #endif /* Saturation */ @@ -534,13 +659,14 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int else if (tmp_out<-32768) tmp_out = -32768; tmp_out = ADD32(tmp_out, EXTEND32(MULT16_16_P15(st->preemph, st->memE))); - /* This is an arbitrary test for saturation */ - if (ref[i] <= -32000 || ref[i] >= 32000) + /* This is an arbitrary test for saturation in the microphone signal */ + if (in[i] <= -32000 || in[i] >= 32000) { tmp_out = 0; - saturated = 1; + if (st->saturated == 0) + st->saturated = 1; } - out[i] = tmp_out; + out[i] = (spx_int16_t)tmp_out; st->memE = tmp_out; } @@ -548,26 +674,57 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int for (i=0;iframe_size;i++) { st->e[i] = 0; - st->e[i+st->frame_size] = st->d[i+st->frame_size] - st->y[i+st->frame_size]; + st->e[i+st->frame_size] = st->input[i] - st->y[i+st->frame_size]; } /* Compute a bunch of correlations */ + Sey = mdf_inner_prod(st->e+st->frame_size, st->y+st->frame_size, st->frame_size); See = mdf_inner_prod(st->e+st->frame_size, st->e+st->frame_size, st->frame_size); - See = ADD32(See, SHR32(EXTEND32(10000),6)); Syy = mdf_inner_prod(st->y+st->frame_size, st->y+st->frame_size, st->frame_size); + Sxx = mdf_inner_prod(st->x+st->frame_size, st->x+st->frame_size, st->frame_size); + Sdd = mdf_inner_prod(st->input, st->input, st->frame_size); + /* Do some sanity check */ + if (!(Syy>=0 && Sxx>=0 && See >= 0) +#ifndef FIXED_POINT + || !(See < N*1e9 && Syy < N*1e9 && Sxx < N*1e9) +#endif + ) + { + /* Things have gone really bad */ + st->screwed_up += 50; + for (i=0;iframe_size;i++) + out[i] = 0; + } else if (SHR32(See, 2) > ADD32(Sdd, SHR32(MULT16_16(N, 100),6))) + { + /* AEC seems to add lots of echo instead of removing it, let's see if it will improve */ + st->screwed_up++; + } else { + /* Everything's fine */ + st->screwed_up=0; + } + if (st->screwed_up>=50) + { + speex_warning("The echo canceller started acting funny and got slapped (reset). It swears it will behave now."); + speex_echo_state_reset(st); + return; + } + + /* Add a small noise floor to make sure not to have problems when dividing */ + 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;iframe_size;i++) st->y[i] = 0; spx_fft(st->fft_table, st->y, st->Y); - /* Compute power spectrum of echo (X), error (E) and filter response (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[(M-1)*N], st->Xf, N); + power_spectrum(st->X, st->Xf, N); - /* Smooth echo energy estimate over time */ + /* 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]); @@ -577,7 +734,7 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int { float scale2 = .5f/M; for (j=0;j<=st->frame_size;j++) - st->power[j] = 0; + st->power[j] = 100; for (i=0;iX[i*N], st->Xf, N); @@ -603,6 +760,9 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int #endif } + Pyy = FLOAT_SQRT(Pyy); + Pey = FLOAT_DIVU(Pey,Pyy); + /* Compute correlation updatete rate */ tmp32 = MULT16_32_Q15(st->beta0,Syy); if (tmp32 > MULT16_32_Q15(st->beta_max,See)) @@ -620,29 +780,41 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int if (FLOAT_GT(st->Pey, st->Pyy)) st->Pey = st->Pyy; /* leak_estimate is the linear regression result */ - leak_estimate = FLOAT_EXTRACT16(FLOAT_SHL(FLOAT_DIVU(st->Pey, st->Pyy),14)); + st->leak_estimate = FLOAT_EXTRACT16(FLOAT_SHL(FLOAT_DIVU(st->Pey, st->Pyy),14)); /* This looks like a stupid bug, but it's right (because we convert from Q14 to Q15) */ - if (leak_estimate > 16383) - leak_estimate = 32767; + if (st->leak_estimate > 16383) + st->leak_estimate = 32767; else - leak_estimate = SHL16(leak_estimate,1); - /*printf ("%f\n", leak_estimate);*/ + st->leak_estimate = SHL16(st->leak_estimate,1); + /*printf ("%f\n", st->leak_estimate);*/ /* Compute Residual to Error Ratio */ #ifdef FIXED_POINT - tmp32 = MULT16_32_Q15(leak_estimate,Syy); - tmp32 = ADD32(tmp32, SHL32(tmp32,1)); + tmp32 = MULT16_32_Q15(st->leak_estimate,Syy); + tmp32 = ADD32(SHR32(Sxx,13), ADD32(tmp32, SHL32(tmp32,1))); + /* Check for y in e (lower bound on RER) */ + { + spx_float_t bound = PSEUDOFLOAT(Sey); + bound = FLOAT_DIVU(FLOAT_MULT(bound, bound), PSEUDOFLOAT(ADD32(1,Syy))); + if (FLOAT_GT(bound, PSEUDOFLOAT(See))) + tmp32 = See; + else if (tmp32 < FLOAT_EXTRACT32(bound)) + tmp32 = FLOAT_EXTRACT32(bound); + } if (tmp32 > SHR32(See,1)) tmp32 = SHR32(See,1); RER = FLOAT_EXTRACT16(FLOAT_SHL(FLOAT_DIV32(tmp32,See),15)); #else - RER = 3.*MULT16_32_Q15(leak_estimate,Syy) / See; + RER = (.0001*Sxx + 3.*MULT16_32_Q15(st->leak_estimate,Syy)) / See; + /* Check for y in e (lower bound on RER) */ + if (RER < Sey*Sey/(1+See*Syy)) + RER = Sey*Sey/(1+See*Syy); if (RER > .5) RER = .5; #endif /* We consider that the filter has had minimal adaptation if the following is true*/ - if (!st->adapted && st->sum_adapt > QCONST32(1,15)) + if (!st->adapted && st->sum_adapt > QCONST32(M,15)) { st->adapted = 1; } @@ -653,7 +825,7 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int { spx_word32_t r, e; /* Compute frequency-domain adaptation mask */ - r = MULT16_32_Q15(leak_estimate,SHL32(st->Yf[i],3)); + r = MULT16_32_Q15(st->leak_estimate,SHL32(st->Yf[i],3)); e = SHL32(st->Rf[i],3)+1; #ifdef FIXED_POINT if (r>SHR32(e,1)) @@ -662,27 +834,26 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int if (r>.5*e) r = .5*e; #endif - r = MULT16_32_Q15(QCONST16(.8,15),r) + MULT16_32_Q15(QCONST16(.2,15),(spx_word32_t)(MULT16_32_Q15(RER,e))); + r = MULT16_32_Q15(QCONST16(.7,15),r) + MULT16_32_Q15(QCONST16(.3,15),(spx_word32_t)(MULT16_32_Q15(RER,e))); /*st->power_1[i] = adapt_rate*r/(e*(1+st->power[i]));*/ - st->power_1[i] = FLOAT_SHL(FLOAT_DIV32_FLOAT(MULT16_32_Q15(M_1,r),FLOAT_MUL32U(e,st->power[i]+10)),WEIGHT_SHIFT+16); + st->power_1[i] = FLOAT_SHL(FLOAT_DIV32_FLOAT(r,FLOAT_MUL32U(e,st->power[i]+10)),WEIGHT_SHIFT+16); } } else { - spx_word32_t Sxx; + /* Temporary adaption rate if filter is not yet adapted enough */ spx_word16_t adapt_rate=0; - Sxx = mdf_inner_prod(st->x+st->frame_size, st->x+st->frame_size, st->frame_size); - /* Temporary adaption rate if filter is not adapted correctly */ - - tmp32 = MULT16_32_Q15(QCONST16(.15f, 15), Sxx); + if (Sxx > SHR32(MULT16_16(N, 1000),6)) + { + tmp32 = MULT16_32_Q15(QCONST16(.25f, 15), Sxx); #ifdef FIXED_POINT - if (Sxx > SHR32(See,2)) - Sxx = SHR32(See,2); + if (tmp32 > SHR32(See,2)) + tmp32 = SHR32(See,2); #else - if (Sxx > .25*See) - Sxx = .25*See; + if (tmp32 > .25*See) + tmp32 = .25*See; #endif - adapt_rate = FLOAT_EXTRACT16(FLOAT_SHL(FLOAT_DIV32(MULT16_32_Q15(M_1,Sxx), See),15)); - + adapt_rate = FLOAT_EXTRACT16(FLOAT_SHL(FLOAT_DIV32(tmp32, See),15)); + } for (i=0;i<=st->frame_size;i++) st->power_1[i] = FLOAT_SHL(FLOAT_DIV32(EXTEND32(adapt_rate),ADD32(st->power[i],10)),WEIGHT_SHIFT+1); @@ -690,100 +861,57 @@ void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *ref, const spx_int /* How much have we adapted so far? */ st->sum_adapt = ADD32(st->sum_adapt,adapt_rate); } - /* Compute weight gradient */ - for (j=0;jpower_1, &st->X[j*N], st->E, st->PHI+N*j, N); - } - if (!saturated) + /* Save residual echo so it can be used by the nonlinear processor */ + if (st->adapted) { - /* Gradient descent */ - for (i=0;iW[i] += st->PHI[i]; - /* Old value of W in PHI */ - st->PHI[i] = st->W[i] - st->PHI[i]; - } + /* If the filter is adapted, take the filtered echo */ + for (i=0;iframe_size;i++) + st->last_y[i] = st->last_y[st->frame_size+i]; + for (i=0;iframe_size;i++) + st->last_y[st->frame_size+i] = in[i]-out[i]; + } else { + /* If filter isn't adapted yet, all we can do is take the far end signal directly */ + for (i=0;ilast_y[i] = st->x[i]; } + +} + +/* Compute spectrum of estimated echo for use in an echo post-filter */ +void speex_echo_get_residual(SpeexEchoState *st, spx_word32_t *residual_echo, int len) +{ + int i; + spx_word16_t leak2; + int N; - /* Update weight to prevent circular convolution (MDF / AUMDF) */ - for (j=0;jcancel_count%(M-1) == j) - { -#ifdef FIXED_POINT - for (i=0;iwtmp2[i] = EXTRACT16(PSHR32(st->W[j*N+i],NORMALIZE_SCALEDOWN+16)); - spx_ifft(st->fft_table, st->wtmp2, st->wtmp); - for (i=0;iframe_size;i++) - { - st->wtmp[i]=0; - } - for (i=st->frame_size;iwtmp[i]=SHL16(st->wtmp[i],NORMALIZE_SCALEUP); - } - spx_fft(st->fft_table, st->wtmp, st->wtmp2); - /* The "-1" in the shift is a sort of kludge that trades less efficient update speed for decrease noise */ - for (i=0;iW[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;iwtmp[i]=0; - } - spx_fft(st->fft_table, st->wtmp, &st->W[j*N]); -#endif - } - } + N = st->window_size; - /* Compute spectrum of estimated echo for use in an echo post-filter (if necessary)*/ - if (Yout) - { - spx_word16_t leak2; - if (st->adapted) - { - /* If the filter is adapted, take the filtered echo */ - for (i=0;iframe_size;i++) - st->last_y[i] = st->last_y[st->frame_size+i]; - for (i=0;iframe_size;i++) - st->last_y[st->frame_size+i] = ref[i]-out[i]; - } else { - /* If filter isn't adapted yet, all we can do is take the echo signal directly */ - for (i=0;ilast_y[i] = st->x[i]; - } - - /* Apply hanning window (should pre-compute it)*/ - for (i=0;iy[i] = MULT16_16_Q15(st->window[i],st->last_y[i]); + /* Apply hanning window (should pre-compute it)*/ + for (i=0;iy[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, st->Yps, N); + /* 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 (leak_estimate > 16383) - leak2 = 32767; - else - leak2 = SHL16(leak_estimate, 1); + if (st->leak_estimate > 16383) + leak2 = 32767; + else + leak2 = SHL16(st->leak_estimate, 1); #else - if (leak_estimate>.5) - leak2 = 1; - else - leak2 = 2*leak_estimate; + if (st->leak_estimate>.5) + leak2 = 1; + else + leak2 = 2*st->leak_estimate; #endif - /* Estimate residual echo */ - for (i=0;i<=st->frame_size;i++) - Yout[i] = MULT16_32_Q15(leak2,st->Yps[i]); - } + /* Estimate residual echo */ + 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) { switch(request) diff --git a/pjmedia/src/pjmedia-codec/speex/medfilter.c b/pjmedia/src/pjmedia-codec/speex/medfilter.c new file mode 100644 index 00000000..83872e54 --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/medfilter.c @@ -0,0 +1,97 @@ +/* Copyright (C) 2004 Jean-Marc Valin + File medfilter.c + Median filter + + + 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. + +*/ + +#include "medfilter.h" +#include "misc.h" + +MedianFilter *median_filter_new(int N) +{ + MedianFilter *f = (MedianFilter*)speex_alloc(sizeof(MedianFilter)); + f->N = N; + f->ids = (int*)speex_alloc(sizeof(int)*N); + f->val = (float*)speex_alloc(sizeof(float)*N); + f->filled = 0; + return f; +} + +void median_filter_update(MedianFilter *f, float val) +{ + int i=0; + int insert = 0; + while (insertfilled && f->val[insert] < val) + { + insert++; + } + if (f->filled == f->N) + { + int remove; + for (remove=0;removeN;remove++) + if (f->ids[remove] == 0) + break; + if (insert>remove) + insert--; + if (insert > remove) + { + for (i=remove;ival[i] = f->val[i+1]; + f->ids[i] = f->ids[i+1]; + } + } else if (insert < remove) + { + for (i=remove;i>insert;i--) + { + f->val[i] = f->val[i-1]; + f->ids[i] = f->ids[i-1]; + } + } + for (i=0;ifilled;i++) + f->ids[i]--; + } else { + for (i=f->filled;i>insert;i--) + { + f->val[i] = f->val[i-1]; + f->ids[i] = f->ids[i-1]; + } + f->filled++; + } + f->val[insert]=val; + f->ids[insert]=f->filled-1; +} + +float median_filter_get(MedianFilter *f) +{ + return f->val[f->filled>>1]; +} + diff --git a/pjmedia/src/pjmedia-codec/speex/medfilter.h b/pjmedia/src/pjmedia-codec/speex/medfilter.h new file mode 100644 index 00000000..718869c3 --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/medfilter.h @@ -0,0 +1,51 @@ +/* Copyright (C) 2004 Jean-Marc Valin */ +/** + @file medfilter.h + @brief Median filter +*/ +/* + 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 MEDFILTER_H +#define MEDFILTER_H + +/** Median filter. */ +typedef struct { + int N; + int filled; + int *ids; + float *val; +} MedianFilter; + +MedianFilter *median_filter_new(int N); +void median_filter_update(MedianFilter *f, float val); +float median_filter_get(MedianFilter *f); + +#endif diff --git a/pjmedia/src/pjmedia-codec/speex/misc.h b/pjmedia/src/pjmedia-codec/speex/misc.h index ea9cc387..2878ece5 100644 --- a/pjmedia/src/pjmedia-codec/speex/misc.h +++ b/pjmedia/src/pjmedia-codec/speex/misc.h @@ -40,7 +40,7 @@ #define SPEEX_MINOR_VERSION 1 /**< Minor Speex version. */ #define SPEEX_MICRO_VERSION 13 /**< Micro Speex version. */ #define SPEEX_EXTRA_VERSION "" /**< Extra Speex version. */ -#define SPEEX_VERSION "speex-1.1.13" /**< Speex version string. */ +#define SPEEX_VERSION "speex-1.2beta1" /**< Speex version string. */ #endif /* A couple test to catch stupid option combinations */ @@ -62,7 +62,7 @@ #error I suppose you can have a [ARM4/ARM5E/Blackfin] that has float instructions? #endif #ifdef FIXED_POINT_DEBUG -#error Dont you think enabling fixed-point is a good thing to do if you want to debug that? +#error "Don't you think enabling fixed-point is a good thing to do if you want to debug that?" #endif diff --git a/pjmedia/src/pjmedia-codec/speex/modes.c b/pjmedia/src/pjmedia-codec/speex/modes.c index 1da236db..97e7d1e3 100644 --- a/pjmedia/src/pjmedia-codec/speex/modes.c +++ b/pjmedia/src/pjmedia-codec/speex/modes.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin +/* Copyright (C) 2002-2006 Jean-Marc Valin File: modes.c Describes the different modes of the codec @@ -211,7 +211,7 @@ static const SpeexSubmode nb_submode8 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_nb_ulbr, - QCONST16(.65,15), + QCONST16(.5,15), 79 }; @@ -232,7 +232,7 @@ static const SpeexSubmode nb_submode2 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_nb_vlbr, - QCONST16(.55,15), + QCONST16(.6,15), 119 }; @@ -253,7 +253,7 @@ static const SpeexSubmode nb_submode3 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_nb_lbr, - QCONST16(.45,15), + QCONST16(.55,15), 160 }; @@ -274,7 +274,7 @@ static const SpeexSubmode nb_submode4 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_nb_med, - QCONST16(.35,15), + QCONST16(.45,15), 220 }; @@ -295,7 +295,7 @@ static const SpeexSubmode nb_submode5 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_nb, - QCONST16(.2,15), + QCONST16(.3,15), 300 }; @@ -316,7 +316,7 @@ static const SpeexSubmode nb_submode6 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_sb, - QCONST16(.1,15), + QCONST16(.2,15), 364 }; @@ -337,7 +337,7 @@ static const SpeexSubmode nb_submode7 = { split_cb_search_shape_sign, split_cb_shape_sign_unquant, &split_cb_nb, - -1, + QCONST16(.1,15), 492 }; diff --git a/pjmedia/src/pjmedia-codec/speex/modes.h b/pjmedia/src/pjmedia-codec/speex/modes.h index 9828ee61..6a632402 100644 --- a/pjmedia/src/pjmedia-codec/speex/modes.h +++ b/pjmedia/src/pjmedia-codec/speex/modes.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin */ +/* Copyright (C) 2002-2006 Jean-Marc Valin */ /** @file modes.h @brief Describes the different modes of the codec @@ -69,7 +69,7 @@ typedef void (*innovation_quant_func)(spx_word16_t *, spx_coef_t *, spx_coef_t * spx_sig_t *, spx_word16_t *, SpeexBits *, char *, int, int); /** Innovation unquantization function */ -typedef void (*innovation_unquant_func)(spx_sig_t *, const void *, int, SpeexBits*, char *); +typedef void (*innovation_unquant_func)(spx_sig_t *, const void *, int, SpeexBits*, char *, spx_int32_t *); /** Description of a Speex sub-mode (wither narrowband or wideband */ typedef struct SpeexSubmode { diff --git a/pjmedia/src/pjmedia-codec/speex/nb_celp.c b/pjmedia/src/pjmedia-codec/speex/nb_celp.c index 2c416499..24b6866e 100644 --- a/pjmedia/src/pjmedia-codec/speex/nb_celp.c +++ b/pjmedia/src/pjmedia-codec/speex/nb_celp.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin +/* Copyright (C) 2002-2006 Jean-Marc Valin File: nb_celp.c Redistribution and use in source and binary forms, with or without @@ -87,14 +87,14 @@ const spx_word16_t exc_gain_quant_scal1[2]={11546, 17224}; #else -const float exc_gain_quant_scal3_bound[7]={0.112338, 0.236980, 0.369316, 0.492054, 0.637471, 0.828874, 1.132784}; -const float exc_gain_quant_scal3[8]={0.061130, 0.163546, 0.310413, 0.428220, 0.555887, 0.719055, 0.938694, 1.326874}; -const float exc_gain_quant_scal1_bound[1]={0.87798}; -const float exc_gain_quant_scal1[2]={0.70469, 1.05127}; +const float exc_gain_quant_scal3_bound[7]={0.112338f, 0.236980f, 0.369316f, 0.492054f, 0.637471f, 0.828874f, 1.132784f}; +const float exc_gain_quant_scal3[8]={0.061130f, 0.163546f, 0.310413f, 0.428220f, 0.555887f, 0.719055f, 0.938694f, 1.326874f}; +const float exc_gain_quant_scal1_bound[1]={0.87798f}; +const float exc_gain_quant_scal1[2]={0.70469f, 1.05127f}; -#define LSP_MARGIN .002 -#define LSP_DELTA1 .2 -#define LSP_DELTA2 .05 +#define LSP_MARGIN .002f +#define LSP_DELTA1 .2f +#define LSP_DELTA2 .05f #endif @@ -150,48 +150,48 @@ void *nb_encoder_init(const SpeexMode *m) #ifdef VORBIS_PSYCHO st->psy = vorbis_psy_init(8000, 256); - st->curve = speex_alloc(128*sizeof(float)); - st->old_curve = speex_alloc(128*sizeof(float)); - st->psy_window = speex_alloc(256*sizeof(float)); + st->curve = (float*)speex_alloc(128*sizeof(float)); + st->old_curve = (float*)speex_alloc(128*sizeof(float)); + st->psy_window = (float*)speex_alloc(256*sizeof(float)); #endif st->cumul_gain = 1024; /* Allocating input buffer */ - st->winBuf = speex_alloc((st->windowSize-st->frameSize)*sizeof(spx_word16_t)); + st->winBuf = (spx_word16_t*)speex_alloc((st->windowSize-st->frameSize)*sizeof(spx_word16_t)); /* Allocating excitation buffer */ - st->excBuf = speex_alloc((mode->frameSize+mode->pitchEnd+2)*sizeof(spx_word16_t)); + st->excBuf = (spx_word16_t*)speex_alloc((mode->frameSize+mode->pitchEnd+2)*sizeof(spx_word16_t)); st->exc = st->excBuf + mode->pitchEnd + 2; - st->swBuf = speex_alloc((mode->frameSize+mode->pitchEnd+2)*sizeof(spx_word16_t)); + st->swBuf = (spx_word16_t*)speex_alloc((mode->frameSize+mode->pitchEnd+2)*sizeof(spx_word16_t)); st->sw = st->swBuf + mode->pitchEnd + 2; st->window= lpc_window; /* Create the window for autocorrelation (lag-windowing) */ - st->lagWindow = speex_alloc((st->lpcSize+1)*sizeof(spx_word16_t)); + st->lagWindow = (spx_word16_t*)speex_alloc((st->lpcSize+1)*sizeof(spx_word16_t)); for (i=0;ilpcSize+1;i++) st->lagWindow[i]=16384*exp(-.5*sqr(2*M_PI*st->lag_factor*i)); - st->old_lsp = speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); - st->old_qlsp = speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); + st->old_lsp = (spx_lsp_t*)speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); + st->old_qlsp = (spx_lsp_t*)speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); st->first = 1; for (i=0;ilpcSize;i++) { st->old_lsp[i]=LSP_SCALING*(M_PI*((float)(i+1)))/(st->lpcSize+1); } - st->mem_sp = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->mem_sw = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->mem_sw_whole = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->mem_exc = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->mem_exc2 = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_sp = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_sw = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_sw_whole = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_exc = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_exc2 = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->pi_gain = speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); + st->pi_gain = (spx_word32_t*)speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); st->innov_save = NULL; - st->pitch = speex_alloc((st->nbSubframes)*sizeof(int)); + st->pitch = (int*)speex_alloc((st->nbSubframes)*sizeof(int)); - st->vbr = speex_alloc(sizeof(VBRState)); + st->vbr = (VBRState*)speex_alloc(sizeof(VBRState)); vbr_init(st->vbr); st->vbr_quality = 8; st->vbr_enabled = 0; @@ -205,7 +205,9 @@ void *nb_encoder_init(const SpeexMode *m) st->complexity=2; st->sampling_rate=8000; st->dtx_count=0; - + st->isWideband = 0; + st->highpass_enabled = 1; + #ifdef ENABLE_VALGRIND VALGRIND_MAKE_READABLE(st, (st->stack-(char*)st)); #endif @@ -278,7 +280,7 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) int pitch_half[2]; int ol_pitch_id=0; #endif - spx_word16_t *in = vin; + spx_word16_t *in = (spx_word16_t*)vin; st=(EncState *)state; stack=st->stack; @@ -297,6 +299,9 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) speex_move(st->excBuf, st->excBuf+st->frameSize, (st->max_pitch+2)*sizeof(spx_word16_t)); speex_move(st->swBuf, st->swBuf+st->frameSize, (st->max_pitch+2)*sizeof(spx_word16_t)); + if (st->highpass_enabled) + highpass(in, in, st->frameSize, (st->isWideband?HIGHPASS_WIDEBAND:HIGHPASS_NARROWBAND)|HIGHPASS_INPUT, st->mem_hp); + { VARDECL(spx_word16_t *w_sig); VARDECL(spx_word16_t *autocorr); @@ -485,7 +490,7 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) /* delta_qual*=.1*(3+st->vbr_quality);*/ if (st->vbr_enabled) { - int mode; + spx_int32_t mode; int choice=0; float min_diff=100; mode = 8; @@ -535,7 +540,7 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) if (st->abr_enabled) { - int bitrate; + spx_int32_t bitrate; speex_encoder_ctl(state, SPEEX_GET_BITRATE, &bitrate); st->abr_drift+=(bitrate-st->abr_enabled); st->abr_drift2 = .95*st->abr_drift2 + .05*(bitrate-st->abr_enabled); @@ -833,9 +838,9 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) for (i=0;ilpcSize;i++) st->mem_sw[i]=mem[i]; - /* Compute target signal */ + /* Compute target signal (saturation prevents overflows on clipped input speech) */ for (i=0;isubframeSize;i++) - target[i]=SUB16(sw[i],PSHR32(ringing[i],1)); + target[i]=EXTRACT16(SATURATE(SUB32(sw[i],PSHR32(ringing[i],1)),32767)); /* Reset excitation */ for (i=0;isubframeSize;i++) @@ -903,8 +908,9 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) for (i=0;isubframeSize;i++) innov[i]=0; + /* FIXME: Make sure this is save from overflows (so far so good) */ for (i=0;isubframeSize;i++) - real_exc[i] = SUB16(real_exc[i], PSHR32(exc32[i],SIG_SHIFT-1)); + real_exc[i] = SUB16(real_exc[i], EXTRACT16(PSHR32(exc32[i],SIG_SHIFT-1))); ener = SHL32(EXTEND32(compute_rms16(real_exc, st->subframeSize)),SIG_SHIFT); @@ -955,7 +961,7 @@ int nb_encode(void *state, void *vin, SpeexBits *bits) signal_mul(innov, innov, ener, st->subframeSize); for (i=0;isubframeSize;i++) - exc[i] = EXTRACT16(PSHR32(ADD32(SHL32(exc32[i],1),innov[i]),SIG_SHIFT)); + exc[i] = EXTRACT16(SATURATE32(PSHR32(ADD32(SHL32(exc32[i],1),innov[i]),SIG_SHIFT),32767)); } else { speex_error("No fixed codebook"); } @@ -1076,15 +1082,15 @@ void *nb_decoder_init(const SpeexMode *m) st->lpc_enh_enabled=1; - st->excBuf = speex_alloc((st->frameSize + 2*st->max_pitch + st->subframeSize + 12)*sizeof(spx_word16_t)); + st->excBuf = (spx_word16_t*)speex_alloc((st->frameSize + 2*st->max_pitch + st->subframeSize + 12)*sizeof(spx_word16_t)); st->exc = st->excBuf + 2*st->max_pitch + st->subframeSize + 6; for (i=0;iframeSize + st->max_pitch + 1;i++) st->excBuf[i]=0; - st->interp_qlpc = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->old_qlsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->mem_sp = speex_alloc(st->lpcSize*sizeof(spx_mem_t)); - st->pi_gain = speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); + st->interp_qlpc = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->old_qlsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->mem_sp = (spx_mem_t*)speex_alloc(st->lpcSize*sizeof(spx_mem_t)); + st->pi_gain = (spx_word32_t*)speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); st->last_pitch = 40; st->count_lost=0; st->pitch_gain_buf[0] = st->pitch_gain_buf[1] = st->pitch_gain_buf[2] = 0; @@ -1102,6 +1108,9 @@ void *nb_decoder_init(const SpeexMode *m) st->voc_m1=st->voc_m2=st->voc_mean=0; st->voc_offset=0; st->dtx_enabled=0; + st->isWideband = 0; + st->highpass_enabled = 1; + #ifdef ENABLE_VALGRIND VALGRIND_MAKE_READABLE(st, (st->stack-(char*)st)); #endif @@ -1137,12 +1146,13 @@ const spx_word16_t attenuation[10] = {1., 0.961, 0.852, 0.698, 0.527, 0.368, 0.2 static void nb_decode_lost(DecState *st, spx_word16_t *out, char *stack) { - int i, sub; + int i; int pitch_val; spx_word16_t pitch_gain; spx_word16_t fact; spx_word16_t gain_med; spx_word16_t innov_gain; + spx_word16_t noise_gain; if (st->count_lost<10) fact = attenuation[st->count_lost]; @@ -1163,49 +1173,31 @@ static void nb_decode_lost(DecState *st, spx_word16_t *out, char *stack) if (pitch_gain>.85) pitch_gain=.85; #endif - pitch_gain = MULT16_16_Q15(fact,pitch_gain) + VERY_SMALL; - + /* FIXME: This was rms of innovation (not exc) */ + innov_gain = compute_rms16(st->exc, st->frameSize); + noise_gain = MULT16_16_Q15(innov_gain, MULT16_16_Q15(fact, SUB16(Q15ONE,MULT16_16_Q15(pitch_gain,pitch_gain)))); /* Shift all buffers by one frame */ speex_move(st->excBuf, st->excBuf+st->frameSize, (2*st->max_pitch + st->subframeSize + 12)*sizeof(spx_word16_t)); - for (sub=0;subnbSubframes;sub++) - { - int offset; - spx_word16_t *sp; - spx_word16_t *exc; - /* Offset relative to start of frame */ - offset = st->subframeSize*sub; - /* Original signal */ - sp=out+offset; - /* Excitation */ - exc=st->exc+offset; - /* Excitation after post-filter*/ - - /* Make up a plausible excitation */ - /* FIXME: THIS CAN BE IMPROVED */ - /*if (pitch_gain>.95) - pitch_gain=.95;*/ - - /* FIXME: This was rms of innovation (not exc) */ - innov_gain = compute_rms16(st->exc, st->frameSize); - pitch_val = st->last_pitch + SHR32((spx_int32_t)speex_rand(1+st->count_lost, &st->seed),SIG_SHIFT); - if (pitch_val > st->max_pitch) - pitch_val = st->max_pitch; - if (pitch_val < st->min_pitch) - pitch_val = st->min_pitch; - for (i=0;isubframeSize;i++) - { - /* FIXME: Second term need to be 16-bit */ - exc[i]= MULT16_16_Q15(pitch_gain, (exc[i-pitch_val]+VERY_SMALL)) + - MULT16_16_Q15(fact, MULT16_16_Q15(SHL(Q15ONE,15)-SHL(MULT16_16(pitch_gain,pitch_gain),1),speex_rand(innov_gain, &st->seed))); - } - for (i=0;isubframeSize;i++) - sp[i]=exc[i-st->subframeSize]; - iir_mem16(sp, st->interp_qlpc, sp, st->subframeSize, st->lpcSize, - st->mem_sp, stack); + - bw_lpc(QCONST16(.98,15), st->interp_qlpc, st->interp_qlpc, st->lpcSize); + pitch_val = st->last_pitch + SHR32((spx_int32_t)speex_rand(1+st->count_lost, &st->seed),SIG_SHIFT); + if (pitch_val > st->max_pitch) + pitch_val = st->max_pitch; + if (pitch_val < st->min_pitch) + pitch_val = st->min_pitch; + for (i=0;iframeSize;i++) + { + st->exc[i]= MULT16_16_Q15(pitch_gain, (st->exc[i-pitch_val]+VERY_SMALL)) + + speex_rand(noise_gain, &st->seed); } + + for (i=0;iframeSize;i++) + out[i]=st->exc[i-st->subframeSize]; + bw_lpc(QCONST16(.98,15), st->interp_qlpc, st->interp_qlpc, st->lpcSize); + iir_mem16(out, st->interp_qlpc, out, st->frameSize, st->lpcSize, + st->mem_sp, stack); + highpass(out, out, st->frameSize, HIGHPASS_NARROWBAND|HIGHPASS_OUTPUT, st->mem_hp); st->first = 0; st->count_lost++; @@ -1238,7 +1230,7 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) int pitch_half[2]; int ol_pitch_id=0; #endif - spx_word16_t *out = vout; + spx_word16_t *out = (spx_word16_t*)vout; VARDECL(spx_lsp_t *interp_qlsp); st=(DecState*)state; @@ -1435,6 +1427,7 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) int qe; qe = speex_bits_unpack_unsigned(bits, 5); #ifdef FIXED_POINT + /* FIXME: Perhaps we could slightly lower the gain here when the output is going to saturate? */ ol_gain = MULT16_32_Q15(28406,ol_gain_table[qe]); #else ol_gain = SIG_SCALING*exp(qe/3.5); @@ -1584,7 +1577,7 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) if (SUBMODE(innovation_unquant)) { /*Fixed codebook contribution*/ - SUBMODE(innovation_unquant)(innov, SUBMODE(innovation_params), st->subframeSize, bits, stack); + SUBMODE(innovation_unquant)(innov, SUBMODE(innovation_params), st->subframeSize, bits, stack, &st->seed); } else { speex_error("No fixed codebook"); } @@ -1628,7 +1621,7 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) } } else { for (i=0;isubframeSize;i++) - exc[i]=PSHR32(ADD32(SHL32(exc32[i],1),innov[i]),SIG_SHIFT); + exc[i]=EXTRACT16(SATURATE32(PSHR32(ADD32(SHL32(exc32[i],1),innov[i]),SIG_SHIFT),32767)); /*print_vec(exc, 40, "innov");*/ } if (innov_save) @@ -1644,10 +1637,10 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) ALLOC(innov2, st->subframeSize, spx_sig_t); for (i=0;isubframeSize;i++) innov2[i]=0; - SUBMODE(innovation_unquant)(innov2, SUBMODE(innovation_params), st->subframeSize, bits, stack); + SUBMODE(innovation_unquant)(innov2, SUBMODE(innovation_params), st->subframeSize, bits, stack, &st->seed); signal_mul(innov2, innov2, MULT16_32_Q15(QCONST16(0.454545,15),ener), st->subframeSize); for (i=0;isubframeSize;i++) - exc[i] = ADD16(exc[i],PSHR32(innov2[i],SIG_SHIFT)); + exc[i] = EXTRACT16(SATURATE32(ADD32(EXTEND32(exc[i]),PSHR32(innov2[i],SIG_SHIFT)),32767)); if (innov_save) { for (i=0;isubframeSize;i++) @@ -1678,8 +1671,8 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) exc_ener = compute_rms16 (st->exc, st->frameSize); gain32 = PDIV32(ol_gain, ADD16(exc_ener,1)); #ifdef FIXED_POINT - if (gain32 > 32768) - gain32 = 32768; + if (gain32 > 32767) + gain32 = 32767; gain = EXTRACT16(gain32); #else if (gain32 > 2) @@ -1734,6 +1727,8 @@ int nb_decode(void *state, SpeexBits *bits, void *vout) } + if (st->highpass_enabled) + highpass(out, out, st->frameSize, (st->isWideband?HIGHPASS_WIDEBAND:HIGHPASS_NARROWBAND)|HIGHPASS_OUTPUT, st->mem_hp); /*for (i=0;iframeSize;i++) printf ("%d\n", (int)st->frame[i]);*/ @@ -1766,40 +1761,40 @@ int nb_encoder_ctl(void *state, int request, void *ptr) switch(request) { case SPEEX_GET_FRAME_SIZE: - (*(int*)ptr) = st->frameSize; + (*(spx_int32_t*)ptr) = st->frameSize; break; case SPEEX_SET_LOW_MODE: case SPEEX_SET_MODE: - st->submodeSelect = st->submodeID = (*(int*)ptr); + st->submodeSelect = st->submodeID = (*(spx_int32_t*)ptr); break; case SPEEX_GET_LOW_MODE: case SPEEX_GET_MODE: - (*(int*)ptr) = st->submodeID; + (*(spx_int32_t*)ptr) = st->submodeID; break; case SPEEX_SET_VBR: - st->vbr_enabled = (*(int*)ptr); + st->vbr_enabled = (*(spx_int32_t*)ptr); break; case SPEEX_GET_VBR: - (*(int*)ptr) = st->vbr_enabled; + (*(spx_int32_t*)ptr) = st->vbr_enabled; break; case SPEEX_SET_VAD: - st->vad_enabled = (*(int*)ptr); + st->vad_enabled = (*(spx_int32_t*)ptr); break; case SPEEX_GET_VAD: - (*(int*)ptr) = st->vad_enabled; + (*(spx_int32_t*)ptr) = st->vad_enabled; break; case SPEEX_SET_DTX: - st->dtx_enabled = (*(int*)ptr); + st->dtx_enabled = (*(spx_int32_t*)ptr); break; case SPEEX_GET_DTX: - (*(int*)ptr) = st->dtx_enabled; + (*(spx_int32_t*)ptr) = st->dtx_enabled; break; case SPEEX_SET_ABR: st->abr_enabled = (*(spx_int32_t*)ptr); st->vbr_enabled = st->abr_enabled!=0; if (st->vbr_enabled) { - int i=10; + spx_int32_t i=10; spx_int32_t rate, target; float vbr_qual; target = (*(spx_int32_t*)ptr); @@ -1832,7 +1827,7 @@ int nb_encoder_ctl(void *state, int request, void *ptr) break; case SPEEX_SET_QUALITY: { - int quality = (*(int*)ptr); + int quality = (*(spx_int32_t*)ptr); if (quality < 0) quality = 0; if (quality > 10) @@ -1841,7 +1836,7 @@ int nb_encoder_ctl(void *state, int request, void *ptr) } break; case SPEEX_SET_COMPLEXITY: - st->complexity = (*(int*)ptr); + st->complexity = (*(spx_int32_t*)ptr); if (st->complexity<0) st->complexity=0; break; @@ -1850,7 +1845,7 @@ int nb_encoder_ctl(void *state, int request, void *ptr) break; case SPEEX_SET_BITRATE: { - int i=10; + spx_int32_t i=10; spx_int32_t rate, target; target = (*(spx_int32_t*)ptr); while (i>=0) @@ -1891,21 +1886,21 @@ int nb_encoder_ctl(void *state, int request, void *ptr) } break; case SPEEX_SET_SUBMODE_ENCODING: - st->encode_submode = (*(int*)ptr); + st->encode_submode = (*(spx_int32_t*)ptr); break; case SPEEX_GET_SUBMODE_ENCODING: - (*(int*)ptr) = st->encode_submode; + (*(spx_int32_t*)ptr) = st->encode_submode; break; case SPEEX_GET_LOOKAHEAD: - (*(int*)ptr)=(st->windowSize-st->frameSize); + (*(spx_int32_t*)ptr)=(st->windowSize-st->frameSize); break; case SPEEX_SET_PLC_TUNING: - st->plc_tuning = (*(int*)ptr); + st->plc_tuning = (*(spx_int32_t*)ptr); if (st->plc_tuning>100) st->plc_tuning=100; break; case SPEEX_GET_PLC_TUNING: - (*(int*)ptr)=(st->plc_tuning); + (*(spx_int32_t*)ptr)=(st->plc_tuning); break; case SPEEX_SET_VBR_MAX_BITRATE: st->vbr_max = (*(spx_int32_t*)ptr); @@ -1913,7 +1908,12 @@ int nb_encoder_ctl(void *state, int request, void *ptr) case SPEEX_GET_VBR_MAX_BITRATE: (*(spx_int32_t*)ptr) = st->vbr_max; break; - + case SPEEX_SET_HIGHPASS: + st->highpass_enabled = (*(spx_int32_t*)ptr); + break; + case SPEEX_GET_HIGHPASS: + (*(spx_int32_t*)ptr) = st->highpass_enabled; + break; /* This is all internal stuff past this point */ case SPEEX_GET_PI_GAIN: @@ -1936,7 +1936,10 @@ int nb_encoder_ctl(void *state, int request, void *ptr) (*(float*)ptr)=st->relative_quality; break; case SPEEX_SET_INNOVATION_SAVE: - st->innov_save = ptr; + st->innov_save = (spx_sig_t*)ptr; + break; + case SPEEX_SET_WIDEBAND: + st->isWideband = *((spx_int32_t*)ptr); break; default: speex_warning_int("Unknown nb_ctl request: ", request); @@ -1953,20 +1956,20 @@ int nb_decoder_ctl(void *state, int request, void *ptr) { case SPEEX_SET_LOW_MODE: case SPEEX_SET_MODE: - st->submodeID = (*(int*)ptr); + st->submodeID = (*(spx_int32_t*)ptr); break; case SPEEX_GET_LOW_MODE: case SPEEX_GET_MODE: - (*(int*)ptr) = st->submodeID; + (*(spx_int32_t*)ptr) = st->submodeID; break; case SPEEX_SET_ENH: - st->lpc_enh_enabled = *((int*)ptr); + st->lpc_enh_enabled = *((spx_int32_t*)ptr); break; case SPEEX_GET_ENH: - *((int*)ptr) = st->lpc_enh_enabled; + *((spx_int32_t*)ptr) = st->lpc_enh_enabled; break; case SPEEX_GET_FRAME_SIZE: - (*(int*)ptr) = st->frameSize; + (*(spx_int32_t*)ptr) = st->frameSize; break; case SPEEX_GET_BITRATE: if (st->submodes[st->submodeID]) @@ -2006,14 +2009,21 @@ int nb_decoder_ctl(void *state, int request, void *ptr) } break; case SPEEX_SET_SUBMODE_ENCODING: - st->encode_submode = (*(int*)ptr); + st->encode_submode = (*(spx_int32_t*)ptr); break; case SPEEX_GET_SUBMODE_ENCODING: - (*(int*)ptr) = st->encode_submode; + (*(spx_int32_t*)ptr) = st->encode_submode; break; case SPEEX_GET_LOOKAHEAD: - (*(int*)ptr)=st->subframeSize; + (*(spx_int32_t*)ptr)=st->subframeSize; break; + case SPEEX_SET_HIGHPASS: + st->highpass_enabled = (*(spx_int32_t*)ptr); + break; + case SPEEX_GET_HIGHPASS: + (*(spx_int32_t*)ptr) = st->highpass_enabled; + break; + case SPEEX_GET_PI_GAIN: { int i; @@ -2031,10 +2041,13 @@ int nb_decoder_ctl(void *state, int request, void *ptr) } break; case SPEEX_GET_DTX_STATUS: - *((int*)ptr) = st->dtx_enabled; + *((spx_int32_t*)ptr) = st->dtx_enabled; break; case SPEEX_SET_INNOVATION_SAVE: - st->innov_save = ptr; + st->innov_save = (spx_sig_t*)ptr; + break; + case SPEEX_SET_WIDEBAND: + st->isWideband = *((spx_int32_t*)ptr); break; default: speex_warning_int("Unknown nb_ctl request: ", request); diff --git a/pjmedia/src/pjmedia-codec/speex/nb_celp.h b/pjmedia/src/pjmedia-codec/speex/nb_celp.h index 92028cb5..9a61d681 100644 --- a/pjmedia/src/pjmedia-codec/speex/nb_celp.h +++ b/pjmedia/src/pjmedia-codec/speex/nb_celp.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin */ +/* Copyright (C) 2002-2006 Jean-Marc Valin */ /** @file nb_celp.h @brief Narrowband CELP encoder/decoder @@ -94,13 +94,14 @@ typedef struct EncState { spx_mem_t *mem_sw_whole; /**< Filter memory for perceptually-weighted signal (whole frame)*/ spx_mem_t *mem_exc; /**< Filter memory for excitation (whole frame) */ spx_mem_t *mem_exc2; /**< Filter memory for excitation (whole frame) */ + spx_mem_t mem_hp[2]; /**< High-pass filter memory */ spx_word32_t *pi_gain; /**< Gain of LPC filter at theta=pi (fe/2) */ spx_sig_t *innov_save; /**< If non-NULL, innovation is copied here */ VBRState *vbr; /**< State of the VBR data */ float vbr_quality; /**< Quality setting for VBR encoding */ float relative_quality; /**< Relative quality that will be needed by VBR */ - int vbr_enabled; /**< 1 for enabling VBR, 0 otherwise */ + spx_int32_t vbr_enabled; /**< 1 for enabling VBR, 0 otherwise */ spx_int32_t vbr_max; /**< Max bit-rate allowed in VBR mode */ int vad_enabled; /**< 1 for enabling VAD, 0 otherwise */ int dtx_enabled; /**< 1 for enabling DTX, 0 otherwise */ @@ -116,6 +117,8 @@ typedef struct EncState { const SpeexSubmode * const *submodes; /**< Sub-mode data */ int submodeID; /**< Activated sub-mode */ int submodeSelect; /**< Mode chosen by the user (may differ from submodeID if VAD is on) */ + int isWideband; /**< Is this used as part of the embedded wideband codec */ + int highpass_enabled; /**< Is the input filter enabled */ } EncState; /**Structure representing the full state of the narrowband decoder*/ @@ -143,6 +146,7 @@ typedef struct DecState { spx_lsp_t *old_qlsp; /**< Quantized LSPs for previous frame */ spx_coef_t *interp_qlpc; /**< Interpolated quantized LPCs */ spx_mem_t *mem_sp; /**< Filter memory for synthesis signal */ + spx_mem_t mem_hp[2]; /**< High-pass filter memory */ spx_word32_t *pi_gain; /**< Gain of LPC filter at theta=pi (fe/2) */ spx_sig_t *innov_save; /** If non-NULL, innovation is copied here */ @@ -168,6 +172,8 @@ typedef struct DecState { int voc_offset; int dtx_enabled; + int isWideband; /**< Is this used as part of the embedded wideband codec */ + int highpass_enabled; /**< Is the input filter enabled */ } DecState; /** Initializes encoder state*/ diff --git a/pjmedia/src/pjmedia-codec/speex/preprocess_spx.c b/pjmedia/src/pjmedia-codec/speex/preprocess_spx.c index ba2e6076..d82230cf 100644 --- a/pjmedia/src/pjmedia-codec/speex/preprocess_spx.c +++ b/pjmedia/src/pjmedia-codec/speex/preprocess_spx.c @@ -1,6 +1,6 @@ -/* Copyright (C) 2003 Epic Games - Written by Jean-Marc Valin - +/* Copyright (C) 2003 Epic Games (written by Jean-Marc Valin) + Copyright (C) 2004-2006 Epic Games + File: preprocess.c Preprocessor with denoising based on the algorithm by Ephraim and Malah @@ -31,96 +31,356 @@ POSSIBILITY OF SUCH DAMAGE. */ + +/* + Recommended papers: + + Y. Ephraim and D. Malah, "Speech enhancement using minimum mean-square error + short-time spectral amplitude estimator". IEEE Transactions on Acoustics, + Speech and Signal Processing, vol. ASSP-32, no. 6, pp. 1109-1121, 1984. + + Y. Ephraim and D. Malah, "Speech enhancement using minimum mean-square error + log-spectral amplitude estimator". IEEE Transactions on Acoustics, Speech and + Signal Processing, vol. ASSP-33, no. 2, pp. 443-445, 1985. + + I. Cohen and B. Berdugo, "Speech enhancement for non-stationary noise environments". + Signal Processing, vol. 81, no. 2, pp. 2403-2418, 2001. + + Stefan Gustafsson, Rainer Martin, Peter Jax, and Peter Vary. "A psychoacoustic + approach to combined acoustic echo cancellation and noise reduction". IEEE + Transactions on Speech and Audio Processing, 2002. + + J.-M. Valin, J. Rouat, and F. Michaud, "Microphone array post-filter for separation + of simultaneous non-stationary sources". In Proceedings IEEE International + Conference on Acoustics, Speech, and Signal Processing, 2004. +*/ + #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include "speex/speex_preprocess.h" +#include "speex/speex_echo.h" #include "misc.h" -#include "smallft.h" - -#define max(a,b) ((a) > (b) ? (a) : (b)) -#define min(a,b) ((a) < (b) ? (a) : (b)) +#include "fftwrap.h" +#include "filterbank.h" +#include "math_approx.h" #ifndef M_PI #define M_PI 3.14159263 #endif -#define SQRT_M_PI_2 0.88623 #define LOUDNESS_EXP 2.5 -#define NB_BANDS 8 +#define NB_BANDS 24 + +#define SPEECH_PROB_START_DEFAULT QCONST16(0.35f,15) +#define SPEECH_PROB_CONTINUE_DEFAULT QCONST16(0.20f,15) +#define NOISE_SUPPRESS_DEFAULT -25 +#define ECHO_SUPPRESS_DEFAULT -45 +#define ECHO_SUPPRESS_ACTIVE_DEFAULT -15 + +#ifndef NULL +#define NULL 0 +#endif + +#define SQR(x) ((x)*(x)) +#define SQR16(x) (MULT16_16((x),(x))) +#define SQR16_Q15(x) (MULT16_16_Q15((x),(x))) + +#ifdef FIXED_POINT +static inline spx_word16_t DIV32_16_Q8(spx_word32_t a, spx_word32_t b) +{ + if (SHR32(a,7) >= b) + { + return 32767; + } else { + if (b>=QCONST32(1,23)) + { + a = SHR32(a,8); + b = SHR32(b,8); + } + if (b>=QCONST32(1,19)) + { + a = SHR32(a,4); + b = SHR32(b,4); + } + if (b>=QCONST32(1,15)) + { + a = SHR32(a,4); + b = SHR32(b,4); + } + a = SHL32(a,8); + return PDIV32_16(a,b); + } + +} +static inline spx_word16_t DIV32_16_Q15(spx_word32_t a, spx_word32_t b) +{ + if (SHR32(a,15) >= b) + { + return 32767; + } else { + if (b>=QCONST32(1,23)) + { + a = SHR32(a,8); + b = SHR32(b,8); + } + if (b>=QCONST32(1,19)) + { + a = SHR32(a,4); + b = SHR32(b,4); + } + if (b>=QCONST32(1,15)) + { + a = SHR32(a,4); + b = SHR32(b,4); + } + a = SHL32(a,15)-a; + return DIV32_16(a,b); + } +} +#define SNR_SCALING 256.f +#define SNR_SCALING_1 0.0039062f +#define SNR_SHIFT 8 + +#define FRAC_SCALING 32767.f +#define FRAC_SCALING_1 3.0518e-05 +#define FRAC_SHIFT 1 -#define SPEEX_PROB_START_DEFAULT 0.35f -#define SPEEX_PROB_CONTINUE_DEFAULT 0.20f +#define EXPIN_SCALING 2048.f +#define EXPIN_SCALING_1 0.00048828f +#define EXPIN_SHIFT 11 +#define EXPOUT_SCALING_1 1.5259e-05 -#define ZMIN .1 -#define ZMAX .316 -#define ZMIN_1 10 -#define LOG_MIN_MAX_1 0.86859 +#define NOISE_SHIFT 7 -static void conj_window(float *w, int len) +#else + +#define DIV32_16_Q8(a,b) ((a)/(b)) +#define DIV32_16_Q15(a,b) ((a)/(b)) +#define SNR_SCALING 1.f +#define SNR_SCALING_1 1.f +#define SNR_SHIFT 0 +#define FRAC_SCALING 1.f +#define FRAC_SCALING_1 1.f +#define FRAC_SHIFT 0 +#define NOISE_SHIFT 0 + +#define EXPIN_SCALING 1.f +#define EXPIN_SCALING_1 1.f +#define EXPOUT_SCALING_1 1.f + +#endif + +/** Speex pre-processor state. */ +struct SpeexPreprocessState_ { + /* Basic info */ + int frame_size; /**< Number of samples processed each time */ + int ps_size; /**< Number of points in the power spectrum */ + int sampling_rate; /**< Sampling rate of the input/output */ + int nbands; + FilterBank *bank; + + /* Parameters */ + int denoise_enabled; + int agc_enabled; + float agc_level; + int vad_enabled; + int dereverb_enabled; + spx_word16_t reverb_decay; + spx_word16_t reverb_level; + spx_word16_t speech_prob_start; + spx_word16_t speech_prob_continue; + int noise_suppress; + int echo_suppress; + int echo_suppress_active; + SpeexEchoState *echo_state; + + /* DSP-related arrays */ + spx_word16_t *frame; /**< Processing frame (2*ps_size) */ + spx_word16_t *ft; /**< Processing frame in freq domain (2*ps_size) */ + spx_word32_t *ps; /**< Current power spectrum */ + spx_word16_t *gain2; /**< Adjusted gains */ + spx_word16_t *gain_floor; /**< Minimum gain allowed */ + spx_word16_t *window; /**< Analysis/Synthesis window */ + spx_word32_t *noise; /**< Noise estimate */ + spx_word32_t *reverb_estimate; /**< Estimate of reverb energy */ + spx_word32_t *old_ps; /**< Power spectrum for last frame */ + spx_word16_t *gain; /**< Ephraim Malah gain */ + spx_word16_t *prior; /**< A-priori SNR */ + spx_word16_t *post; /**< A-posteriori SNR */ + + spx_word32_t *S; /**< Smoothed power spectrum */ + spx_word32_t *Smin; /**< See Cohen paper */ + spx_word32_t *Stmp; /**< See Cohen paper */ + int *update_prob; /**< Propability of speech presence for noise update */ + + spx_word16_t *zeta; /**< Smoothed a priori SNR */ + spx_word32_t *echo_noise; + spx_word32_t *residual_echo; + + /* Misc */ + spx_word16_t *inbuf; /**< Input buffer (overlapped analysis) */ + spx_word16_t *outbuf; /**< Output buffer (for overlap and add) */ + +#ifndef FIXED_POINT + float *loudness_weight; /**< Perceptual loudness curve */ + float loudness; /**< loudness estimate */ + float loudness2; /**< loudness estimate */ + int nb_loudness_adapt; /**< Number of frames used for loudness adaptation so far */ +#endif + int nb_adapt; /**< Number of frames used for adaptation so far */ + int was_speech; + int min_count; /**< Number of frames processed so far */ + void *fft_lookup; /**< Lookup table for the FFT */ +#ifdef FIXED_POINT + int frame_shift; +#endif +}; + + +static void conj_window(spx_word16_t *w, int len) { int i; for (i=0;i19) + return ADD32(EXTEND32(Q15_ONE),EXTEND32(DIV32_16(QCONST32(.1296,23), SHR32(xx,EXPIN_SHIFT-SNR_SHIFT)))); + frac = SHL32(xx-SHL32(ind,10),5); + return SHL32(DIV32_16(PSHR32(MULT16_16(Q15_ONE-frac,table[ind]) + MULT16_16(frac,table[ind+1]),7),(spx_sqrt(SHL32(xx,15)+6711))),7); +} + +static inline spx_word16_t qcurve(spx_word16_t x) +{ + x = MAX16(x, 1); + return DIV32_16(SHL32(EXTEND32(32767),9),ADD16(512,MULT16_16_Q15(QCONST16(.60f,15),DIV32_16(32767,x)))); +} + +/* Compute the gain floor based on different floors for the background noise and residual echo */ +static void compute_gain_floor(int noise_suppress, int effective_echo_suppress, spx_word32_t *noise, spx_word32_t *echo, spx_word16_t *gain_floor, int len) +{ + int i; + + if (noise_suppress > effective_echo_suppress) + { + spx_word16_t noise_gain, gain_ratio; + noise_gain = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(0.11513,11),noise_suppress)),1))); + gain_ratio = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(.2302585f,11),effective_echo_suppress-noise_suppress)),1))); + + /* gain_floor = sqrt [ (noise*noise_floor + echo*echo_floor) / (noise+echo) ] */ + for (i=0;i9.5) - return 1+.1296/x; - - integer = floor(2*x); - frac = 2*x-integer; - ind = (int)integer; - - return ((1-frac)*table[ind] + frac*table[ind+1])/sqrt(x+.0001f); + x = EXPIN_SCALING_1*xx; + integer = floor(2*x); + ind = (int)integer; + if (ind<0) + return FRAC_SCALING; + if (ind>19) + return FRAC_SCALING*(1+.1296/x); + frac = 2*x-integer; + return FRAC_SCALING*((1-frac)*table[ind] + frac*table[ind+1])/sqrt(x+.0001f); } -static inline float qcurve(float x) +static inline spx_word16_t qcurve(spx_word16_t x) { - return 1.f/(1.f+.1f/(x*x)); + return 1.f/(1.f+.15f/(SNR_SCALING_1*x)); } +static void compute_gain_floor(int noise_suppress, int effective_echo_suppress, spx_word32_t *noise, spx_word32_t *echo, spx_word16_t *gain_floor, int len) +{ + int i; + float echo_floor; + float noise_floor; + + noise_floor = exp(.2302585f*noise_suppress); + echo_floor = exp(.2302585f*effective_echo_suppress); + + /* Compute the gain floor based on different floors for the background noise and residual echo */ + for (i=0;iframe_size = frame_size; @@ -157,45 +417,49 @@ SpeexPreprocessState *speex_preprocess_state_init(int frame_size, int sampling_r st->agc_level = 8000; st->vad_enabled = 0; st->dereverb_enabled = 0; - st->reverb_decay = .5; - st->reverb_level = .2; - - st->speech_prob_start = SPEEX_PROB_START_DEFAULT; - st->speech_prob_continue = SPEEX_PROB_CONTINUE_DEFAULT; - - st->frame = (float*)speex_alloc(2*N*sizeof(float)); - st->ps = (float*)speex_alloc(N*sizeof(float)); - st->gain2 = (float*)speex_alloc(N*sizeof(float)); - st->window = (float*)speex_alloc(2*N*sizeof(float)); - st->noise = (float*)speex_alloc(N*sizeof(float)); - st->reverb_estimate = (float*)speex_alloc(N*sizeof(float)); - st->old_ps = (float*)speex_alloc(N*sizeof(float)); - st->gain = (float*)speex_alloc(N*sizeof(float)); - st->prior = (float*)speex_alloc(N*sizeof(float)); - st->post = (float*)speex_alloc(N*sizeof(float)); - st->loudness_weight = (float*)speex_alloc(N*sizeof(float)); - st->inbuf = (float*)speex_alloc(N3*sizeof(float)); - st->outbuf = (float*)speex_alloc(N3*sizeof(float)); - st->echo_noise = (float*)speex_alloc(N*sizeof(float)); - - st->S = (float*)speex_alloc(N*sizeof(float)); - st->Smin = (float*)speex_alloc(N*sizeof(float)); - st->Stmp = (float*)speex_alloc(N*sizeof(float)); - st->update_prob = (float*)speex_alloc(N*sizeof(float)); + st->reverb_decay = 0; + st->reverb_level = 0; + st->noise_suppress = NOISE_SUPPRESS_DEFAULT; + st->echo_suppress = ECHO_SUPPRESS_DEFAULT; + st->echo_suppress_active = ECHO_SUPPRESS_ACTIVE_DEFAULT; - st->zeta = (float*)speex_alloc(N*sizeof(float)); - st->Zpeak = 0; - st->Zlast = 0; + st->speech_prob_start = SPEECH_PROB_START_DEFAULT; + st->speech_prob_continue = SPEECH_PROB_CONTINUE_DEFAULT; - st->noise_bands = (float*)speex_alloc(NB_BANDS*sizeof(float)); - st->noise_bands2 = (float*)speex_alloc(NB_BANDS*sizeof(float)); - st->speech_bands = (float*)speex_alloc(NB_BANDS*sizeof(float)); - st->speech_bands2 = (float*)speex_alloc(NB_BANDS*sizeof(float)); - st->noise_bandsN = st->speech_bandsN = 1; + st->echo_state = NULL; + + st->nbands = NB_BANDS; + M = st->nbands; + st->bank = filterbank_new(M, sampling_rate, N, 1); + + st->frame = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t)); + st->window = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t)); + st->ft = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t)); + + st->ps = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t)); + st->noise = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t)); + st->echo_noise = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t)); + st->residual_echo = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t)); + st->reverb_estimate = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t)); + st->old_ps = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t)); + st->prior = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t)); + st->post = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t)); + st->gain = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t)); + st->gain2 = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t)); + st->gain_floor = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t)); + st->zeta = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t)); + + st->S = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); + st->Smin = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); + st->Stmp = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t)); + st->update_prob = (int*)speex_alloc(N*sizeof(int)); + + st->inbuf = (spx_word16_t*)speex_alloc(N3*sizeof(spx_word16_t)); + st->outbuf = (spx_word16_t*)speex_alloc(N3*sizeof(spx_word16_t)); conj_window(st->window, 2*N3); for (i=2*N3;i<2*st->ps_size;i++) - st->window[i]=1; + st->window[i]=Q15_ONE; if (N4>0) { @@ -205,22 +469,25 @@ SpeexPreprocessState *speex_preprocess_state_init(int frame_size, int sampling_r st->window[i+N3]=1; } } - for (i=0;inoise[i]=1e4; - st->reverb_estimate[i]=0.; - st->old_ps[i]=1e4; - st->gain[i]=1; - st->post[i]=1; - st->prior[i]=1; + st->noise[i]=QCONST32(1.f,NOISE_SHIFT); + st->reverb_estimate[i]=0; + st->old_ps[i]=1; + st->gain[i]=Q15_ONE; + st->post[i]=SHL16(1, SNR_SHIFT); + st->prior[i]=SHL16(1, SNR_SHIFT); } + for (i=0;iupdate_prob[i] = 1; for (i=0;iinbuf[i]=0; st->outbuf[i]=0; } - +#ifndef FIXED_POINT + st->loudness_weight = (float*)speex_alloc(N*sizeof(float)); for (i=0;iloudness_weight[i]=.01f; st->loudness_weight[i] *= st->loudness_weight[i]; } - - st->speech_prob = 0; - st->last_speech = 1000; st->loudness = pow(6000,LOUDNESS_EXP); st->loudness2 = 6000; st->nb_loudness_adapt = 0; +#endif + st->was_speech = 0; - st->fft_lookup = (struct drft_lookup*)speex_alloc(sizeof(struct drft_lookup)); - spx_drft_init(st->fft_lookup,2*N); + st->fft_lookup = spx_fft_init(2*N); st->nb_adapt=0; - st->consec_noise=0; - st->nb_preprocess=0; + st->min_count=0; return st; } void speex_preprocess_state_destroy(SpeexPreprocessState *st) { speex_free(st->frame); + speex_free(st->ft); speex_free(st->ps); speex_free(st->gain2); + speex_free(st->gain_floor); speex_free(st->window); speex_free(st->noise); speex_free(st->reverb_estimate); @@ -257,8 +523,11 @@ void speex_preprocess_state_destroy(SpeexPreprocessState *st) speex_free(st->gain); speex_free(st->prior); speex_free(st->post); +#ifndef FIXED_POINT speex_free(st->loudness_weight); +#endif speex_free(st->echo_noise); + speex_free(st->residual_echo); speex_free(st->S); speex_free(st->Smin); @@ -266,241 +535,17 @@ void speex_preprocess_state_destroy(SpeexPreprocessState *st) speex_free(st->update_prob); speex_free(st->zeta); - speex_free(st->noise_bands); - speex_free(st->noise_bands2); - speex_free(st->speech_bands); - speex_free(st->speech_bands2); - speex_free(st->inbuf); speex_free(st->outbuf); - spx_drft_clear(st->fft_lookup); - speex_free(st->fft_lookup); - + spx_fft_destroy(st->fft_lookup); + filterbank_destroy(st->bank); speex_free(st); } -static void update_noise(SpeexPreprocessState *st, float *ps, spx_int32_t *echo) -{ - int i; - float beta; - st->nb_adapt++; - beta=1.0f/st->nb_adapt; - if (beta < .05f) - beta=.05f; - - if (!echo) - { - for (i=0;ips_size;i++) - st->noise[i] = (1.f-beta)*st->noise[i] + beta*ps[i]; - } else { - for (i=0;ips_size;i++) - st->noise[i] = (1.f-beta)*st->noise[i] + beta*max(1.f,ps[i]-st->frame_size*st->frame_size*4.0*echo[i]); -#if 0 - for (i=0;ips_size;i++) - st->noise[i] = 0; -#endif - } -} - -static int speex_compute_vad(SpeexPreprocessState *st, float *ps, float mean_prior, float mean_post) -{ - int i, is_speech=0; - int N = st->ps_size; - float scale=.5f/N; - - /* FIXME: Clean this up a bit */ - { - float bands[NB_BANDS]; - int j; - float p0, p1; - float tot_loudness=0; - float x = sqrt(mean_post); - - for (i=5;ips[i] * st->loudness_weight[i]; - } - - for (i=0;ispeech_prob + .01*(1-st->speech_prob); - p1 *= .01*st->speech_prob + .99*(1-st->speech_prob); - - st->speech_prob = p0/(p1+p0); - */ - - if (st->noise_bandsN < 50 || st->speech_bandsN < 50) - { - if (mean_post > 5.f) - { - float adapt = 1./st->speech_bandsN++; - if (adapt<.005f) - adapt = .005f; - for (i=0;ispeech_bands[i] = (1.f-adapt)*st->speech_bands[i] + adapt*bands[i]; - /*st->speech_bands2[i] = (1-adapt)*st->speech_bands2[i] + adapt*bands[i]*bands[i];*/ - st->speech_bands2[i] = (1.f-adapt)*st->speech_bands2[i] + adapt*(bands[i]-st->speech_bands[i])*(bands[i]-st->speech_bands[i]); - } - } else { - float adapt = 1./st->noise_bandsN++; - if (adapt<.005f) - adapt = .005f; - for (i=0;inoise_bands[i] = (1.f-adapt)*st->noise_bands[i] + adapt*bands[i]; - /*st->noise_bands2[i] = (1-adapt)*st->noise_bands2[i] + adapt*bands[i]*bands[i];*/ - st->noise_bands2[i] = (1.f-adapt)*st->noise_bands2[i] + adapt*(bands[i]-st->noise_bands[i])*(bands[i]-st->noise_bands[i]); - } - } - } - p0=p1=1; - for (i=0;inoise_bands2[i] - st->noise_bands[i]*st->noise_bands[i]; - speech_var = 1.01*st->speech_bands2[i] - st->speech_bands[i]*st->speech_bands[i];*/ - noise_var = st->noise_bands2[i]; - speech_var = st->speech_bands2[i]; - if (noise_var < .1f) - noise_var = .1f; - if (speech_var < .1f) - speech_var = .1f; - - /*speech_var = sqrt(speech_var*noise_var); - noise_var = speech_var;*/ - if (speech_var < .05f*speech_var) - noise_var = .05f*speech_var; - if (speech_var < .05f*noise_var) - speech_var = .05f*noise_var; - - if (bands[i] < st->noise_bands[i]) - speech_var = noise_var; - if (bands[i] > st->speech_bands[i]) - noise_var = speech_var; - - speech_mean = st->speech_bands[i]; - noise_mean = st->noise_bands[i]; - if (noise_mean < speech_mean - 5.f) - noise_mean = speech_mean - 5.f; - - tmp1 = exp(-.5f*(bands[i]-speech_mean)*(bands[i]-speech_mean)/speech_var)/sqrt(2.f*M_PI*speech_var); - tmp2 = exp(-.5f*(bands[i]-noise_mean)*(bands[i]-noise_mean)/noise_var)/sqrt(2.f*M_PI*noise_var); - /*fprintf (stderr, "%f ", (float)(p0/(.01+p0+p1)));*/ - /*fprintf (stderr, "%f ", (float)(bands[i]));*/ - pr = tmp1/(1e-25+tmp1+tmp2); - /*if (bands[i] < st->noise_bands[i]) - pr=.01; - if (bands[i] > st->speech_bands[i] && pr < .995) - pr=.995;*/ - if (pr>.999f) - pr=.999f; - if (pr<.001f) - pr=.001f; - /*fprintf (stderr, "%f ", pr);*/ - p0 *= pr; - p1 *= (1-pr); - } - - p0 = pow(p0,.2); - p1 = pow(p1,.2); - -#if 1 - p0 *= 2.f; - p0=p0/(p1+p0); - if (st->last_speech>20) - { - float tmp = sqrt(tot_loudness)/st->loudness2; - tmp = 1.f-exp(-10.f*tmp); - if (p0>tmp) - p0=tmp; - } - p1=1-p0; -#else - if (sqrt(tot_loudness) < .6f*st->loudness2 && p0>15.f*p1) - p0=15.f*p1; - if (sqrt(tot_loudness) < .45f*st->loudness2 && p0>7.f*p1) - p0=7.f*p1; - if (sqrt(tot_loudness) < .3f*st->loudness2 && p0>3.f*p1) - p0=3.f*p1; - if (sqrt(tot_loudness) < .15f*st->loudness2 && p0>p1) - p0=p1; - /*fprintf (stderr, "%f %f ", (float)(sqrt(tot_loudness) /( .25*st->loudness2)), p0/(p1+p0));*/ -#endif - - p0 *= .99f*st->speech_prob + .01f*(1-st->speech_prob); - p1 *= .01f*st->speech_prob + .99f*(1-st->speech_prob); - - st->speech_prob = p0/(1e-25f+p1+p0); - /*fprintf (stderr, "%f %f %f ", tot_loudness, st->loudness2, st->speech_prob);*/ - - if (st->speech_prob > st->speech_prob_start - || (st->last_speech < 20 && st->speech_prob > st->speech_prob_continue)) - { - is_speech = 1; - st->last_speech = 0; - } else { - st->last_speech++; - if (st->last_speech<20) - is_speech = 1; - } - - if (st->noise_bandsN > 50 && st->speech_bandsN > 50) - { - if (mean_post > 5) - { - float adapt = 1./st->speech_bandsN++; - if (adapt<.005f) - adapt = .005f; - for (i=0;ispeech_bands[i] = (1-adapt)*st->speech_bands[i] + adapt*bands[i]; - /*st->speech_bands2[i] = (1-adapt)*st->speech_bands2[i] + adapt*bands[i]*bands[i];*/ - st->speech_bands2[i] = (1-adapt)*st->speech_bands2[i] + adapt*(bands[i]-st->speech_bands[i])*(bands[i]-st->speech_bands[i]); - } - } else { - float adapt = 1./st->noise_bandsN++; - if (adapt<.005f) - adapt = .005f; - for (i=0;inoise_bands[i] = (1-adapt)*st->noise_bands[i] + adapt*bands[i]; - /*st->noise_bands2[i] = (1-adapt)*st->noise_bands2[i] + adapt*bands[i]*bands[i];*/ - st->noise_bands2[i] = (1-adapt)*st->noise_bands2[i] + adapt*(bands[i]-st->noise_bands[i])*(bands[i]-st->noise_bands[i]); - } - } - } - - - } - - return is_speech; -} - -static void speex_compute_agc(SpeexPreprocessState *st, float mean_prior) +/* FIXME: The AGC doesn't work yet with fixed-point*/ +#ifndef FIXED_POINT +static void speex_compute_agc(SpeexPreprocessState *st) { int i; int N = st->ps_size; @@ -535,7 +580,7 @@ static void speex_compute_agc(SpeexPreprocessState *st, float mean_prior) for (i=2;ips[i] * st->gain2[i] * st->gain2[i] * st->loudness_weight[i]; + loudness += scale*st->ps[i] * FRAC_SCALING_1*FRAC_SCALING_1*st->gain2[i] * st->gain2[i] * st->loudness_weight[i]; } loudness=sqrt(loudness); /*if (loudness < 2*pow(st->loudness, 1.0/LOUDNESS_EXP) && @@ -558,6 +603,7 @@ static void speex_compute_agc(SpeexPreprocessState *st, float mean_prior) st->gain2[i] *= agc_gain; } +#endif static void preprocess_analysis(SpeexPreprocessState *st, spx_int16_t *x) { @@ -565,7 +611,7 @@ static void preprocess_analysis(SpeexPreprocessState *st, spx_int16_t *x) int N = st->ps_size; int N3 = 2*N - st->frame_size; int N4 = st->frame_size - N3; - float *ps=st->ps; + spx_word32_t *ps=st->ps; /* 'Build' input frame */ for (i=0;iframe[i] *= st->window[i]; + st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]); +#ifdef FIXED_POINT + { + spx_word16_t max_val=0; + for (i=0;i<2*N;i++) + max_val = MAX16(max_val, ABS16(st->frame[i])); + st->frame_shift = 14-spx_ilog2(EXTEND32(max_val)); + for (i=0;i<2*N;i++) + st->frame[i] = SHL16(st->frame[i], st->frame_shift); + } +#endif + /* Perform FFT */ - spx_drft_forward(st->fft_lookup, st->frame); - + spx_fft(st->fft_lookup, st->frame, st->ft); + /* Power spectrum */ - ps[0]=1; + ps[0]=MULT16_16(st->ft[0],st->ft[0]); for (i=1;iframe[2*i-1]*st->frame[2*i-1] + st->frame[2*i]*st->frame[2*i]; + ps[i]=MULT16_16(st->ft[2*i-1],st->ft[2*i-1]) + MULT16_16(st->ft[2*i],st->ft[2*i]); + for (i=0;ips[i] = PSHR32(st->ps[i], 2*st->frame_shift); + filterbank_compute_bank32(st->bank, ps, ps+N); } static void update_noise_prob(SpeexPreprocessState *st) { int i; + int min_range; int N = st->ps_size; for (i=1;iS[i] = 100.f+ .8f*st->S[i] + .05f*st->ps[i-1]+.1f*st->ps[i]+.05f*st->ps[i+1]; + st->S[i] = MULT16_32_Q15(QCONST16(.8f,15),st->S[i]) + MULT16_32_Q15(QCONST16(.05f,15),st->ps[i-1]) + + MULT16_32_Q15(QCONST16(.1f,15),st->ps[i]) + MULT16_32_Q15(QCONST16(.05f,15),st->ps[i+1]); + st->S[0] = MULT16_32_Q15(QCONST16(.8f,15),st->S[0]) + MULT16_32_Q15(QCONST16(.2f,15),st->ps[0]); + st->S[N-1] = MULT16_32_Q15(QCONST16(.8f,15),st->S[N-1]) + MULT16_32_Q15(QCONST16(.2f,15),st->ps[N-1]); - if (st->nb_preprocess<1) + if (st->nb_adapt==1) { - for (i=1;iSmin[i] = st->Stmp[i] = st->S[i]+100.f; + for (i=0;iSmin[i] = st->Stmp[i] = 0; } - if (st->nb_preprocess%200==0) + if (st->nb_adapt < 100) + min_range = 15; + else if (st->nb_adapt < 1000) + min_range = 50; + else if (st->nb_adapt < 10000) + min_range = 150; + else + min_range = 300; + if (st->min_count > min_range) { - for (i=1;imin_count = 0; + for (i=0;iSmin[i] = min(st->Stmp[i], st->S[i]); + st->Smin[i] = MIN32(st->Stmp[i], st->S[i]); st->Stmp[i] = st->S[i]; } } else { - for (i=1;iSmin[i] = min(st->Smin[i], st->S[i]); - st->Stmp[i] = min(st->Stmp[i], st->S[i]); + st->Smin[i] = MIN32(st->Smin[i], st->S[i]); + st->Stmp[i] = MIN32(st->Stmp[i], st->S[i]); } } - for (i=1;iupdate_prob[i] *= .2f; - if (st->S[i] > 2.5*st->Smin[i]) - st->update_prob[i] += .8f; + if (MULT16_32_Q15(QCONST16(.4f,15),st->S[i]) > ADD32(st->Smin[i],EXTEND32(20))) + st->update_prob[i] = 1; + else + st->update_prob[i] = 0; /*fprintf (stderr, "%f ", st->S[i]/st->Smin[i]);*/ /*fprintf (stderr, "%f ", st->update_prob[i]);*/ } } -#define NOISE_OVERCOMPENS 1.4 +#define NOISE_OVERCOMPENS 1. + +void speex_echo_get_residual(SpeexEchoState *st, spx_word32_t *Yout, int len); int speex_preprocess(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo) +{ + return speex_preprocess_run(st, x); +} + +int speex_preprocess_run(SpeexPreprocessState *st, spx_int16_t *x) { int i; - int is_speech=1; - float mean_post=0; - float mean_prior=0; + int M; int N = st->ps_size; int N3 = 2*N - st->frame_size; int N4 = st->frame_size - N3; - float scale=.5f/N; - float *ps=st->ps; - float Zframe=0, Pframe; - + spx_word32_t *ps=st->ps; + spx_word32_t Zframe; + spx_word16_t Pframe; + spx_word16_t beta, beta_1; + spx_word16_t effective_echo_suppress; + + st->nb_adapt++; + st->min_count++; + + beta = MAX16(QCONST16(.03,15),DIV32_16(Q15_ONE,st->nb_adapt)); + beta_1 = Q15_ONE-beta; + M = st->nbands; + /* Deal with residual echo if provided */ + if (st->echo_state) + { + speex_echo_get_residual(st->echo_state, st->residual_echo, N); +#ifndef FIXED_POINT + /* If there are NaNs or ridiculous values, it'll show up in the DC and we just reset everything to zero */ + if (!(st->residual_echo[0] >=0 && st->residual_echo[0]residual_echo[i] = 0; + } +#endif + for (i=0;iecho_noise[i] = MAX32(MULT16_32_Q15(QCONST16(.6f,15),st->echo_noise[i]), st->residual_echo[i]); + filterbank_compute_bank32(st->bank, st->echo_noise, st->echo_noise+N); + } else { + for (i=0;iecho_noise[i] = 0; + } preprocess_analysis(st, x); update_noise_prob(st); - st->nb_preprocess++; - - /* Noise estimation always updated for the 20 first times */ - if (st->nb_adapt<10) + /* Noise estimation always updated for the 10 first frames */ + /*if (st->nb_adapt<10) { - update_noise(st, ps, echo); + for (i=1;iupdate_prob[i] = 0; } - - /* Deal with residual echo if provided */ - if (echo) - for (i=1;iecho_noise[i] = (.3f*st->echo_noise[i] + st->frame_size*st->frame_size*4.0*echo[i]); - - /* Compute a posteriori SNR */ - for (i=1;inoise[i] + st->echo_noise[i] + st->reverb_estimate[i]; - st->post[i] = ps[i]/tot_noise - 1.f; - if (st->post[i]>100.f) - st->post[i]=100.f; - /*if (st->post[i]<0) - st->post[i]=0;*/ - mean_post+=st->post[i]; + if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i], NOISE_SHIFT)) + st->noise[i] = MAX32(EXTEND32(0),MULT16_32_Q15(beta_1,st->noise[i]) + MULT16_32_Q15(beta,SHL32(st->ps[i],NOISE_SHIFT))); } - mean_post /= N; - if (mean_post<0.f) - mean_post=0.f; + filterbank_compute_bank32(st->bank, st->noise, st->noise+N); /* Special case for first frame */ if (st->nb_adapt==1) - for (i=1;iold_ps[i] = ps[i]; - /* Compute a priori SNR */ - { - /* A priori update rate */ - for (i=1;iprior[i]*st->prior[i]/((1+st->prior[i])*(1+st->prior[i])); - float tot_noise = 1.f+ NOISE_OVERCOMPENS*st->noise[i] + st->echo_noise[i] + st->reverb_estimate[i]; - /* A priori SNR update */ - st->prior[i] = gamma*max(0.0f,st->post[i]) + - (1.f-gamma)* (.8*st->gain[i]*st->gain[i]*st->old_ps[i]/tot_noise + .2*st->prior[i]); - - if (st->prior[i]>100.f) - st->prior[i]=100.f; - - mean_prior+=st->prior[i]; - } - } - mean_prior /= N; - -#if 0 - for (i=0;iprior[i]); - } - fprintf (stderr, "\n"); -#endif - /*fprintf (stderr, "%f %f\n", mean_prior,mean_post);*/ - - if (st->nb_preprocess>=20) + /* Compute a posteriori SNR */ + for (i=0;inoise[i]; - sig_ener += ps[i]; - } - if (noise_ener > 3.f*sig_ener) - do_update = 1; - /*do_update = 0;*/ - if (do_update) - { - st->consec_noise++; - } else { - st->consec_noise=0; - } + spx_word16_t gamma; + + /* Total noise estimate including residual echo and reverberation */ + spx_word32_t tot_noise = ADD32(ADD32(ADD32(EXTEND32(1), PSHR32(st->noise[i],NOISE_SHIFT)) , st->echo_noise[i]) , st->reverb_estimate[i]); + + /* A posteriori SNR = ps/noise - 1*/ + st->post[i] = SUB16(DIV32_16_Q8(ps[i],tot_noise), QCONST16(1.f,SNR_SHIFT)); + st->post[i]=MIN16(st->post[i], QCONST16(100.f,SNR_SHIFT)); + + /* Computing update gamma = .1 + .9*(old/(old+noise))^2 */ + gamma = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.89f,15),SQR16_Q15(DIV32_16_Q15(st->old_ps[i],ADD32(st->old_ps[i],tot_noise)))); + + /* A priori SNR update = gamma*max(0,post) + (1-gamma)*old/noise */ + st->prior[i] = EXTRACT16(PSHR32(ADD32(MULT16_16(gamma,MAX16(0,st->post[i])), MULT16_16(Q15_ONE-gamma,DIV32_16_Q8(st->old_ps[i],tot_noise))), 15)); + st->prior[i]=MIN16(st->prior[i], QCONST16(100.f,SNR_SHIFT)); } - if (st->vad_enabled) - is_speech = speex_compute_vad(st, ps, mean_prior, mean_post); - + /*print_vec(st->post, N+M, "");*/ - if (st->consec_noise>=3) + /* Recursive average of the a priori SNR. A bit smoothed for the psd components */ + st->zeta[0] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[0]), MULT16_16(QCONST16(.3f,15),st->prior[0])),15); + for (i=1;izeta[i] = PSHR32(ADD32(ADD32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.15f,15),st->prior[i])), + MULT16_16(QCONST16(.075f,15),st->prior[i-1])), MULT16_16(QCONST16(.075f,15),st->prior[i+1])),15); + for (i=N-1;izeta[i] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.3f,15),st->prior[i])),15); + + /* Speech probability of presence for the entire frame is based on the average filterbank a priori SNR */ + Zframe = 0; + for (i=N;izeta[i])); + Pframe = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.899f,15),qcurve(DIV32_16(Zframe,st->nbands))); + + effective_echo_suppress = EXTRACT16(PSHR32(ADD32(MULT16_16(SUB16(Q15_ONE,Pframe), st->echo_suppress), MULT16_16(Pframe, st->echo_suppress_active)),15)); + + compute_gain_floor(st->noise_suppress, effective_echo_suppress, st->noise+N, st->echo_noise+N, st->gain_floor+N, M); + + /* Compute Ephraim & Malah gain speech probability of presence for each critical band (Bark scale) + Technically this is actually wrong because the EM gaim assumes a slightly different probability + distribution */ + for (i=N;iold_ps, echo); - } else { - for (i=1;iupdate_prob[i]<.5f/* || st->ps[i] < st->noise[i]*/) - { - if (echo) - st->noise[i] = .95f*st->noise[i] + .05f*max(1.0f,st->ps[i]-st->frame_size*st->frame_size*4.0*echo[i]); - else - st->noise[i] = .95f*st->noise[i] + .05f*st->ps[i]; - } - } - } + /* See EM and Cohen papers*/ + spx_word32_t theta; + /* Gain from hypergeometric function */ + spx_word32_t MM; + /* Weiner filter gain */ + spx_word16_t prior_ratio; + /* a priority probability of speech presence based on Bark sub-band alone */ + spx_word16_t P1; + /* Speech absence a priori probability (considering sub-band and frame) */ + spx_word16_t q; +#ifdef FIXED_POINT + spx_word16_t tmp; +#endif + + prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT))); + theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT)); - for (i=1;izeta[i] = .7f*st->zeta[i] + .3f*st->prior[i]; + MM = hypergeom_gain(theta); + /* Gain with bound */ + st->gain[i] = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM))); + /* Save old Bark power spectrum */ + st->old_ps[i] = MULT16_32_P15(QCONST16(.2f,15),st->old_ps[i]) + MULT16_32_P15(MULT16_16_P15(QCONST16(.8f,15),SQR16_Q15(st->gain[i])),ps[i]); + + P1 = QCONST16(.199f,15)+MULT16_16_Q15(QCONST16(.8f,15),qcurve (st->zeta[i])); + q = Q15_ONE-MULT16_16_Q15(Pframe,P1); +#ifdef FIXED_POINT + theta = MIN32(theta, EXTEND32(32767)); +/*Q8*/tmp = MULT16_16_Q15((SHL32(1,SNR_SHIFT)+st->prior[i]),EXTRACT16(MIN32(Q15ONE,SHR32(spx_exp(-EXTRACT16(theta)),1)))); + tmp = MIN16(QCONST16(3.,SNR_SHIFT), tmp); /* Prevent overflows in the next line*/ +/*Q8*/tmp = PSHR(MULT16_16(PDIV32_16(SHL32(EXTEND32(q),8),(Q15_ONE-q)),tmp),8); + st->gain2[i]=DIV32_16(SHL(EXTEND32(32767),SNR_SHIFT), ADD16(256,tmp)); +#else + st->gain2[i]=1/(1.f + (q/(1.f-q))*(1+st->prior[i])*exp(-theta)); +#endif } - + /* Convert the EM gains and speech prob to linear frequency */ + filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2); + filterbank_compute_psd16(st->bank,st->gain+N, st->gain); + + /* Use 1 for linear gain resolution (best) or 0 for Bark gain resolution (faster) */ + if (1) { - int freq_start = (int)(300.0f*2.f*N/st->sampling_rate); - int freq_end = (int)(2000.0f*2.f*N/st->sampling_rate); - for (i=freq_start;ibank,st->gain_floor+N, st->gain_floor); + + /* Compute gain according to the Ephraim-Malah algorithm -- linear frequency */ + for (i=0;izeta[i]; - } - Zframe /= (freq_end-freq_start); - } - st->Zlast = Zframe; - - Pframe = qcurve(Zframe); + spx_word32_t MM; + spx_word32_t theta; + spx_word16_t prior_ratio; + spx_word16_t tmp; + spx_word16_t p; + spx_word16_t g; + + /* Wiener filter gain */ + prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT))); + theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT)); + + /* Optimal estimator for loudness domain */ + MM = hypergeom_gain(theta); + /* EM gain with bound */ + g = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM))); + /* Interpolated speech probability of presence */ + p = st->gain2[i]; + + /* Constrain the gain to be close to the Bark scale gain */ + if (MULT16_16_Q15(QCONST16(.333f,15),g) > st->gain[i]) + g = MULT16_16(3,st->gain[i]); + st->gain[i] = g; + + /* Save old power spectrum */ + st->old_ps[i] = MULT16_32_P15(QCONST16(.2f,15),st->old_ps[i]) + MULT16_32_P15(MULT16_16_P15(QCONST16(.8f,15),SQR16_Q15(st->gain[i])),ps[i]); + + /* Apply gain floor */ + if (st->gain[i] < st->gain_floor[i]) + st->gain[i] = st->gain_floor[i]; - /*fprintf (stderr, "%f\n", Pframe);*/ - /* Compute gain according to the Ephraim-Malah algorithm */ - for (i=1;iprior[i]/(1.0001f+st->prior[i]); - theta = (1.f+st->post[i])*prior_ratio; - - if (i==1 || i==N-1) - zeta1 = st->zeta[i]; - else - zeta1 = .25f*st->zeta[i-1] + .5f*st->zeta[i] + .25f*st->zeta[i+1]; - P1 = qcurve (zeta1); - - /* FIXME: add global prob (P2) */ - q = 1-Pframe*P1; - q = 1-P1; - if (q>.95f) - q=.95f; - p=1.f/(1.f + (q/(1.f-q))*(1.f+st->prior[i])*exp(-theta)); - /*p=1;*/ - - /* Optimal estimator for loudness domain */ - MM = hypergeom_gain(theta); + /* Exponential decay model for reverberation (unused) */ + /*st->reverb_estimate[i] = st->reverb_decay*st->reverb_estimate[i] + st->reverb_decay*st->reverb_level*st->gain[i]*st->gain[i]*st->ps[i];*/ + + /* Take into account speech probability of presence (loudness domain MMSE estimator) */ + /* gain2 = [p*sqrt(gain)+(1-p)*sqrt(gain _floor) ]^2 */ + tmp = MULT16_16_P15(p,spx_sqrt(SHL32(EXTEND32(st->gain[i]),15))) + MULT16_16_P15(SUB16(Q15_ONE,p),spx_sqrt(SHL32(EXTEND32(st->gain_floor[i]),15))); + st->gain2[i]=SQR16_Q15(tmp); - st->gain[i] = prior_ratio * MM; - /*Put some (very arbitraty) limit on the gain*/ - if (st->gain[i]>2.f) - { - st->gain[i]=2.f; + /* Use this if you want a log-domain MMSE estimator instead */ + /*st->gain2[i] = pow(st->gain[i], p) * pow(st->gain_floor[i],1.f-p);*/ } - - st->reverb_estimate[i] = st->reverb_decay*st->reverb_estimate[i] + st->reverb_decay*st->reverb_level*st->gain[i]*st->gain[i]*st->ps[i]; - if (st->denoise_enabled) + } else { + for (i=N;igain2[i] = p*p*st->gain[i]; - /*st->gain2[i]=(p*sqrt(st->gain[i])+.05*(1-p))*(p*sqrt(st->gain[i])+.05*(1-p));*/ - /*st->gain2[i] = pow(st->gain[i], p) * pow(.2f,1.f-p);*/ - } else { - st->gain2[i]=1.f; + spx_word16_t tmp; + spx_word16_t p = st->gain2[i]; + st->gain[i] = MAX16(st->gain[i], st->gain_floor[i]); + tmp = MULT16_16_P15(p,spx_sqrt(SHL32(EXTEND32(st->gain[i]),15))) + MULT16_16_P15(SUB16(Q15_ONE,p),spx_sqrt(SHL32(EXTEND32(st->gain_floor[i]),15))); + st->gain2[i]=SQR16_Q15(tmp); } + filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2); } - st->gain2[0]=st->gain[0]=0.f; - st->gain2[N-1]=st->gain[N-1]=0.f; - /* - for (i=30;idenoise_enabled) { - st->gain[i] = st->gain2[i]*st->gain2[i] + (1-st->gain2[i])*.333*(.6*st->gain2[i-1]+st->gain2[i]+.6*st->gain2[i+1]+.4*st->gain2[i-2]+.4*st->gain2[i+2]); + for (i=0;igain2[i]=Q15_ONE; } - for (i=30;igain2[i] = st->gain[i]; - */ + + /*FIXME: This *will* not work for fixed-point */ +#ifndef FIXED_POINT if (st->agc_enabled) - speex_compute_agc(st, mean_prior); - -#if 0 - if (!is_speech) - { - for (i=0;igain2[i] = 0; - } -#if 0 - else { - for (i=0;igain2[i] = 1; - } -#endif + speex_compute_agc(st); #endif - + /* Apply computed gain */ for (i=1;iframe[2*i-1] *= st->gain2[i]; - st->frame[2*i] *= st->gain2[i]; + st->ft[2*i-1] = MULT16_16_P15(st->gain2[i],st->ft[2*i-1]); + st->ft[2*i] = MULT16_16_P15(st->gain2[i],st->ft[2*i]); } - - /* Get rid of the DC and very low frequencies */ - st->frame[0]=0; - st->frame[1]=0; - st->frame[2]=0; - /* Nyquist frequency is mostly useless too */ - st->frame[2*N-1]=0; - + st->ft[0] = MULT16_16_P15(st->gain2[0],st->ft[0]); + st->ft[2*N-1] = MULT16_16_P15(st->gain2[N-1],st->ft[2*N-1]); + /* Inverse FFT with 1/N scaling */ - spx_drft_backward(st->fft_lookup, st->frame); - + spx_ifft(st->fft_lookup, st->ft, st->frame); + /* Scale back to original (lower) amplitude */ for (i=0;i<2*N;i++) - st->frame[i] *= scale; + st->frame[i] = PSHR16(st->frame[i], st->frame_shift); + /*FIXME: This *will* not work for fixed-point */ +#ifndef FIXED_POINT + if (st->agc_enabled) { float max_sample=0; for (i=0;i<2*N;i++) @@ -880,9 +964,11 @@ int speex_preprocess(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo st->frame[i] *= damp; } } - +#endif + + /* Synthesis window (for WOLA) */ for (i=0;i<2*N;i++) - st->frame[i] *= st->window[i]; + st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]); /* Perform overlap and add */ for (i=0;ioutbuf[i] = st->frame[st->frame_size+i]; - /* Save old power spectrum */ - for (i=1;iold_ps[i] = ps[i]; - - return is_speech; + /* FIXME: This VAD is a kludge */ + if (st->vad_enabled) + { + if (Pframe > st->speech_prob_start || (st->was_speech && Pframe > st->speech_prob_continue)) + { + st->was_speech=1; + return 1; + } else + { + st->was_speech=0; + return 0; + } + } else { + return 1; + } } -void speex_preprocess_estimate_update(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo) +void speex_preprocess_estimate_update(SpeexPreprocessState *st, spx_int16_t *x) { int i; int N = st->ps_size; int N3 = 2*N - st->frame_size; + int M; + spx_word32_t *ps=st->ps; - float *ps=st->ps; - + M = st->nbands; + st->min_count++; + preprocess_analysis(st, x); update_noise_prob(st); - - st->nb_preprocess++; for (i=1;iupdate_prob[i]<.5f || st->ps[i] < st->noise[i]) + if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i],NOISE_SHIFT)) { - if (echo) - st->noise[i] = .95f*st->noise[i] + .1f*max(1.0f,st->ps[i]-st->frame_size*st->frame_size*4.0*echo[i]); - else - st->noise[i] = .95f*st->noise[i] + .1f*st->ps[i]; + st->noise[i] = MULT16_32_Q15(QCONST16(.95f,15),st->noise[i]) + MULT16_32_Q15(QCONST16(.05f,15),SHL32(st->ps[i],NOISE_SHIFT)); } } for (i=0;ioutbuf[i] = x[st->frame_size-N3+i]*st->window[st->frame_size+i]; + st->outbuf[i] = MULT16_16_Q15(x[st->frame_size-N3+i],st->window[st->frame_size+i]); /* Save old power spectrum */ - for (i=1;iold_ps[i] = ps[i]; - for (i=1;ireverb_estimate[i] *= st->reverb_decay; + for (i=0;ireverb_estimate[i] = MULT16_32_Q15(st->reverb_decay, st->reverb_estimate[i]); } @@ -946,17 +1040,17 @@ int speex_preprocess_ctl(SpeexPreprocessState *state, int request, void *ptr) switch(request) { case SPEEX_PREPROCESS_SET_DENOISE: - st->denoise_enabled = (*(int*)ptr); + st->denoise_enabled = (*(spx_int32_t*)ptr); break; case SPEEX_PREPROCESS_GET_DENOISE: - (*(int*)ptr) = st->denoise_enabled; + (*(spx_int32_t*)ptr) = st->denoise_enabled; break; - +#ifndef FIXED_POINT case SPEEX_PREPROCESS_SET_AGC: - st->agc_enabled = (*(int*)ptr); + st->agc_enabled = (*(spx_int32_t*)ptr); break; case SPEEX_PREPROCESS_GET_AGC: - (*(int*)ptr) = st->agc_enabled; + (*(spx_int32_t*)ptr) = st->agc_enabled; break; case SPEEX_PREPROCESS_SET_AGC_LEVEL: @@ -969,21 +1063,22 @@ int speex_preprocess_ctl(SpeexPreprocessState *state, int request, void *ptr) case SPEEX_PREPROCESS_GET_AGC_LEVEL: (*(float*)ptr) = st->agc_level; break; - +#endif case SPEEX_PREPROCESS_SET_VAD: - st->vad_enabled = (*(int*)ptr); + speex_warning("The VAD has been replaced by a hack pending a complete rewrite"); + st->vad_enabled = (*(spx_int32_t*)ptr); break; case SPEEX_PREPROCESS_GET_VAD: - (*(int*)ptr) = st->vad_enabled; + (*(spx_int32_t*)ptr) = st->vad_enabled; break; case SPEEX_PREPROCESS_SET_DEREVERB: - st->dereverb_enabled = (*(int*)ptr); + st->dereverb_enabled = (*(spx_int32_t*)ptr); for (i=0;ips_size;i++) st->reverb_estimate[i]=0; break; case SPEEX_PREPROCESS_GET_DEREVERB: - (*(int*)ptr) = st->dereverb_enabled; + (*(spx_int32_t*)ptr) = st->dereverb_enabled; break; case SPEEX_PREPROCESS_SET_DEREVERB_LEVEL: @@ -1001,24 +1096,47 @@ int speex_preprocess_ctl(SpeexPreprocessState *state, int request, void *ptr) break; case SPEEX_PREPROCESS_SET_PROB_START: - st->speech_prob_start = (*(int*)ptr) / 100.0; - if ( st->speech_prob_start > 1 || st->speech_prob_start < 0 ) - st->speech_prob_start = SPEEX_PROB_START_DEFAULT; + *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr)); + st->speech_prob_start = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100); break; case SPEEX_PREPROCESS_GET_PROB_START: - (*(int*)ptr) = st->speech_prob_start * 100; + (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_start, 100); break; case SPEEX_PREPROCESS_SET_PROB_CONTINUE: - st->speech_prob_continue = (*(int*)ptr) / 100.0; - if ( st->speech_prob_continue > 1 || st->speech_prob_continue < 0 ) - st->speech_prob_continue = SPEEX_PROB_CONTINUE_DEFAULT; + *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr)); + st->speech_prob_continue = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100); break; case SPEEX_PREPROCESS_GET_PROB_CONTINUE: - (*(int*)ptr) = st->speech_prob_continue * 100; + (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_continue, 100); + break; + + case SPEEX_PREPROCESS_SET_NOISE_SUPPRESS: + st->noise_suppress = -ABS(*(spx_int32_t*)ptr); + break; + case SPEEX_PREPROCESS_GET_NOISE_SUPPRESS: + (*(spx_int32_t*)ptr) = st->noise_suppress; + break; + case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS: + st->echo_suppress = -ABS(*(spx_int32_t*)ptr); + break; + case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS: + (*(spx_int32_t*)ptr) = st->echo_suppress; + break; + case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS_ACTIVE: + st->echo_suppress_active = -ABS(*(spx_int32_t*)ptr); + break; + case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS_ACTIVE: + (*(spx_int32_t*)ptr) = st->echo_suppress_active; + break; + case SPEEX_PREPROCESS_SET_ECHO_STATE: + st->echo_state = (SpeexEchoState*)ptr; + break; + case SPEEX_PREPROCESS_GET_ECHO_STATE: + ptr = (void*)st->echo_state; break; - default: + default: speex_warning_int("Unknown speex_preprocess_ctl request: ", request); return -1; } diff --git a/pjmedia/src/pjmedia-codec/speex/pseudofloat.h b/pjmedia/src/pjmedia-codec/speex/pseudofloat.h index e85f60e4..67f01b33 100644 --- a/pjmedia/src/pjmedia-codec/speex/pseudofloat.h +++ b/pjmedia/src/pjmedia-codec/speex/pseudofloat.h @@ -36,6 +36,7 @@ #define PSEUDOFLOAT_H #include "misc.h" +#include "math_approx.h" #include #ifdef FIXED_POINT @@ -64,18 +65,8 @@ static inline spx_float_t PSEUDOFLOAT(spx_int32_t x) spx_float_t r = {0,0}; return r; } - while (x>32767) - { - x >>= 1; - /*x *= .5;*/ - e++; - } - while (x<16383) - { - x <<= 1; - /*x *= 2;*/ - e--; - } + e = spx_ilog2(ABS32(x))-14; + x = VSHR32(x, e); if (sign) { spx_float_t r; @@ -166,9 +157,9 @@ static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b) static inline int FLOAT_LT(spx_float_t a, spx_float_t b) { if (a.m==0) - return b.m<0; + return b.m>0; else if (b.m==0) - return a.m>0; + return a.m<0; if ((a).e > (b).e) return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1)); else @@ -204,6 +195,14 @@ static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b) return r; } +static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b) +{ + spx_float_t r; + r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15); + r.e = (a).e+(b).e+15; + return r; +} + static inline spx_float_t FLOAT_SHL(spx_float_t a, int b) { @@ -216,68 +215,53 @@ static inline spx_float_t FLOAT_SHL(spx_float_t a, int b) static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a) { if (a.e<0) - return EXTRACT16((EXTEND32(a.m)+(1<<(-a.e-1)))>>-a.e); + return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e); else return a.m<>-a.e; else - return SHL32(MULT16_32_Q15(a.m, b),15+a.e); + return EXTEND32(a.m)<32767) - { - a >>= 1; - e++; - } - while (a<16384) - { - a <<= 1; - e--; - } - while (b>32767) - { - b >>= 1; - e++; - } - while (b<16384) - { - b <<= 1; - e--; - } + e1 = spx_ilog2(ABS32(a)); + a = VSHR32(a, e1-14); + e2 = spx_ilog2(ABS32(b)); + b = VSHR32(b, e2-14); r.m = MULT16_16_Q15(a,b); - r.e = e+15; + r.e = e1+e2-13; return r; } +/* Do NOT attempt to divide by a negative number */ static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b) { int e=0; spx_float_t r; - /* FIXME: Handle the sign */ if (a==0) { return FLOAT_ZERO; } - while (a=SHL32(EXTEND32(b.m-1),15)) + e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15; + a = VSHR32(a, e); + if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15)) { a >>= 1; e++; @@ -288,41 +272,47 @@ static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b) } +/* Do NOT attempt to divide by a negative number */ static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b) { - int e=0; + int e0=0,e=0; spx_float_t r; - /* FIXME: Handle the sign */ if (a==0) { return FLOAT_ZERO; } - while (b>32767) + if (b>32767) { - b >>= 1; - e--; + e0 = spx_ilog2(b)-14; + b = VSHR32(b, e0); + e0 = -e0; } - while (a=SHL32(b-1,15)) + e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15; + a = VSHR32(a, e); + if (ABS32(a)>=SHL32(EXTEND32(b-1),15)) { a >>= 1; e++; } + e += e0; r.m = DIV32_16(a,b); r.e = e; return r; } +/* Do NOT attempt to divide by a negative number */ static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b) { int e=0; spx_int32_t num; spx_float_t r; + if (b.m<=0) + { + speex_warning_int("Attempted to divide by", b.m); + return FLOAT_ONE; + } num = a.m; + a.m = ABS16(a.m); while (a.m >= b.m) { e++; @@ -334,6 +324,22 @@ static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b) return r; } +static inline spx_float_t FLOAT_SQRT(spx_float_t a) +{ + spx_float_t r; + spx_int32_t m; + m = SHL32(EXTEND32(a.m), 14); + r.e = a.e - 14; + if (r.e & 1) + { + r.e -= 1; + m <<= 1; + } + r.e >>= 1; + r.m = spx_sqrt(m); + return r; +} + #else #define spx_float_t float @@ -342,9 +348,11 @@ static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b) #define FLOAT_HALF 0.5f #define PSEUDOFLOAT(x) (x) #define FLOAT_MULT(a,b) ((a)*(b)) +#define FLOAT_AMULT(a,b) ((a)*(b)) #define FLOAT_MUL32(a,b) ((a)*(b)) #define FLOAT_DIV32(a,b) ((a)/(b)) #define FLOAT_EXTRACT16(a) (a) +#define FLOAT_EXTRACT32(a) (a) #define FLOAT_ADD(a,b) ((a)+(b)) #define FLOAT_SUB(a,b) ((a)-(b)) #define REALFLOAT(x) (x) @@ -354,6 +362,7 @@ static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b) #define FLOAT_LT(a,b) ((a)<(b)) #define FLOAT_GT(a,b) ((a)>(b)) #define FLOAT_DIVU(a,b) ((a)/(b)) +#define FLOAT_SQRT(a) (spx_sqrt(a)) #endif diff --git a/pjmedia/src/pjmedia-codec/speex/quant_lsp_bfin.h b/pjmedia/src/pjmedia-codec/speex/quant_lsp_bfin.h new file mode 100644 index 00000000..c884078e --- /dev/null +++ b/pjmedia/src/pjmedia-codec/speex/quant_lsp_bfin.h @@ -0,0 +1,165 @@ +/* Copyright (C) 2006 David Rowe */ +/** + @file quant_lsp_bfin.h + @author David Rowe + @brief Various compatibility routines for Speex (Blackfin version) +*/ +/* + 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 OVERRIDE_LSP_QUANT +#ifdef OVERRIDE_LSP_QUANT + +/* + Note http://gcc.gnu.org/onlinedocs/gcc/Machine-Constraints.html + well tell you all the magic resgister constraints used below + for gcc in-line asm. +*/ + +static int lsp_quant( + spx_word16_t *x, + const signed char *cdbk, + int nbVec, + int nbDim +) +{ + int j; + spx_word32_t best_dist=1<<30; + int best_id=0; + + __asm__ __volatile__ + ( +" %0 = 1 (X);\n\t" /* %0: best_dist */ +" %0 <<= 30;\n\t" +" %1 = 0 (X);\n\t" /* %1: best_i */ +" P2 = %3\n\t" /* P2: ptr to cdbk */ +" R5 = 0;\n\t" /* R5: best cb entry */ + +" R0 = %5;\n\t" /* set up circ addr */ +" R0 <<= 1;\n\t" +" L0 = R0;\n\t" +" I0 = %2;\n\t" /* %2: &x[0] */ +" B0 = %2;\n\t" + +" R2.L = W [I0++];\n\t" +" LSETUP (lq1, lq2) LC0 = %4;\n\t" +"lq1: R3 = 0;\n\t" /* R3: dist */ +" LSETUP (lq3, lq4) LC1 = %5;\n\t" +"lq3: R1 = B [P2++] (X);\n\t" +" R1 <<= 5;\n\t" +" R0.L = R2.L - R1.L || R2.L = W [I0++];\n\t" +" R0 = R0.L*R0.L;\n\t" +"lq4: R3 = R3 + R0;\n\t" + +" cc =R3<%0;\n\t" +" if cc %0=R3;\n\t" +" if cc %1=R5;\n\t" +"lq2: R5 += 1;\n\t" +" L0 = 0;\n\t" + : "=&d" (best_dist), "=&d" (best_id) + : "a" (x), "b" (cdbk), "a" (nbVec), "a" (nbDim) + : "I0", "P2", "R0", "R1", "R2", "R3", "R5", "L0", "B0", "A0" + ); + + for (j=0;j>> 16;\n\t" +" R1 = (A1 += R2.L*R0.H) (IS);\n\t" +"lwq4: R3 = R3 + R1;\n\t" + +" cc =R3<%0;\n\t" +" if cc %0=R3;\n\t" +" if cc %1=R5;\n\t" +"lwq2: R5 += 1;\n\t" +" L0 = 0;\n\t" +" L1 = 0;\n\t" + : "=&d" (best_dist), "=&d" (best_id) + : "a" (x), "a" (weight), "b" (cdbk), "a" (nbVec), "a" (nbDim) + : "I0", "I1", "P2", "R0", "R1", "R2", "R3", "R5", "A1", + "L0", "L1", "B0", "B1" + ); + + for (j=0;jsubmodes=mode->submodes; st->submodeSelect = st->submodeID=mode->defaultSubmode; - i=9; - speex_encoder_ctl(st->st_low, SPEEX_SET_QUALITY, &i); + tmp=9; + speex_encoder_ctl(st->st_low, SPEEX_SET_QUALITY, &tmp); + tmp=1; + speex_encoder_ctl(st->st_low, SPEEX_SET_WIDEBAND, &tmp); st->lag_factor = mode->lag_factor; st->lpc_floor = mode->lpc_floor; @@ -263,48 +266,48 @@ void *sb_encoder_init(const SpeexMode *m) st->gamma2=mode->gamma2; st->first=1; - st->x0d=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); - st->x1d=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); - st->high=speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); - st->y0=speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); - st->y1=speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); + st->x0d=(spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->x1d=(spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->high=(spx_sig_t*)speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); + st->y0=(spx_sig_t*)speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); + st->y1=(spx_sig_t*)speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); - st->h0_mem=speex_alloc((QMF_ORDER)*sizeof(spx_word16_t)); - st->h1_mem=speex_alloc((QMF_ORDER)*sizeof(spx_word16_t)); - st->g0_mem=speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); - st->g1_mem=speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); + st->h0_mem=(spx_word16_t*)speex_alloc((QMF_ORDER)*sizeof(spx_word16_t)); + st->h1_mem=(spx_word16_t*)speex_alloc((QMF_ORDER)*sizeof(spx_word16_t)); + st->g0_mem=(spx_word32_t*)speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); + st->g1_mem=(spx_word32_t*)speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); - st->excBuf=speex_alloc((st->bufSize)*sizeof(spx_sig_t)); + st->excBuf=(spx_sig_t*)speex_alloc((st->bufSize)*sizeof(spx_sig_t)); st->exc = st->excBuf + st->bufSize - st->windowSize; - st->res=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); - st->sw=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->res=(spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->sw=(spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); st->window= lpc_window; - st->lagWindow = speex_alloc((st->lpcSize+1)*sizeof(spx_word16_t)); + st->lagWindow = (spx_word16_t*)speex_alloc((st->lpcSize+1)*sizeof(spx_word16_t)); for (i=0;ilpcSize+1;i++) st->lagWindow[i]=16384*exp(-.5*sqr(2*M_PI*st->lag_factor*i)); - st->autocorr = speex_alloc((st->lpcSize+1)*sizeof(spx_word16_t)); - st->lpc = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->bw_lpc1 = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->bw_lpc2 = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->lsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->qlsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->old_lsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->old_qlsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->interp_lsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->interp_qlsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->interp_lpc = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->interp_qlpc = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->pi_gain = speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); - st->low_innov = speex_alloc((st->frame_size)*sizeof(spx_word32_t)); + st->autocorr = (spx_word16_t*)speex_alloc((st->lpcSize+1)*sizeof(spx_word16_t)); + st->lpc = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->bw_lpc1 = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->bw_lpc2 = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->lsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->qlsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->old_lsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->old_qlsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->interp_lsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->interp_qlsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->interp_lpc = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->interp_qlpc = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->pi_gain = (spx_word32_t*)speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); + st->low_innov = (spx_word32_t*)speex_alloc((st->frame_size)*sizeof(spx_word32_t)); speex_encoder_ctl(st->st_low, SPEEX_SET_INNOVATION_SAVE, st->low_innov); st->innov_save = NULL; - st->mem_sp = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->mem_sp2 = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); - st->mem_sw = speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_sp = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_sp2 = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); + st->mem_sw = (spx_mem_t*)speex_alloc((st->lpcSize)*sizeof(spx_mem_t)); st->vbr_quality = 8; st->vbr_enabled = 0; @@ -384,8 +387,8 @@ int sb_encode(void *state, void *vin, SpeexBits *bits) VARDECL(spx_word32_t *low_pi_gain); VARDECL(spx_word16_t *low_exc); const SpeexSBMode *mode; - int dtx; - spx_word16_t *in = vin; + spx_int32_t dtx; + spx_word16_t *in = (spx_word16_t*)vin; st = (SBEncState*)state; stack=st->stack; @@ -499,7 +502,7 @@ int sb_encode(void *state, void *vin, SpeexBits *bits) /*if (ratio>-2)*/ if (st->vbr_enabled) { - int modeid; + spx_int32_t modeid; modeid = mode->nb_modes-1; st->relative_quality+=1.0*(ratio+2); if (st->relative_quality<-1) @@ -521,7 +524,7 @@ int sb_encode(void *state, void *vin, SpeexBits *bits) speex_encoder_ctl(state, SPEEX_SET_HIGH_MODE, &modeid); if (st->abr_enabled) { - int bitrate; + spx_int32_t bitrate; speex_encoder_ctl(state, SPEEX_GET_BITRATE, &bitrate); st->abr_drift+=(bitrate-st->abr_enabled); st->abr_drift2 = .95*st->abr_drift2 + .05*(bitrate-st->abr_enabled); @@ -765,7 +768,7 @@ int sb_encode(void *state, void *vin, SpeexBits *bits) /*print_vec(target, st->subframeSize, "\ntarget");*/ SUBMODE(innovation_quant)(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2, SUBMODE(innovation_params), st->lpcSize, st->subframeSize, - innov, syn_resp, bits, stack, (st->complexity+1)>>1, SUBMODE(double_codebook)); + innov, syn_resp, bits, stack, st->complexity, SUBMODE(double_codebook)); /*print_vec(target, st->subframeSize, "after");*/ signal_mul(innov, innov, scale, st->subframeSize); @@ -789,7 +792,7 @@ int sb_encode(void *state, void *vin, SpeexBits *bits) target[i]*=2.5; SUBMODE(innovation_quant)(target, st->interp_qlpc, st->bw_lpc1, st->bw_lpc2, SUBMODE(innovation_params), st->lpcSize, st->subframeSize, - innov2, syn_resp, bits, stack, (st->complexity+1)>>1, 0); + innov2, syn_resp, bits, stack, st->complexity, 0); for (i=0;isubframeSize;i++) innov2[i]*=scale*(1/2.5)/SIG_SCALING; for (i=0;isubframeSize;i++) @@ -834,6 +837,7 @@ int sb_encode(void *state, void *vin, SpeexBits *bits) void *sb_decoder_init(const SpeexMode *m) { + spx_int32_t tmp; SBDecState *st; const SpeexSBMode *mode; st = (SBDecState*)speex_alloc(sizeof(SBDecState)); @@ -860,6 +864,8 @@ void *sb_decoder_init(const SpeexMode *m) st->lpcSize=mode->lpcSize; speex_decoder_ctl(st->st_low, SPEEX_GET_SAMPLING_RATE, &st->sampling_rate); st->sampling_rate*=2; + tmp=1; + speex_decoder_ctl(st->st_low, SPEEX_SET_WIDEBAND, &tmp); st->submodes=mode->submodes; st->submodeID=mode->defaultSubmode; @@ -867,27 +873,27 @@ void *sb_decoder_init(const SpeexMode *m) st->first=1; - st->x0d=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); - st->x1d=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); - st->high=speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); - st->y0=speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); - st->y1=speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); + st->x0d = (spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->x1d = (spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->high = (spx_sig_t*)speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); + st->y0 = (spx_sig_t*)speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); + st->y1 = (spx_sig_t*)speex_alloc((st->full_frame_size)*sizeof(spx_sig_t)); - st->g0_mem=speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); - st->g1_mem=speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); + st->g0_mem = (spx_word32_t*)speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); + st->g1_mem = (spx_word32_t*)speex_alloc((QMF_ORDER)*sizeof(spx_word32_t)); - st->exc=speex_alloc((st->frame_size)*sizeof(spx_sig_t)); - st->excBuf=speex_alloc((st->subframeSize)*sizeof(spx_sig_t)); + st->exc = (spx_sig_t*)speex_alloc((st->frame_size)*sizeof(spx_sig_t)); + st->excBuf = (spx_sig_t*)speex_alloc((st->subframeSize)*sizeof(spx_sig_t)); - st->qlsp = speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); - st->old_qlsp = speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); - st->interp_qlsp = speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); - st->interp_qlpc = speex_alloc(st->lpcSize*sizeof(spx_coef_t)); + st->qlsp = (spx_lsp_t*)speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); + st->old_qlsp = (spx_lsp_t*)speex_alloc((st->lpcSize)*sizeof(spx_lsp_t)); + st->interp_qlsp = (spx_lsp_t*)speex_alloc(st->lpcSize*sizeof(spx_lsp_t)); + st->interp_qlpc = (spx_coef_t*)speex_alloc(st->lpcSize*sizeof(spx_coef_t)); - st->pi_gain = speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); - st->mem_sp = speex_alloc((2*st->lpcSize)*sizeof(spx_mem_t)); + st->pi_gain = (spx_word32_t*)speex_alloc((st->nbSubframes)*sizeof(spx_word32_t)); + st->mem_sp = (spx_mem_t*)speex_alloc((2*st->lpcSize)*sizeof(spx_mem_t)); - st->low_innov = speex_alloc((st->frame_size)*sizeof(spx_word32_t)); + st->low_innov = (spx_word32_t*)speex_alloc((st->frame_size)*sizeof(spx_word32_t)); speex_decoder_ctl(st->st_low, SPEEX_SET_INNOVATION_SAVE, st->low_innov); st->innov_save = NULL; @@ -924,8 +930,8 @@ void sb_decoder_destroy(void *state) speex_free(st->interp_qlsp); speex_free(st->interp_qlpc); speex_free(st->pi_gain); - speex_free(st->mem_sp); speex_free(st->low_innov); + speex_free(st->mem_sp); speex_free(state); } @@ -986,9 +992,9 @@ int sb_decode(void *state, SpeexBits *bits, void *vout) VARDECL(spx_word32_t *low_pi_gain); VARDECL(spx_word16_t *low_exc); VARDECL(spx_coef_t *ak); - int dtx; + spx_int32_t dtx; const SpeexSBMode *mode; - spx_word16_t *out = vout; + spx_word16_t *out = (spx_word16_t*)vout; st = (SBDecState*)state; stack=st->stack; @@ -1192,7 +1198,7 @@ int sb_decode(void *state, SpeexBits *bits, void *vout) scale = SHL(MULT16_16(PDIV32_16(SHL(gc,SIG_SHIFT-6),filter_ratio),(1+el)),6); SUBMODE(innovation_unquant)(exc, SUBMODE(innovation_params), st->subframeSize, - bits, stack); + bits, stack, &st->seed); signal_mul(exc,exc,scale,st->subframeSize); @@ -1203,7 +1209,7 @@ int sb_decode(void *state, SpeexBits *bits, void *vout) for (i=0;isubframeSize;i++) innov2[i]=0; SUBMODE(innovation_unquant)(innov2, SUBMODE(innovation_params), st->subframeSize, - bits, stack); + bits, stack, &st->seed); for (i=0;isubframeSize;i++) innov2[i]*=scale/(float)SIG_SCALING*(1/2.5); for (i=0;isubframeSize;i++) @@ -1251,10 +1257,10 @@ int sb_encoder_ctl(void *state, int request, void *ptr) switch(request) { case SPEEX_GET_FRAME_SIZE: - (*(int*)ptr) = st->full_frame_size; + (*(spx_int32_t*)ptr) = st->full_frame_size; break; case SPEEX_SET_HIGH_MODE: - st->submodeSelect = st->submodeID = (*(int*)ptr); + st->submodeSelect = st->submodeID = (*(spx_int32_t*)ptr); break; case SPEEX_SET_LOW_MODE: speex_encoder_ctl(st->st_low, SPEEX_SET_LOW_MODE, ptr); @@ -1272,22 +1278,22 @@ int sb_encoder_ctl(void *state, int request, void *ptr) speex_encoder_ctl(st, SPEEX_SET_QUALITY, ptr); break; case SPEEX_SET_VBR: - st->vbr_enabled = (*(int*)ptr); + st->vbr_enabled = (*(spx_int32_t*)ptr); speex_encoder_ctl(st->st_low, SPEEX_SET_VBR, ptr); break; case SPEEX_GET_VBR: - (*(int*)ptr) = st->vbr_enabled; + (*(spx_int32_t*)ptr) = st->vbr_enabled; break; case SPEEX_SET_VAD: - st->vad_enabled = (*(int*)ptr); + st->vad_enabled = (*(spx_int32_t*)ptr); speex_encoder_ctl(st->st_low, SPEEX_SET_VAD, ptr); break; case SPEEX_GET_VAD: - (*(int*)ptr) = st->vad_enabled; + (*(spx_int32_t*)ptr) = st->vad_enabled; break; case SPEEX_SET_VBR_QUALITY: { - int q; + spx_int32_t q; float qual = (*(float*)ptr)+.6; st->vbr_quality = (*(float*)ptr); if (qual>10) @@ -1308,7 +1314,7 @@ int sb_encoder_ctl(void *state, int request, void *ptr) speex_encoder_ctl(st->st_low, SPEEX_SET_VBR, &st->vbr_enabled); if (st->vbr_enabled) { - int i=10, rate, target; + spx_int32_t i=10, rate, target; float vbr_qual; target = (*(spx_int32_t*)ptr); while (i>=0) @@ -1334,8 +1340,8 @@ int sb_encoder_ctl(void *state, int request, void *ptr) break; case SPEEX_SET_QUALITY: { - int nb_qual; - int quality = (*(int*)ptr); + spx_int32_t nb_qual; + int quality = (*(spx_int32_t*)ptr); if (quality < 0) quality = 0; if (quality > 10) @@ -1347,16 +1353,16 @@ int sb_encoder_ctl(void *state, int request, void *ptr) break; case SPEEX_SET_COMPLEXITY: speex_encoder_ctl(st->st_low, SPEEX_SET_COMPLEXITY, ptr); - st->complexity = (*(int*)ptr); + st->complexity = (*(spx_int32_t*)ptr); if (st->complexity<1) st->complexity=1; break; case SPEEX_GET_COMPLEXITY: - (*(int*)ptr) = st->complexity; + (*(spx_int32_t*)ptr) = st->complexity; break; case SPEEX_SET_BITRATE: { - int i=10; + spx_int32_t i=10; spx_int32_t rate, target; target = (*(spx_int32_t*)ptr); while (i>=0) @@ -1404,15 +1410,15 @@ int sb_encoder_ctl(void *state, int request, void *ptr) } break; case SPEEX_SET_SUBMODE_ENCODING: - st->encode_submode = (*(int*)ptr); + st->encode_submode = (*(spx_int32_t*)ptr); speex_encoder_ctl(st->st_low, SPEEX_SET_SUBMODE_ENCODING, &ptr); break; case SPEEX_GET_SUBMODE_ENCODING: - (*(int*)ptr) = st->encode_submode; + (*(spx_int32_t*)ptr) = st->encode_submode; break; case SPEEX_GET_LOOKAHEAD: speex_encoder_ctl(st->st_low, SPEEX_GET_LOOKAHEAD, ptr); - (*(int*)ptr) = 2*(*(int*)ptr) + QMF_ORDER - 1; + (*(spx_int32_t*)ptr) = 2*(*(spx_int32_t*)ptr) + QMF_ORDER - 1; break; case SPEEX_SET_PLC_TUNING: speex_encoder_ctl(st->st_low, SPEEX_SET_PLC_TUNING, ptr); @@ -1429,7 +1435,6 @@ int sb_encoder_ctl(void *state, int request, void *ptr) st->vbr_max_high = 17600; } else { spx_int32_t low_rate; - /* FIXME: Need to adapt that to ultra-wideband */ if (st->vbr_max >= 42200) { st->vbr_max_high = 17600; @@ -1442,6 +1447,8 @@ int sb_encoder_ctl(void *state, int request, void *ptr) } else { st->vbr_max_high = 1800; } + if (st->subframeSize==80) + st->vbr_max_high = 1800; low_rate = st->vbr_max - st->vbr_max_high; speex_encoder_ctl(st->st_low, SPEEX_SET_VBR_MAX_BITRATE, &low_rate); } @@ -1450,6 +1457,12 @@ int sb_encoder_ctl(void *state, int request, void *ptr) case SPEEX_GET_VBR_MAX_BITRATE: (*(spx_int32_t*)ptr) = st->vbr_max; break; + case SPEEX_SET_HIGHPASS: + speex_encoder_ctl(st->st_low, SPEEX_SET_HIGHPASS, ptr); + break; + case SPEEX_GET_HIGHPASS: + speex_encoder_ctl(st->st_low, SPEEX_GET_HIGHPASS, ptr); + break; /* This is all internal stuff past this point */ @@ -1485,8 +1498,12 @@ int sb_encoder_ctl(void *state, int request, void *ptr) (*(float*)ptr)=st->relative_quality; break; case SPEEX_SET_INNOVATION_SAVE: - st->innov_save = ptr; + st->innov_save = (spx_sig_t*)ptr; + break; + case SPEEX_SET_WIDEBAND: + speex_encoder_ctl(st->st_low, SPEEX_SET_WIDEBAND, ptr); break; + default: speex_warning_int("Unknown nb_ctl request: ", request); return -1; @@ -1501,7 +1518,7 @@ int sb_decoder_ctl(void *state, int request, void *ptr) switch(request) { case SPEEX_SET_HIGH_MODE: - st->submodeID = (*(int*)ptr); + st->submodeID = (*(spx_int32_t*)ptr); break; case SPEEX_SET_LOW_MODE: speex_decoder_ctl(st->st_low, SPEEX_SET_LOW_MODE, ptr); @@ -1510,20 +1527,20 @@ int sb_decoder_ctl(void *state, int request, void *ptr) speex_decoder_ctl(st->st_low, SPEEX_GET_LOW_MODE, ptr); break; case SPEEX_GET_FRAME_SIZE: - (*(int*)ptr) = st->full_frame_size; + (*(spx_int32_t*)ptr) = st->full_frame_size; break; case SPEEX_SET_ENH: speex_decoder_ctl(st->st_low, request, ptr); - st->lpc_enh_enabled = *((int*)ptr); + st->lpc_enh_enabled = *((spx_int32_t*)ptr); break; case SPEEX_GET_ENH: - *((int*)ptr) = st->lpc_enh_enabled; + *((spx_int32_t*)ptr) = st->lpc_enh_enabled; break; case SPEEX_SET_MODE: case SPEEX_SET_QUALITY: { - int nb_qual; - int quality = (*(int*)ptr); + spx_int32_t nb_qual; + int quality = (*(spx_int32_t*)ptr); if (quality < 0) quality = 0; if (quality > 10) @@ -1567,16 +1584,23 @@ int sb_decoder_ctl(void *state, int request, void *ptr) } break; case SPEEX_SET_SUBMODE_ENCODING: - st->encode_submode = (*(int*)ptr); + st->encode_submode = (*(spx_int32_t*)ptr); speex_decoder_ctl(st->st_low, SPEEX_SET_SUBMODE_ENCODING, &ptr); break; case SPEEX_GET_SUBMODE_ENCODING: - (*(int*)ptr) = st->encode_submode; + (*(spx_int32_t*)ptr) = st->encode_submode; break; case SPEEX_GET_LOOKAHEAD: speex_decoder_ctl(st->st_low, SPEEX_GET_LOOKAHEAD, ptr); - (*(int*)ptr) = 2*(*(int*)ptr); + (*(spx_int32_t*)ptr) = 2*(*(spx_int32_t*)ptr); + break; + case SPEEX_SET_HIGHPASS: + speex_decoder_ctl(st->st_low, SPEEX_SET_HIGHPASS, ptr); + break; + case SPEEX_GET_HIGHPASS: + speex_decoder_ctl(st->st_low, SPEEX_GET_HIGHPASS, ptr); break; + case SPEEX_GET_PI_GAIN: { int i; @@ -1609,8 +1633,12 @@ int sb_decoder_ctl(void *state, int request, void *ptr) speex_decoder_ctl(st->st_low, SPEEX_GET_DTX_STATUS, ptr); break; case SPEEX_SET_INNOVATION_SAVE: - st->innov_save = ptr; + st->innov_save = (spx_sig_t*)ptr; + break; + case SPEEX_SET_WIDEBAND: + speex_decoder_ctl(st->st_low, SPEEX_SET_WIDEBAND, ptr); break; + default: speex_warning_int("Unknown nb_ctl request: ", request); return -1; diff --git a/pjmedia/src/pjmedia-codec/speex/sb_celp.h b/pjmedia/src/pjmedia-codec/speex/sb_celp.h index 194cdd93..4da03e43 100644 --- a/pjmedia/src/pjmedia-codec/speex/sb_celp.h +++ b/pjmedia/src/pjmedia-codec/speex/sb_celp.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2002 Jean-Marc Valin */ +/* Copyright (C) 2002-2006 Jean-Marc Valin */ /** @file sb_celp.h @brief Sub-band CELP mode used for wideband encoding diff --git a/pjmedia/src/pjmedia-codec/speex/speex.c b/pjmedia/src/pjmedia-codec/speex/speex.c index 94829e61..846e0218 100644 --- a/pjmedia/src/pjmedia-codec/speex/speex.c +++ b/pjmedia/src/pjmedia-codec/speex/speex.c @@ -86,7 +86,7 @@ int speex_decode_native(void *state, SpeexBits *bits, spx_word16_t *out) int speex_encode(void *state, float *in, SpeexBits *bits) { int i; - int N; + spx_int32_t N; spx_int16_t short_in[MAX_IN_SAMPLES]; speex_encoder_ctl(state, SPEEX_GET_FRAME_SIZE, &N); for (i=0;idec(state, bits, short_out); @@ -136,7 +136,7 @@ int speex_encode(void *state, float *in, SpeexBits *bits) int speex_encode_int(void *state, spx_int16_t *in, SpeexBits *bits) { int i; - int N; + spx_int32_t N; float float_in[MAX_IN_SAMPLES]; speex_encoder_ctl(state, SPEEX_GET_FRAME_SIZE, &N); for (i=0;i +#include +#include + +/* psychoacoustic setup ********************************************/ + +static VorbisPsyInfo example_tuning = { + + .5,.5, + 3,3,25, + + /*63 125 250 500 1k 2k 4k 8k 16k*/ + // vorbis mode 4 style + //{-32,-32,-32,-32,-28,-24,-22,-20,-20, -20, -20, -8, -6, -6, -6, -6, -6}, + { -4, -6, -6, -6, -6, -6, -6, -6, -8, -8,-10,-10, -8, -6, -4, -4, -2}, + + { + 0, 1, 2, 3, 4, 5, 5, 5, /* 7dB */ + 6, 6, 6, 5, 4, 4, 4, 4, /* 15dB */ + 4, 4, 5, 5, 5, 6, 6, 6, /* 23dB */ + 7, 7, 7, 8, 8, 8, 9, 10, /* 31dB */ + 11,12,13,14,15,16,17, 18, /* 39dB */ + } + +}; + + + +/* there was no great place to put this.... */ +#include +static void _analysis_output(char *base,int i,float *v,int n,int bark,int dB){ + int j; + FILE *of; + char buffer[80]; + + sprintf(buffer,"%s_%d.m",base,i); + of=fopen(buffer,"w"); + + if(!of)perror("failed to open data dump file"); + + for(j=0;j> 16; + if( lo>=0 ) break; + hi = b[i] & 0xffff; + + tN = N[hi] + N[-lo]; + tX = X[hi] - X[-lo]; + tXX = XX[hi] + XX[-lo]; + tY = Y[hi] + Y[-lo]; + tXY = XY[hi] - XY[-lo]; + + A = tY * tXX - tX * tXY; + B = tN * tXY - tX * tY; + D = tN * tXX - tX * tX; + R = (A + x * B) / D; + if (R < 0.f) + R = 0.f; + + noise[i] = R - offset; + } + + for ( ;; i++, x += 1.f) { + + lo = b[i] >> 16; + hi = b[i] & 0xffff; + if(hi>=n)break; + + tN = N[hi] - N[lo]; + tX = X[hi] - X[lo]; + tXX = XX[hi] - XX[lo]; + tY = Y[hi] - Y[lo]; + tXY = XY[hi] - XY[lo]; + + A = tY * tXX - tX * tXY; + B = tN * tXY - tX * tY; + D = tN * tXX - tX * tX; + R = (A + x * B) / D; + if (R < 0.f) R = 0.f; + + noise[i] = R - offset; + } + for ( ; i < n; i++, x += 1.f) { + + R = (A + x * B) / D; + if (R < 0.f) R = 0.f; + + noise[i] = R - offset; + } + + if (fixed <= 0) return; + + for (i = 0, x = 0.f;; i++, x += 1.f) { + hi = i + fixed / 2; + lo = hi - fixed; + if(lo>=0)break; + + tN = N[hi] + N[-lo]; + tX = X[hi] - X[-lo]; + tXX = XX[hi] + XX[-lo]; + tY = Y[hi] + Y[-lo]; + tXY = XY[hi] - XY[-lo]; + + + A = tY * tXX - tX * tXY; + B = tN * tXY - tX * tY; + D = tN * tXX - tX * tX; + R = (A + x * B) / D; + + if (R - offset < noise[i]) noise[i] = R - offset; + } + for ( ;; i++, x += 1.f) { + + hi = i + fixed / 2; + lo = hi - fixed; + if(hi>=n)break; + + tN = N[hi] - N[lo]; + tX = X[hi] - X[lo]; + tXX = XX[hi] - XX[lo]; + tY = Y[hi] - Y[lo]; + tXY = XY[hi] - XY[lo]; + + A = tY * tXX - tX * tXY; + B = tN * tXY - tX * tY; + D = tN * tXX - tX * tX; + R = (A + x * B) / D; + + if (R - offset < noise[i]) noise[i] = R - offset; + } + for ( ; i < n; i++, x += 1.f) { + R = (A + x * B) / D; + if (R - offset < noise[i]) noise[i] = R - offset; + } +} + +static void _vp_noisemask(VorbisPsy *p, + float *logfreq, + float *logmask){ + + int i,n=p->n/2; + float *work=alloca(n*sizeof(*work)); + + bark_noise_hybridmp(n,p->bark,logfreq,logmask, + 140.,-1); + + for(i=0;ibark,work,logmask,0., + p->vi->noisewindowfixed); + + for(i=0;i=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1; + if(dB<0)dB=0; + logmask[i]= work[i]+p->vi->noisecompand[dB]+p->noiseoffset[i]; + } + +} + +VorbisPsy *vorbis_psy_init(int rate, int n) +{ + long i,j,lo=-99,hi=1; + VorbisPsy *p = speex_alloc(sizeof(VorbisPsy)); + memset(p,0,sizeof(*p)); + + p->n = n; + spx_drft_init(&p->lookup, n); + p->bark = speex_alloc(n*sizeof(*p->bark)); + p->rate=rate; + p->vi = &example_tuning; + + /* BH4 window */ + p->window = speex_alloc(sizeof(*p->window)*n); + float a0 = .35875f; + float a1 = .48829f; + float a2 = .14128f; + float a3 = .01168f; + for(i=0;iwindow[i] = //a0 - a1*cos(2.*M_PI/n*(i+.5)) + a2*cos(4.*M_PI/n*(i+.5)) - a3*cos(6.*M_PI/n*(i+.5)); + sin((i+.5)/n * M_PI)*sin((i+.5)/n * M_PI); + /* bark scale lookups */ + for(i=0;ivi->noisewindowlominvi->noisewindowlo);lo++); + + for(;hi<=n && (hivi->noisewindowhimin || + toBARK(rate/(2*n)*hi)<(bark+p->vi->noisewindowhi));hi++); + + p->bark[i]=((lo-1)<<16)+(hi-1); + + } + + /* set up rolling noise median */ + p->noiseoffset=speex_alloc(n*sizeof(*p->noiseoffset)); + + for(i=0;i=P_BANDS-1)halfoc=P_BANDS-1; + inthalfoc=(int)halfoc; + del=halfoc-inthalfoc; + + p->noiseoffset[i]= + p->vi->noiseoff[inthalfoc]*(1.-del) + + p->vi->noiseoff[inthalfoc+1]*del; + + } +#if 0 + _analysis_output_always("noiseoff0",ls,p->noiseoffset,n,1,0,0); +#endif + + return p; +} + +void vorbis_psy_destroy(VorbisPsy *p) +{ + if(p){ + spx_drft_clear(&p->lookup); + if(p->bark) + speex_free(p->bark); + if(p->noiseoffset) + speex_free(p->noiseoffset); + if(p->window) + speex_free(p->window); + memset(p,0,sizeof(*p)); + speex_free(p); + } +} + +void compute_curve(VorbisPsy *psy, float *audio, float *curve) +{ + int i; + float work[psy->n]; + + float scale=4.f/psy->n; + float scale_dB; + + scale_dB=todB(scale); + + /* window the PCM data; use a BH4 window, not vorbis */ + for(i=0;in;i++) + work[i]=audio[i] * psy->window[i]; + + { + static int seq=0; + + //_analysis_output("win",seq,work,psy->n,0,0); + + seq++; + } + + /* FFT yields more accurate tonal estimation (not phase sensitive) */ + spx_drft_forward(&psy->lookup,work); + + /* magnitudes */ + work[0]=scale_dB+todB(work[0]); + for(i=1;in-1;i+=2){ + float temp = work[i]*work[i] + work[i+1]*work[i+1]; + work[(i+1)>>1] = scale_dB+.5f * todB(temp); + } + + /* derive a noise curve */ + _vp_noisemask(psy,work,curve); +#define SIDEL 12 + for (i=0;in>>1)-i-1]=curve[(psy->n>>1)-SIDEH]; + } + for(i=0;i<((psy->n)>>1);i++) + curve[i] = fromdB(1.2*curve[i]+.2*i); + //curve[i] = fromdB(0.8*curve[i]+.35*i); + //curve[i] = fromdB(0.9*curve[i])*pow(1.0*i+45,1.3); +} + +/* Transform a masking curve (power spectrum) into a pole-zero filter */ +void curve_to_lpc(VorbisPsy *psy, float *curve, float *awk1, float *awk2, int ord) +{ + int i; + float ac[psy->n]; + float tmp; + int len = psy->n >> 1; + for (i=0;i<2*len;i++) + ac[i] = 0; + for (i=1;ilookup, ac); + _spx_lpc(awk1, ac, ord); + tmp = 1.; + for (i=0;ilookup, ac); + /* Compute (power) response of awk1 (all zero) */ + ac[0] *= ac[0]; + for (i=1;ilookup, ac); + _spx_lpc(awk2, ac, ord); + tmp = 1; + for (i=0;i +#include + +#define ORDER 10 +#define CURVE_SIZE 24 + +int main() +{ + int i; + float curve[CURVE_SIZE]; + float awk1[ORDER], awk2[ORDER]; + for (i=0;i